😸
単体法入門：Pythonで学ぶ数理最適化アルゴリズム
はじめに
この記事で使用したソースコード、データはこちらを参照ください。
単体法の流れ
今回は不等式標準形等の最大化問題を考えます。
ステップ1：最適性判定
- 目的関数の係数が全て負⇒最適解
- ある係数が非負⇒基底に入る非基底変数
x_s
ステップ2：非有界性判定、ピボット演算
- 変数
x_s - 無限に増やせる⇒非有界
- それ以外⇒
x_s
- 基底変数を基底から出る変数
x_r - 非基底変数と基底変数を入れ替えてステップ1へ
単体法のピボット演算規則
単体法では、ピボット演算の際に基底に入れる非基底変数と基底から出る基底変数の組み合わせを選び、目的関数の値を変化させて最適解を探します。この非基底変数と基底変数を選ぶ規則にはいくつかの種類があります。
この記事では、以下の規則についてPythonで実装し、その解説を行います。
- 最大係数規則
- 最小添え字規則
- 最大改善規則
- 最急降下(Steeppest Edge)規則
最大係数規則
説明
目的関数の係数が最大（絶対値が最大）の非基底変数を選びます。
全体のコード
最大係数規則のコード
simplex_max_rc.py
import numpy as np
error = 1.0e-10 # 許容誤差
def print_simplex_detail(
iter=None,
m=None,
n=None,
basis=None,
nonbasis=None,
x_b=None,
y=None,
rc=None,
s=None,
d=None,
r=None,
):
if iter is not None:
print(f"iter = {iter}")
if m is not None and n is not None:
print(f"(m, n) = {(m, n)}")
if basis is not None:
print(f"basis = {basis}")
if nonbasis is not None:
print(f"nonbasis = {nonbasis}")
if y is not None:
print(f"y = {y}")
if rc is not None:
print(f"rc = {rc}")
if s is not None:
print(f"s = {s}")
if d is not None:
print(f"d = {d}")
if x_b is not None:
print(f"x_basis = {x_b}")
if r is not None:
print(f"r = {r}")
def lp_simplex(A, b, c):
(m, n) = A.shape # m:制約条件数, n:変数数
print("m = {}, n = {}".format(m, n))
# 初期化
Ai = np.hstack((A, np.identity(m)))
c0 = np.r_[c, np.zeros(m)]
basis = [n + i for i in range(m)]
nonbasis = [j for j in range(n)]
iter = 0
while True:
# 巡回の確認の際に使用
#while iter < 10:
# 基本解の計算
x = np.zeros(n + m)
x[basis] = np.linalg.solve(Ai[:, basis], b)
# 双対変数の計算
y = np.linalg.solve(Ai[:, basis].T, c0[basis])
# 現在の目的関数の係数を計算
rc = c0[nonbasis] - y @ Ai[:, nonbasis]
# 最適性のチェック（双対可能性）
if np.all(rc <= error):
print("--optimal solution found---------------------------------------------------")
print_simplex_detail(iter=iter)
print("obj.val. = {}".format(c0[basis] @ x[basis]))
print("x = {}".format(x))
print("---------------------------------------------------------------------------")
break
# 入る変数の決定(最大係数)
s = np.argmax(rc)
d = np.linalg.solve(
Ai[:, basis], Ai[:, nonbasis[s]]
)
# 非有界性のチェック
if np.all(d <= error):
print("problem is unbounded")
break
# 出る変数の決定
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
r = np.argmin(ratio)
nonbasis[s], basis[r] = basis[r], nonbasis[s]
if iter % 50 == 0:
print("iter: {}".format(iter))
print("current obj.val. = {}".format(x[basis] @ c0[basis]))
# 巡回の確認の際に使用
#print_simplex_detail(iter=iter, basis=basis, nonbasis=nonbasis, rc=rc, d=d, s=s, r=r, x_b=x[basis])
iter += 1
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 巡回する例
#A = np.array([[1, 12, -2, -12], [0.25, 1, -0.25, -2], [1, -4, 0, -8]])
#b = np.array([0, 0, 1])
#c = np.array([1, -4, 0, -8])
#lp_simplex(A, b, c)
コードの解説
今回は行列計算をnumpyで解きます。許容誤差は
import numpy as np
error = 1.0e-10 # 許容誤差
単体法の各ステップで詳細をプリントするための関数。
def print_simplex_detail(
iter=None,
m=None,
n=None,
basis=None,
nonbasis=None,
x_b=None,
y=None,
rc=None,
s=None,
d=None,
r=None,
):
if iter is not None:
print(f"iter = {iter}")
if m is not None and n is not None:
print(f"(m, n) = {(m, n)}")
if basis is not None:
print(f"basis = {basis}")
if nonbasis is not None:
print(f"nonbasis = {nonbasis}")
if y is not None:
print(f"y = {y}")
if rc is not None:
print(f"rc = {rc}")
if s is not None:
print(f"s = {s}")
if d is not None:
print(f"d = {d}")
if x_b is not None:
print(f"x_basis = {x_b}")
if r is not None:
print(f"r = {r}")
Aはmは制約条件の数、nは変数の数を表します。
bは制約条件の値です。cは目的関数の係数です。
Aにスラック変数を導入したAiとします。
目的関数関数にもスラック変数を導入し、その係数をc0とします。
basisとnonbasisには基底変数と非基底変数の添え字を保存しときます。
def lp_simplex(A, b, c):
(m, n) = A.shape # m:制約条件数, n:変数数
print("m = {}, n = {}".format(m, n))
# 初期化
Ai = np.hstack((A, np.identity(m)))
c0 = np.r_[c, np.zeros(m)]
basis = [n + i for i in range(m)]
nonbasis = [j for j in range(n)]
iter = 0
ループの中では毎ループごとに基本解(x[basis])、双対変数(y)、目的関数の係数(rc)を計算し、最適性のチェックをします。目的関数の係数が全て負であれば最適解が得られます。
while True:
# 巡回の確認の際に使用
#while iter < 10:
# 基本解の計算
x = np.zeros(n + m)
x[basis] = np.linalg.solve(Ai[:, basis], b)
# 双対変数の計算
y = np.linalg.solve(Ai[:, basis].T, c0[basis])
# 現在の目的関数の係数を計算
rc = c0[nonbasis] - y @ Ai[:, nonbasis]
# 最適性のチェック（双対可能性）
if np.all(rc <= error):
print("--optimal solution found---------------------------------------------------")
print_simplex_detail(iter=iter, basis=basis, nonbasis=nonbasis, rc=rc, x_b=x[basis])
print("obj.val. = {}".format(c0[basis] @ x[basis]))
print("x = {}".format(x))
print("---------------------------------------------------------------------------")
break
最大係数規則では目的関数の係数rcの中から最大の要素を持つインデックスsを選び、対応する列dを計算します。
dの要素がすべて許容誤差以下であれば、問題は非有界であると判断し、ループを終了します。そうでない場合、比率ratioを計算し、最小の比率を持つインデックス
# 入る変数の決定(最大係数)
s = np.argmax(rc)
d = np.linalg.solve(
Ai[:, basis], Ai[:, nonbasis[s]]
)
# 非有界性のチェック
if np.all(d <= error):
print("problem is unbounded")
break
# 出る変数の決定
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
r = np.argmin(ratio)
nonbasis[s], basis[r] = basis[r], nonbasis[s]
if iter % 50 == 0:
print("iter: {}".format(iter))
print("current obj.val. = {}".format(x[basis] @ c0[basis]))
# 巡回の確認の際に使用
#print_simplex_detail(iter=iter, basis=basis, nonbasis=nonbasis, rc=rc, d=d, s=s, r=r, x_b=x[basis])
iter += 1
実行結果
例題として以下のような問題を用意しました。巡回に関しては後述。
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 巡回する例
#A = np.array([[1, 12, -2, -12], [0.25, 1, -0.25, -2], [1, -4, 0, -8]])
#b = np.array([0, 0, 1])
#c = np.array([1, -4, 0, -8])
#lp_simplex(A, b, c)
実行結果は以下のようになりました。目的関数値は16.0、最適化解はx = [2. 1. 3. 0. 0. 0.]です。
最大係数規則の実行結果
m = 3, n = 3
iter: 0
current obj.val. = 0.0
--optimal solution found---------------------------------------------------
iter = 3
basis = [0, 2, 1]
nonbasis = [3, 5, 4]
rc = [-1.33333333 -1.33333333 -0.33333333]
x_basis = [2. 3. 1.]
obj.val. = 16.0
x = [2. 1. 3. 0. 0. 0.]
---------------------------------------------------------------------------
最小添え字規則
説明
単体法では退化が起こると、巡回が発生して、解が求まらない場合があります。巡回が発生しないようにするための規則として、最小添え字規則が知られています。
最小添え字規則では、入れ替え可能な非基底変数および基底変数の中から、元の添え字が最も小さいものを選びます。
退化:
基底解の値に0が含まれる状態です。具体的には、基本変数のいずれかが0であるときに退化が発生します。退化が起こると、単体法の進行が停滞しやすくなります。
巡回:
基底変数と非基底変数の入れ替え操作を行っても、同じ基底解を繰り返し訪れる状態です。これにより、アルゴリズムが無限ループに陥り、最適解にたどり着けないことがあります。
最大係数規則での巡回の確認
実際に最大係数規則で巡回する例を実行してみます。
最初10回を表示すると結果は以下のようになります。
結果から、確かに同じ基底解を選んでしまっているようです。
最大係数規則で巡回する例の実行結果(10ループ回目まで)
m = 3, n = 4
iter: 0
current obj.val. = 0.0
iter = 0
basis = [0, 5, 6]
nonbasis = [4, 1, 2, 3]
rc = [ 1. -4. 0. -8.]
s = 0
d = [1. 0.25 1. ]
x_basis = [0. 0. 1.]
r = 0
iter = 1
basis = [0, 3, 6]
nonbasis = [4, 1, 2, 5]
rc = [ -1. -16. 2. 4.]
s = 3
d = [-12. 1. 4.]
x_basis = [0. 0. 1.]
r = 1
iter = 2
basis = [2, 3, 6]
nonbasis = [4, 1, 0, 5]
rc = [ 0. -8. 1. -4.]
s = 2
d = [1. 0.25 1. ]
x_basis = [0. 0. 1.]
r = 0
iter = 3
basis = [2, 1, 6]
nonbasis = [4, 3, 0, 5]
rc = [ 2. 4. -1. -16.]
s = 1
d = [-12. 1. 4.]
x_basis = [-0. 0. 1.]
r = 1
iter = 4
basis = [4, 1, 6]
nonbasis = [2, 3, 0, 5]
rc = [ 1.00000000e+00 -4.00000000e+00 2.22044605e-16 -8.00000000e+00]
s = 0
d = [1. 0.25 1. ]
x_basis = [ 0. -0. 1.]
r = 0
iter = 5
basis = [4, 5, 6]
nonbasis = [2, 3, 0, 1]
rc = [ -1. -16. 2. 4.]
s = 3
d = [-12. 1. 4.]
x_basis = [0. 0. 1.]
r = 1
iter = 6
basis = [0, 5, 6]
nonbasis = [2, 3, 4, 1]
rc = [ 0. -8. 1. -4.]
s = 2
d = [1. 0.25 1. ]
x_basis = [0. 0. 1.]
r = 0
iter = 7
basis = [0, 3, 6]
nonbasis = [2, 5, 4, 1]
rc = [ 2. 4. -1. -16.]
s = 1
d = [-12. 1. 4.]
x_basis = [0. 0. 1.]
r = 1
iter = 8
basis = [2, 3, 6]
nonbasis = [0, 5, 4, 1]
rc = [ 1. -4. 0. -8.]
s = 0
d = [1. 0.25 1. ]
x_basis = [0. 0. 1.]
r = 0
iter = 9
basis = [2, 1, 6]
nonbasis = [0, 5, 4, 3]
rc = [ -1. -16. 2. 4.]
s = 3
d = [-12. 1. 4.]
x_basis = [-0. 0. 1.]
r = 1
全体のコード
最小添え字規則のコード
simplex_min_index.py
import numpy as np
error = 1.0e-10 # 許容誤差
def print_simplex_detail(
iter=None,
m=None,
n=None,
basis=None,
nonbasis=None,
x_b=None,
y=None,
rc=None,
s=None,
d=None,
r=None,
):
if iter is not None:
print(f"iter = {iter}")
if m is not None and n is not None:
print(f"(m, n) = {(m, n)}")
if basis is not None:
print(f"basis = {basis}")
if nonbasis is not None:
print(f"nonbasis = {nonbasis}")
if y is not None:
print(f"y = {y}")
if rc is not None:
print(f"rc = {rc}")
if s is not None:
print(f"s = {s}")
if d is not None:
print(f"d = {d}")
if x_b is not None:
print(f"x_basis = {x_b}")
if r is not None:
print(f"r = {r}")
def lp_simplex(A, b, c):
(m, n) = A.shape # m:制約条件数, n:変数数
# 初期化
Ai = np.hstack((A, np.identity(m)))
c0 = np.r_[c, np.zeros(m)]
basis = [n + i for i in range(m)]
nonbasis = [j for j in range(n)]
iter = 0
print_simplex_detail(m=m, n=n)
while True:
# 基本解の計算
x = np.zeros(n + m)
x[basis] = np.linalg.solve(Ai[:, basis], b)
# 双対変数の計算
y = np.linalg.solve(Ai[:, basis].T, c0[basis])
# 現在の目的関数の係数を計算
rc = c0[nonbasis] - y @ Ai[:, nonbasis]
# 最適性のチェック（双対可能性）
if np.all(rc <= error):
print("--optimal solution found---------------------------------------------------")
print_simplex_detail(iter=iter, basis=basis, nonbasis=nonbasis, rc=rc, x_b=x[basis])
print("obj.val. = {}".format(c0[basis] @ x[basis]))
print("x = {}".format(x))
print("---------------------------------------------------------------------------")
break
# 入る変数の決定(最小添え字)
s = -1
candidates_nonbasis_index = []
for i in range(len(rc)):
if(rc[i] >= error):
candidates_nonbasis_index.append(i)
s = min(candidates_nonbasis_index)
if s == -1:
print("no entering variable found, optimal solution found")
break
d = np.linalg.solve(Ai[:, basis], Ai[:, nonbasis[s]])
# 非有界性のチェック
if np.all(d <= error):
print("problem is unbounded")
break
# 出る変数の決定(最小添え字)
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
# 最小の比率を持つ変数のインデックスを取得
min_ratio = min(ratio)
candidates_basis_index = [i for i, val in enumerate(ratio) if val == min_ratio]
# 最小添え字を持つインデックスを選択
r = min(candidates_basis_index)
print_simplex_detail(iter=iter, basis=basis, nonbasis=nonbasis, rc=rc, d=d, s=s, r=r, x_b=x[basis])
nonbasis[s], basis[r] = basis[r], nonbasis[s]
if iter % 50 == 0:
print("iter: {}".format(iter))
print("current obj.val. = {}".format(x[basis] @ c0[basis]))
iter += 1
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 巡回する例
A = np.array([[1, 12, -2, -12], [0.25, 1, -0.25, -2], [1, -4, 0, -8]])
b = np.array([0, 0, 1])
c = np.array([1, -4, 0, -8])
lp_simplex(A, b, c)
コードの解説
最大係数規則から最小添え字規則に変更された部分について解説します。
主な変更点はピボット規則くらいです。
ピボット規則
-
candidates_nonbasis_indexリストを初期化します。 -
rc[i]が許容誤差error以上の非基底変数のインデックスを追加します。 -
candidates_nonbasis_indexリストから最も小さい添え字を持つ非基底変数を選択し、そのインデックスをsに設定します。
# 入る変数の決定(最小添え字)
s = -1
candidates_nonbasis_index = []
for i in range(len(rc)):
if(rc[i] >= error):
candidates_nonbasis_index.append(i)
s = min(candidates_nonbasis_index)
- 各基底変数に対して、改善量ベクトル
d[i]が許容誤差errorより大きい場合、比率x[basis][i] / d[i]を計算し、ratioリストに追加します。 -
ratioリストから最小の比率を取得し、その比率を持つ変数のインデックスをcandidates_basis_indexリストに追加します。 -
candidates_basis_indexリストから最も小さい添え字を持つ基底変数を選択し、そのインデックスをrに設定します。
# 出る変数の決定(最小添え字)
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
# 最小の比率を持つ変数のインデックスを取得
min_ratio = min(ratio)
candidates_basis_index = [i for i, val in enumerate(ratio) if val == min_ratio]
# 最小添え字を持つインデックスを選択
r = min(candidates_basis_index)
実行結果
今回は同じ例題と巡回する例を解きます。
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 巡回する例
A = np.array([[1, 12, -2, -12], [0.25, 1, -0.25, -2], [1, -4, 0, -8]])
b = np.array([0, 0, 1])
c = np.array([1, -4, 0, -8])
lp_simplex(A, b, c)
実行結果は以下のようになりました。
目的関数値は16.0、最適化解はx = [2. 1. 3. 0. 0. 0.]です。
目的関数値は1.0、最適化解はx = [1. 0. 1. 0. 1. 0. 0.]です。
最小添え字規則の実行結果
(m, n) = (3, 3)
iter = 0
basis = [3, 4, 5]
nonbasis = [0, 1, 2]
rc = [3. 4. 2.]
s = 0
d = [2. 1. 0.]
x_basis = [4. 8. 6.]
r = 0
iter: 0
current obj.val. = 0.0
iter = 1
basis = [0, 4, 5]
nonbasis = [3, 1, 2]
rc = [-1.5 4. 2. ]
s = 1
d = [0. 0. 3.]
x_basis = [2. 6. 6.]
r = 2
iter = 2
basis = [0, 4, 1]
nonbasis = [3, 5, 2]
rc = [-1.5 -1.33333333 0.66666667]
s = 2
d = [0. 2. 0.33333333]
x_basis = [2. 6. 2.]
r = 1
--optimal solution found---------------------------------------------------
iter = 3
basis = [0, 2, 1]
nonbasis = [3, 5, 4]
rc = [-1.33333333 -1.33333333 -0.33333333]
x_basis = [2. 3. 1.]
obj.val. = 16.0
x = [2. 1. 3. 0. 0. 0.]
---------------------------------------------------------------------------
(m, n) = (3, 4)
iter = 0
basis = [4, 5, 6]
nonbasis = [0, 1, 2, 3]
rc = [ 1. -4. 0. -8.]
s = 0
d = [1. 0.25 1. ]
x_basis = [0. 0. 1.]
r = 0
iter: 0
current obj.val. = 0.0
iter = 1
basis = [0, 5, 6]
nonbasis = [4, 1, 2, 3]
rc = [ -1. -16. 2. 4.]
s = 2
d = [-2. 0.25 2. ]
x_basis = [0. 0. 1.]
r = 1
iter = 2
basis = [0, 2, 6]
nonbasis = [4, 1, 5, 3]
rc = [ 1. 0. -8. -4.]
s = 0
d = [-1. -1. 1.]
x_basis = [0. 0. 1.]
r = 2
--optimal solution found---------------------------------------------------
iter = 3
basis = [0, 2, 4]
nonbasis = [6, 1, 5, 3]
rc = [-1. 0. 0. 0.]
x_basis = [1. 1. 1.]
obj.val. = 1.0
x = [1. 0. 1. 0. 1. 0. 0.]
---------------------------------------------------------------------------
最大改善規則
説明
改善量が最大となるような非基底変数を選びます。
改善量
ここで、
-
rc_i x_i -
d_i x_i -
x_{\text{basis}}[j] x_{\text{basis}}
全体のコード
最大改善規則のコード
simplex_max_improvement.py
import numpy as np
from scipy.linalg import lu_factor, lu_solve
from settings import time_dec
error = 1.0e-10 # 許容誤差
def print_simplex_detail(
iter=None,
m=None,
n=None,
basis=None,
nonbasis=None,
x_b=None,
y=None,
rc=None,
s=None,
d=None,
r=None,
):
if iter is not None:
print(f"iter = {iter}")
if m is not None and n is not None:
print(f"(m, n) = {(m, n)}")
if basis is not None:
print(f"basis = {basis}")
if nonbasis is not None:
print(f"nonbasis = {nonbasis}")
if y is not None:
print(f"y = {y}")
if rc is not None:
print(f"rc = {rc}")
if s is not None:
print(f"s = {s}")
if d is not None:
print(f"d = {d}")
if x_b is not None:
print(f"x_basis = {x_b}")
if r is not None:
print(f"r = {r}")
def check_degeneracy(x_basis):
return np.where(x_basis == 0)[0]
@time_dec
def lp_simplex(A, b, c):
(m, n) = A.shape # m:制約条件数, n:変数数
print("m = {}, n = {}".format(m, n))
# 初期化
Ai = np.hstack((A, np.identity(m)))
c0 = np.r_[c, np.zeros(m)]
basis = [n + i for i in range(m)]
nonbasis = [j for j in range(n)]
iter = 0
while True:
# 基本解の計算
x = np.zeros(n + m)
x[basis] = np.linalg.solve(Ai[:, basis], b)
# 双対変数の計算
y = np.linalg.solve(Ai[:, basis].T, c0[basis])
# 現在の目的関数の係数を計算
rc = c0[nonbasis] - y @ Ai[:, nonbasis]
# LU分解
LU, piv = lu_factor(Ai[:, basis])
# 最適性のチェック（双対可能性）
if np.all(rc <= error):
print("--optimal solution found---------------------------------------------------")
print_simplex_detail(iter=iter)
print("obj.val. = {}".format(c0[basis] @ x[basis]))
print("x = {}".format(x))
print("---------------------------------------------------------------------------")
break
degenerate_indices = check_degeneracy(x[basis])
#退化の確認
if degenerate_indices.size > 0:
# 入る変数の決定(最小添え字)
s = -1
candidates_nonbasis_index = []
for i in range(len(rc)):
if(rc[i] >= error):
candidates_nonbasis_index.append(i)
s = min(candidates_nonbasis_index)
if s == -1:
print("no entering variable found, optimal solution found")
break
d = lu_solve((LU, piv), Ai[:, nonbasis[s]])
# 非有界性のチェック
if np.all(d <= error):
print("problem is unbounded")
break
# 出る変数の決定
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
# 最小の比率を持つ変数のインデックスを取得
min_ratio = min(ratio)
candidates_basis_index = [i for i, val in enumerate(ratio) if val == min_ratio]
# 最小添え字を持つインデックスを選択
r = min(candidates_basis_index)
nonbasis[s], basis[r] = basis[r], nonbasis[s]
else:
# 入る変数の決定(最大改善)
r_max_improvement = -1
s_max_improvement = -1
d_max_improvement = None
max_improvement = -np.inf
for i in range(len(rc)):
if rc[i] > error:
try:
tmp = lu_solve((LU, piv), Ai[:, nonbasis[i]])
except np.linalg.LinAlgError:
print(f"Skipping singular matrix for nonbasis index {i}")
continue
if np.all(tmp <= error):
continue
ratio = []
for j in range(len(tmp)):
if tmp[j] > error:
ratio.append(x[basis][j] / tmp[j])
else:
ratio.append(np.inf)
min_ratio_index = np.argmin(ratio)
improvement = rc[i] * ratio[min_ratio_index]
if improvement > max_improvement:
max_improvement = improvement
s_max_improvement = i
d_max_improvement = tmp
r_max_improvement = min_ratio_index
if s_max_improvement == -1:
print("no entering variable found by maximum improvement, optimal solution found")
break
if r_max_improvement == -1:
print("no entering variable found by maximum improvement, optimal solution found")
break
if np.all(d_max_improvement <= error):
print("problem is unbounded")
break
nonbasis[s_max_improvement], basis[r_max_improvement] = basis[r_max_improvement], nonbasis[s_max_improvement]
if iter % 50 == 0:
print("iter: {}".format(iter))
print("current obj.val. = {}".format(x[basis] @ c0[basis]))
iter += 1
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 計測(karate)
A = np.loadtxt("./test_data/A_karate.txt", delimiter="\t")
m = len(A)
b = np.ones(m)
c = np.loadtxt("./test_data/c_karate.txt", delimiter="\t")
lp_simplex(A, b, c)
# 計測(dolphins)
A = np.loadtxt("./test_data/A_dolphins.txt", delimiter="\t")
m = len(A)
b = np.ones(m)
c = np.loadtxt("./test_data/c_dolphins.txt", delimiter="\t")
lp_simplex(A, b, c)
# 計測(ランダム)
n = 300
m = 500
a = 0.1
b = 0.7
from matrix_generator import generate_random_lp_problem
A, b_vec, c_vec = generate_random_lp_problem(n, m, a, b)
lp_simplex(A, b_vec, c_vec)
コードの解説
最大係数規則から最大改善規則に変更された部分について解説します。
主な変更点以下の通りです。
- LU分解の導入
- 巡回対策
- ピボット規則
LU分解の導入
LU分解をしたいので、scipyを新たにインポートしています。
import numpy as np
from scipy.linalg import lu_factor, lu_solve
from .settings import time_dec
error = 1.0e-10 # 許容誤差
LU分解をしています。dを求めるときは以下のようにして、はやめに計算できます。
# LU分解
LU, piv = lu_factor(Ai[:, basis])
# ... 略
d = lu_solve((LU, piv), Ai[:, nonbasis[s]])
# ... 略
巡回対策
巡回対策は非退化だったら最大改善規則、退化だったら最小添え字規則を行うピボット規則を行います。
退化の判定をします。x_basisに0が含まれていたら退化、いなかったら非退化です。
def check_degeneracy(x_basis):
return np.where(x_basis == 0)[0]
# ... 略
degenerate_indices = check_degeneracy(x[basis])
if degenerate_indices.size > 0:
# 最小添え字規則
else:
# 最大改善規則 + 最急降下(Stepeset Edge)規則
# ... 略
ピボット規則
目的関数値の減少量が最大となるような非基底変数を選びます。各非基底変数に対して、双対解を用いて改善量を計算し、最大の改善量を持つ変数を選びます。
-
sr_max_improvement、s_max_improvement、d_max_improvementを初期化し、max_improvementを負の無限大に設定します。 - 各非基底変数に対して、対応する列ベクトル
tmpを計算します。 - 改善量を計算し、最大の改善量を持つ変数を更新します。
# 入る変数の決定(最大改善)
r_max_improvement = -1
s_max_improvement = -1
d_max_improvement = None
max_improvement = -np.inf
for i in range(len(rc)):
if rc[i] > error:
try:
tmp = lu_solve((LU, piv), Ai[:, nonbasis[i]])
except np.linalg.LinAlgError:
print(f"Skipping singular matrix for nonbasis index {i}")
continue
if np.all(tmp <= error):
continue
ratio = []
for j in range(len(tmp)):
if tmp[j] > error:
ratio.append(x[basis][j] / tmp[j])
else:
ratio.append(np.inf)
min_ratio_index = np.argmin(ratio)
improvement = rc[i] * ratio[min_ratio_index]
if improvement > max_improvement:
max_improvement = improvement
s_max_improvement = i
d_max_improvement = tmp
r_max_improvement = min_ratio_index
# ... 略
nonbasis[s_max_improvement], basis[r_max_improvement] = basis[r_max_improvement], nonbasis[s_max_improvement]
実行結果
最大係数規則と同じ例題を解きます。
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
実行結果は以下のようになりました。目的関数値は16.0、最適化解はx = [2. 1. 3. 0. 0. 0.]です。
m = 3, n = 3
iter: 0
current obj.val. = 0.0
--optimal solution found---------------------------------------------------
iter = 3
basis = [0, 2, 1]
nonbasis = [3, 5, 4]
rc = [-1.33333333 -1.33333333 -0.33333333]
x_basis = [2. 3. 1.]
obj.val. = 16.0
x = [2. 1. 3. 0. 0. 0.]
---------------------------------------------------------------------------
最急降下(Steeppest Edge)規則
説明
最急降下規則では、目的関数の減少量が最大となるような非基底変数を選びます。この規則では、各非基底変数に対して、改善量をその変数に対応する列ベクトルの調整ノルムで割った値として計算し、その中で最大のものを選びます。これにより、改善が最も急な方向を選択することができます。(あまり自信ない;;)
改善量
ここで、
-
rc_i x_i -
d_{ij} x_i
調整ノルムは次のように定義されます：
全体のコード
最急降下規則のコード
simplex_steepest_edge.py
import numpy as np
from scipy.linalg import lu_factor, lu_solve
from settings import time_dec
error = 1.0e-10 # 許容誤差
def print_simplex_detail(
iter=None,
m=None,
n=None,
basis=None,
nonbasis=None,
x_b=None,
y=None,
rc=None,
s=None,
d=None,
r=None,
):
if iter is not None:
print(f"iter = {iter}")
if m is not None and n is not None:
print(f"(m, n) = {(m, n)}")
if basis is not None:
print(f"basis = {basis}")
if nonbasis is not None:
print(f"nonbasis = {nonbasis}")
if y is not None:
print(f"y = {y}")
if rc is not None:
print(f"rc = {rc}")
if s is not None:
print(f"s = {s}")
if d is not None:
print(f"d = {d}")
if x_b is not None:
print(f"x_basis = {x_b}")
if r is not None:
print(f"r = {r}")
def check_degeneracy(x_basis):
return np.where(x_basis == 0)[0]
@time_dec
def lp_simplex(A, b, c):
(m, n) = A.shape # m:制約条件数, n:変数数
print("m = {}, n = {}".format(m, n))
# 初期化
Ai = np.hstack((A, np.identity(m)))
c0 = np.r_[c, np.zeros(m)]
basis = [n + i for i in range(m)]
nonbasis = [j for j in range(n)]
iter = 0
while True:
# 基本解の計算
x = np.zeros(n + m)
x[basis] = np.linalg.solve(Ai[:, basis], b)
# 双対変数の計算
y = np.linalg.solve(Ai[:, basis].T, c0[basis])
# 現在の目的関数の係数を計算
rc = c0[nonbasis] - y @ Ai[:, nonbasis]
# LU分解
LU, piv = lu_factor(Ai[:, basis])
# 最適性のチェック（双対可能性）
if np.all(rc <= error):
print("--optimal solution found---------------------------------------------------")
print_simplex_detail(iter=iter)
print("obj.val. = {}".format(c0[basis] @ x[basis]))
print("x = {}".format(x))
print("---------------------------------------------------------------------------")
break
degenerate_indices = check_degeneracy(x[basis])
if degenerate_indices.size > 0:
# 入る変数の決定(最小添え字)
s = -1
candidates_nonbasis_index = []
for i in range(len(rc)):
if(rc[i] >= error):
candidates_nonbasis_index.append(i)
s = min(candidates_nonbasis_index)
if s == -1:
print("no entering variable found, optimal solution found")
break
#d = np.linalg.solve(Ai[:, basis], Ai[:, nonbasis[s]])
d = lu_solve((LU, piv), Ai[:, nonbasis[s]])
# 非有界性のチェック
if np.all(d <= error):
print("problem is unbounded")
break
# 出る変数の決定
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
# 最小の比率を持つ変数のインデックスを取得
min_ratio = min(ratio)
candidates_basis_index = [i for i, val in enumerate(ratio) if val == min_ratio]
# 最小添え字を持つインデックスを選択
r = min(candidates_basis_index)
nonbasis[s], basis[r] = basis[r], nonbasis[s]
else:
# 入る変数の決定(最急降下)
s_steepest_edge = -1
d_steepest_edge = None
max_ratio = -np.inf
for i in range(len(rc)):
if rc[i] > error:
tmp = lu_solve((LU, piv), Ai[:, nonbasis[i]])
norm = np.sqrt(1 + np.sum(tmp**2))
ratio = rc[i] / norm
if ratio > max_ratio:
max_ratio = ratio
s_steepest_edge = i
d_steepest_edge = tmp
if s_steepest_edge == -1:
print("no entering variable found by steepest edge, optimal solution found")
break
if np.all(d_steepest_edge <= error):
print("problem is unbounded")
break
# 最急辺規則で選んだ変数の比率を計算
ratio = []
for i in range(len(d_steepest_edge)):
if d_steepest_edge[i] > error:
ratio.append(x[basis][i] / d_steepest_edge[i])
else:
ratio.append(np.inf)
# 出る変数の決定
r_steepest_edge = np.argmin(ratio)
nonbasis[s_steepest_edge], basis[r_steepest_edge] = basis[r_steepest_edge], nonbasis[s_steepest_edge]
if iter % 50 == 0:
print("iter: {}".format(iter))
print("current obj.val. = {}".format(x[basis] @ c0[basis]))
iter += 1
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 計測(karate)
A = np.loadtxt("./test_data/A_karate.txt", delimiter="\t")
m = len(A)
b = np.ones(m)
c = np.loadtxt("./test_data/c_karate.txt", delimiter="\t")
lp_simplex(A, b, c)
# 計測(dolphins)
A = np.loadtxt("./test_data/A_dolphins.txt", delimiter="\t")
m = len(A)
b = np.ones(m)
c = np.loadtxt("./test_data/c_dolphins.txt", delimiter="\t")
lp_simplex(A, b, c)
# 計測(ランダム)
n = 300
m = 500
a = 0.1
b = 0.7
from matrix_generator import generate_random_lp_problem
A, b_vec, c_vec = generate_random_lp_problem(n, m, a, b)
lp_simplex(A, b_vec, c_vec)
コードの解説
最大改善規則から最急降下規則に変更された部分について解説します。
主な変更点はピボット規則くらいです。
ピボット規則
最急降下規則では、目的関数の減少量が最大となるような非基底変数を選びます。この規則では、各非基底変数に対して、改善量
-
sとdを初期化し、max_ratioを負の無限大に設定します。 - 各非基底変数に対して、対応する列ベクトル
tmpを計算します。 - 列ベクトル
tmpの調整ノルムを計算します。 - 相対コストを調整ノルムで割った値として改善量を計算し、最大の改善量を持つ変数を更新します。
# 入る変数の決定(最急降下)
s_steepest_edge = -1
d_steepest_edge = None
max_ratio = -np.inf
for i in range(len(rc)):
if rc[i] > error:
tmp = lu_solve((LU, piv), Ai[:, nonbasis[i]])
norm = np.sqrt(1 + np.sum(tmp**2))
ratio = rc[i] / norm
if ratio > max_ratio:
max_ratio = ratio
s_steepest_edge = i
d_steepest_edge = tmp
# ... 略
nonbasis[s_steepest_edge], basis[r_steepest_edge] = basis[r_steepest_edge], nonbasis[s_steepest_edge]
実行結果
今回も同じ例題を解きます。
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 巡回する例
#A = np.array([[1, 12, -2, -12], [0.25, 1, -0.25, -2], [1, -4, 0, -8]])
#b = np.array([0, 0, 1])
#c = np.array([1, -4, 0, -8])
#lp_simplex(A, b, c)
実行結果は以下のようになりました。目的関数値は16.0、最適化解はx = [2. 1. 3. 0. 0. 0.]です。
m = 3, n = 3
iter: 0
current obj.val. = 0.0
--optimal solution found---------------------------------------------------
iter = 3
basis = [0, 2, 1]
nonbasis = [3, 5, 4]
rc = [-1.33333333 -1.33333333 -0.33333333]
x_basis = [2. 3. 1.]
obj.val. = 16.0
x = [2. 1. 3. 0. 0. 0.]
---------------------------------------------------------------------------
計測
いくつかのデータで今回実装した最大係数規則、最大改善規則、最急降下規則(Steeppest Edge)、そしてPython-MIPを比較してみようと思います。
最大係数規則は巡回対策を行っていないので巡回対策した以下のコードを使用します。
巡回対策した最大係数規則のコード
simplex_revised_max_rc.py
import numpy as np
from scipy.linalg import lu_factor, lu_solve
from settings import time_dec
error = 1.0e-10 # 許容誤差
def print_simplex_detail(
iter=None,
m=None,
n=None,
basis=None,
nonbasis=None,
x_b=None,
y=None,
rc=None,
s=None,
d=None,
r=None,
):
if iter is not None:
print(f"iter = {iter}")
if m is not None and n is not None:
print(f"(m, n) = {(m, n)}")
if basis is not None:
print(f"basis = {basis}")
if nonbasis is not None:
print(f"nonbasis = {nonbasis}")
if y is not None:
print(f"y = {y}")
if rc is not None:
print(f"rc = {rc}")
if s is not None:
print(f"s = {s}")
if d is not None:
print(f"d = {d}")
if x_b is not None:
print(f"x_basis = {x_b}")
if r is not None:
print(f"r = {r}")
def check_degeneracy(x_basis):
return np.where(x_basis == 0)[0]
@time_dec
def lp_simplex(A, b, c):
(m, n) = A.shape # m:制約条件数, n:変数数
print("m = {}, n = {}".format(m, n))
# 初期化
Ai = np.hstack((A, np.identity(m)))
c0 = np.r_[c, np.zeros(m)]
basis = [n + i for i in range(m)]
nonbasis = [j for j in range(n)]
iter = 0
while True:
# 基本解の計算
x = np.zeros(n + m)
x[basis] = np.linalg.solve(Ai[:, basis], b)
# 双対変数の計算
y = np.linalg.solve(Ai[:, basis].T, c0[basis])
# 現在の目的関数の係数を計算
rc = c0[nonbasis] - y @ Ai[:, nonbasis]
# LU分解
LU, piv = lu_factor(Ai[:, basis])
# 最適性のチェック（双対可能性）
if np.all(rc <= error):
print("--optimal solution found---------------------------------------------------")
print_simplex_detail(iter=iter)
print("obj.val. = {}".format(c0[basis] @ x[basis]))
print("x = {}".format(x))
print("---------------------------------------------------------------------------")
break
degenerate_indices = check_degeneracy(x[basis])
if degenerate_indices.size > 0:
# 入る変数の決定(最小添え字)
s = -1
candidates_nonbasis_index = []
for i in range(len(rc)):
if(rc[i] >= error):
candidates_nonbasis_index.append(i)
s = min(candidates_nonbasis_index)
if s == -1:
print("no entering variable found, optimal solution found")
break
#d = np.linalg.solve(Ai[:, basis], Ai[:, nonbasis[s]])
d = lu_solve((LU, piv), Ai[:, nonbasis[s]])
# 非有界性のチェック
if np.all(d <= error):
print("problem is unbounded")
break
# 出る変数の決定
ratio = []
for i in range(len(d)):
if d[i] > error:
ratio.append(x[basis][i] / d[i])
else:
ratio.append(np.inf)
# 最小の比率を持つ変数のインデックスを取得
min_ratio = min(ratio)
candidates_basis_index = [i for i, val in enumerate(ratio) if val == min_ratio]
# 最小添え字を持つインデックスを選択
r = min(candidates_basis_index)
nonbasis[s], basis[r] = basis[r], nonbasis[s]
else:
# 入る変数の決定(最大係数)
s_max_rc = np.argmax(rc)
d_max_rc = np.linalg.solve(
Ai[:, basis], Ai[:, nonbasis[s_max_rc]]
)
if s_max_rc == -1:
print("no entering variable found by steepest edge, optimal solution found")
break
if np.all(d_max_rc <= error):
print("problem is unbounded")
break
# 最急辺規則で選んだ変数の比率を計算
ratio = []
for i in range(len(d_max_rc)):
if d_max_rc[i] > error:
ratio.append(x[basis][i] / d_max_rc[i])
else:
ratio.append(np.inf)
# 出る変数の決定
r_max_rc = np.argmin(ratio)
nonbasis[s_max_rc], basis[r_max_rc] = basis[r_max_rc], nonbasis[s_max_rc]
if iter % 50 == 0:
print("iter: {}".format(iter))
print("current obj.val. = {}".format(x[basis] @ c0[basis]))
iter += 1
if __name__ == "__main__":
# 例題
A = np.array([[2, 0, 0], [1, 0, 2], [0, 3, 1]])
b = np.array([4, 8, 6])
c = np.array([3, 4, 2])
lp_simplex(A, b, c)
# 計測(karate)
A = np.loadtxt("./test_data/A_karate.txt", delimiter="\t")
m = len(A)
b = np.ones(m)
c = np.loadtxt("./test_data/c_karate.txt", delimiter="\t")
lp_simplex(A, b, c)
# 計測(dolphins)
A = np.loadtxt("./test_data/A_dolphins.txt", delimiter="\t")
m = len(A)
b = np.ones(m)
c = np.loadtxt("./test_data/c_dolphins.txt", delimiter="\t")
lp_simplex(A, b, c)
# 計測(ランダム)
n = 300
m = 500
a = 0.1
b = 0.7
from matrix_generator import generate_random_lp_problem
A, b_vec, c_vec = generate_random_lp_problem(n, m, a, b)
lp_simplex(A, b_vec, c_vec)
環境
今回の計測には"AMD Ryzen 5 5625U with Radeon Graphics 2.30 GHz"を使用しています。
計測結果
空手データ
| obj.val. | 反復回数(回) | 実行時間(秒) | 実行時間/反復回数 | |
|---|---|---|---|---|
| Python-MIP | 0.4197896120986266 | 29 | 0.02504110336303711 | 0.00086348632 |
| 最大係数 | 0.4197896120973046 | 57 | 0.04264235496520996 | 0.00074811149 |
| 最大改善 | 0.4197896120973046 | 57 | 0.0494387149810791 | 0.00086734587 |
| 最急降下 | 0.4197896120973046 | 58 | 0.047457218170166016 | 0.00081822789 |
ドルフィンデータ
| obj.val. | 反復回数(回) | 実行時間(秒) | 実行時間/反復回数 | |
|---|---|---|---|---|
| Python-MIP | 0.5285194414840361 | 63 | 0.057343482971191406 | 0.00091021401 |
| 最大係数 | 0.52851944147779 | 4905 | 8.825235366821289 | 0.00179923249 |
| 最大改善 | 0.5285194414777898 | 2822 | 5.089014053344727 | 0.00180333595 |
| 最急降下 | 0.5285194414777898 | 2311 | 4.076169729232788 | 0.00176381208 |
ランダムデータ(500 \times 300
| obj.val. | 反復回数(回) | 実行時間(秒) | 実行時間/反復回数 | |
|---|---|---|---|---|
| Python-MIP | 2.199968233370979 | 246 | 1.1082000732421875 | 0.00450487834 |
| 最大係数 | 2.1999467451415566 | 979 | 34.45111870765686 | 0.03519011103 |
| 最大改善 | 2.1999467451416312 | 404 | 217.79493355751038 | 0.53909637019 |
| 最急降下 | 2.1999467451414434 | 203 | 6.295671224594116 | 0.03101315874 |
Discussion