trifle

技術メモ

線形計画ソルバー PuLP で強双対性・相補性の確認

線形計画問題はなんとなくエクセルとかで解くイメージがあったけれど, Python にもちゃんとライブラリがあった

github.com

ちょっと使いたくなったので触ってみます
線形計画問題において,

  •  \min c^{T} x \ \ s.t. \ \ A x = b, x \geq 0 (主問題)
  •  \max b^{T} y \ \ s.t. \ \ A^{T} y \leq c(双対問題)

のそれぞれの実行可能な最適解を x_t, y_t とすると, 強双対定理より  c^{T} x_t = b^{T} y_t が, 相補性定理より  (c - A^{T} y_t) \cdot x_t = 0 が成り立つ. これを PuLP で確認しましょう

from pulp import *

def solve_primal_problem(A, b, c, m, n):
    prob = LpProblem("primal problem", LpMinimize)

    x = []
    for i in range(n):
        x.append(LpVariable("x{}".format(i), 0))

    prob += lpDot(c, x) # objective score

    for i in range(m):
        prob += lpDot(A[i], x) == b[i]

    prob.solve()
    return [x[i].varValue for i in range(n)]

def solve_dual_problem(A_T, b, c, m, n):
    prob = LpProblem("dual problem", LpMaximize)

    y = []
    for i in range(m):
        y.append(LpVariable("y{}".format(i)))

    prob += lpDot(b, y) # objective score

    for i in range(n):
        prob += lpDot(A_T[i], y) <= c[i]

    prob.solve()
    return [y[i].varValue for i in range(m)]


# 転置行列を求めるやつ
# まあ numpy にあるメソッドを使えばいいんだけどこのためだけに numpy を import するのもあれなので
def transpose(A, m, n):
    newA = [[0 for i in range(m)] for j in range(n)]
    for i in range(n):
        for j in range(m):
            newA[i][j] = A[j][i]
    return newA

上から順に主問題の最適解を求める関数, 双対問題の最適解を求める関数で, まあコードを見れば使い方がなんとなく分かると思います

# Example

m = 2
n = 4
A = [[2, 1, -1, 0],
     [1, 3, 0, -1]]
A_T = transpose(A, m, n)
b = [3, 4]
c = [1, 2, 0, 0]

x_t = solve_primal_problem(A, b, c, m, n)
print(x_t)
y_t = solve_dual_problem(A_T, b, c, m, n)
print(y_t)

# 強双対性
print(lpDot(c, x_t))
print(lpDot(b, y_t))
# が一致するはず

# 相補性
z_t = [c[i] - lpDot(A_T[i], y_t).value() for i in range(n)]
print(lpDot(x_t, z_t))
# が0になるはず

実行すると順に

[1.0, 1.0, 0.0, 0.0]
[0.2, 0.6]
3.0
3.0
2.220446049250313e-16

で, 強双対性の方はバッチリ2つのスコアが合っていて, 相補性の方は本当はぴったり0になってほしいがまあ浮動小数点なので無理ですよね, という感じ


もう少し複雑な例でやると,

# Example2

m = 3
n = 6
A = [[4, 6, 3, 1, 0, 0],
     [2, 3, 6, 0, 1, 0],
     [3, 1, 0, 0, 0, 1]]
A_T = transpose(A, m, n)
b = [24, 24, 12]
c = [-1, -12, -10, 0, 0, 0]

x_t = solve_primal_problem(A, b, c, m, n)
print(x_t)
y_t = solve_dual_problem(A_T, b, c, m, n)
print(y_t)

# 強双対性
print(lpDot(c, x_t))
print(lpDot(b, y_t))
# が一致するはず

# 相補性
z_t = [c[i] - lpDot(A_T[i], y_t).value() for i in range(n)]
print(lpDot(x_t, z_t))
# が0になるはず

実行すると順に

[0.0, 2.6666667, 2.6666667, 0.0, 0.0, 9.3333333]
[-1.5555556, -0.88888889, 0.0]
-58.666667399999994
-58.666667759999996
1.0933333474607257e-06