1. 目的
業務で使用することが多い線形計画法について理解した内容をゆるふわにまとめつつ、numpyとモデリングツールで実装し理解を深めます。間違いがあればご指摘頂ければ幸いです。
2. 線形計画法について
線形計画法の標準形は下式で定義される。
\begin{eqnarray}
\min & &c^Tx \\
s.t. & &Ax=b\\
& &x\geq 0
\end{eqnarray}
\tag{1}
双対問題は下式で与えられる。
\begin{eqnarray}
\max & &b^Ty \\
s.t. & &A^Ty-c+z=0\\
& &z\geq 0
\end{eqnarray}
\tag{2}
双対問題の導出は下記を参照ください。
主問題と双対問題の間には、弱双対定理、双対定理、相補性定理などの重要な関係がある。
2.1 弱双対定理
弱双対定理は、主問題である最小化問題と双対問題である最大化問題の大小関係を示したものである。$x$と$y$が主問題と双対問題の実行可能解であるならば、
\begin{eqnarray}
f_p = c^Tx \geq b^Ty = f_d
\end{eqnarray}
\tag{3}
が成り立つ。
- 証明
\begin{aligned}
c^Tx - b^Ty &=c^Tx - (Ax)^Ty \\
&=c^Tx - x^T A^Ty \\
&=x^Tc - x^T A^Ty \\
&=x^T(c - A^Ty) \\
&=x^Tz \\
&\geq 0
\end{aligned}
- 弱双対定理からわかること
大小関係から、主問題の目的関数の最小値は双対問題で抑えられることがわかる。さらに、双対問題が上に有界でないとき、つまり$b^Ty\rightarrow\infty$のとき、主問題に実行可能解が存在しないことがわかる。
大小関係から、双対問題の目的関数の最大値は主問題で抑えられることがわかる。さらに、主問題が下に有界でないとき、つまり$c^Tx\rightarrow-\infty$のとき、双対問題に実行可能解が存在しないことがわかる。
大小関係から、主問題の目的関数の最小値と双対問題の目的関数の最大値が存在するとき目的関数値は一致することがわかる。
2.2 双対定理
主問題が最適解をもつならば、双対問題も解をもち、$b^Ty = c^Tx$が成り立つ。
-
証明
Farkasの定理を用いて$w=b^Ty - c^Tx \geq 0$となる$w$が存在することを証明する。 -
双対定理からわかること
双対定理の結果から、双対ギャップ$w=b^Ty - c^Tx=x^Tz$に着目することで、数値解析結果が最適解に近いか判断できる。
2.3 相補性定理
$x$と$y$がそれぞれ主問題と双対問題の最適解になるための必要十分条件は、次の三つの条件が成り立つことである。この三つの条件はKKT条件に相当する。1つ目の条件式は双対問題の実行可能性、２つ目の条件式は主問題の実行可能性、３つ目の条件式は相補性条件である。
\begin{eqnarray}
&& A^Ty-c+z=0&, z \geq 0 \\
&& Ax = b &, x \geq 0 \\
&& x_i z_i=0
\end{eqnarray}
- KKT条件からの導出
一般の制約付き最小化問題に関する最適解の条件を考える。
\begin{eqnarray}
\min & &f(x) \\
s.t. & &h_i(x)\leq 0\\
& &g_j(x) = 0
\end{eqnarray}
このときのKKT条件は5つの条件で与えられる。
\begin{aligned}
\nabla f(x) + \sum_{i=1}^m y_i &\nabla g_i(x) + \sum_{i=1}^l z_i \nabla h_i(x)=0,\\
g_i(x) = 0,\\
h_i(x) \leq 0,\\
z_i \geq 0, \\
z_i h_i (x) = 0
\end{aligned}
線形計画のKKT条件は次式で与えられる。線形計画のKKT条件と相補性定理と一致することがわかる。
\begin{aligned}
\nabla f(x) + \sum_{i=1}^m y_i &\nabla g_i(x) + \sum_{i=1}^l z_i \nabla h_i(x) \\
= \nabla(c^Tx) + \nabla &y^T(b - Ax) + \nabla z^T(-x) \\
= c - A^Ty -z &= 0,\\[7pt]
Ax&=b,\\
-x &\leq 0, \\
z &\geq 0,\\
z_i x_i&=0
\end{aligned}
KKT条件の導出については下記の記事を参照ください。
3. 数値解析について
数値解析は主に単体法と内点法がある。
単体法は、代数的な構造に着目した数値解析である。単体法は、実行可能領域の境界を移動することで最適解を求める。内点法は、非線形計画法の考え方を援用し、実行可能解領域に内点を生成し最適解を求める数値解析である。
3.1. 単体法について
単体法は、代数的な構造に着目した数値解析である。
単体法の派生的な数値解析として二段階法や罰金法などがある。
単体法は、実行可能領域の境界を移動することで最適解を求める。初期段階で実行可能基底解をまず定める必要があるのだが、必ずしも見つかるわけではないという課題があった。
この課題を解決するため、2段階法と罰金法が提案された。2段階法は名前からもわかるように、人工変数を導入することで1段階目で初期実行可能基底解を算出し、2段階目で単体法を適用することで、単体法の抱える課題を解決した。
2段階法では2段階を踏む数値解析であったが、罰金法は目的関数に人口変数をペナルティ項として導入することで、2段階を踏まずに連続的に最適解を算出できるようにした数値解析である。
詳細は下記の記事を参照ください。
3.2. 内点法について
内点法は、非線形計画法の考え方を援用し、実行可能解領域に内点を生成し最適解を求める数値解析である。内点法には主変数と双対変数を対等に扱う主双対内点法、主問題に適用される主内点法、双対問題に適用される双対内点法がある。また、主双対内点法には探索ベクトル$[\Delta x,\Delta y,\Delta z]$の算出方法やステップ幅$\alpha_k \in (0,1]$の算出方法によって種類がある。主双対アフィンスケーリング法、主双対パス追跡法、主双対ポテンシャル減少法などが挙げられる。各主双対内点法をできるだけnumpyで実装した。GitHubにアップロードしたのでよかったらご参考ください。
下記の主問題を代表として、実行結果を考察する。
\begin{eqnarray}
\min & & -5x_1 -4x_2 \\
s.t. & &5x_1 + 2x_2 \leq 30\\
& &x_1 + 2x_2 \leq 14\\
& &x_1\geq 0, x_2\geq 0
\end{eqnarray}
不等式標準形から等式標準形に変換する。変換するためにスラック変数$x_3\geq0$と$x_4\geq0$を導入する。
\begin{eqnarray}
\min & & -5x_1 - 4x_2 + 0 x_3 + 0 x_4\\
s.t. & & 5x_1 + 2x_2 +1x_3+0x_4= 30\\
& & \ \ x_1 + 2x_2 +0x_3+1x_4= 14\\
& & \ \ x_1\geq 0, x_2\geq 0, x_3\geq 0, x_4\geq 0
\end{eqnarray}
主双対アフィンスケーリング法、主双対パス追跡法、主双対ポテンシャル減少法における双対ギャップの収束結果と１イテレーション毎の実行時間を示す。10回程度の反復で双対ギャップが1e-7まで収束している。主双対内点法のなかで主双対ポテンシャル減少法の実行時間が長い。ポテンシャル関数に基づくステップ幅の最適化に時間を要したためと考えられる。
各節で、各主双対内点法の実行結果を示す。
@taka_horibe様のデザインを踏襲しました。
3.2.1 主双対パス追跡法
主双対パス追跡法の実行結果を示す。左図は実行可能領域における各イテレーションにおける内点の位置を表している。右上図は各イテレーションにおける主問題と双対問題の目的関数値を示す。右下図は各イテレーションにおける双対ギャップを示す。
3.2.2 主双対アフィンスケーリング法
主双対アフィンスケーリング法の実行結果を示す。左図は実行可能領域における各イテレーションにおける内点の位置を表している。右上図は各イテレーションにおける主問題と双対問題の目的関数値を示す。右下図は各イテレーションにおける双対ギャップを示す。
3.2.3 主双対ポテンシャル減少法
主双対ポテンシャル減少法の実行結果を示す。左図は実行可能領域における各イテレーションにおける内点の位置を表している。右上図は各イテレーションにおける主問題と双対問題の目的関数値を示す。右下図は各イテレーションにおける双対ギャップを示す。
3.2.1 主双対パス追跡法、再び
本記事では、主双対パス追跡法を例に主双対内点法の理解を深める。さて、線形計画の最適性条件はKKT条件から与えられるのであった。
\begin{aligned}
&A^Ty-c+z=0 \\
&Ax = b \\
&x_i z_i=0 \\
& x \geq 0, z \geq 0
\end{aligned}
\tag{4}
内点法では最適解に近づく前に非負制約の境界に近づくと、そこで収束が停滞してしまう可能性がある。できるだけ非負制約の境界境界から遠ざけるため、解析的中心という概念を導入する。下式を中心化KKT条件という。
\begin{aligned}
&A^Ty-c+z=0 \\
&Ax = b \\
&x_i z_i=\mu \\
& x > 0, z > 0
\end{aligned}
\tag{5}
$\mu$の定め方として、$||XZe - \mu e||$が最小となるようにとることで、現在の反復に近い解析点を求める。$||\cdot||$をユークリッドノルムとすると、$\mu$は下式で与えられる。
\begin{aligned}
&\dfrac{\partial}{\partial \mu}\sum_{i=1}^{n} (x_i z_i - \mu)^2 = 0 \\
& \sum_{i=1}^{n} 2 (x_i z_i - \mu)(-1) = 0 \\
& \sum_{i=1}^{n} x_i z_i = \sum_{i=1}^{n} \mu = n \mu \\
& \mu = \dfrac{1}{n}\sum_{i=1}^{n} x_i z_i = \dfrac{1}{n}x^Tz
\end{aligned}
\tag{6}
ここで、注目すべきは最適解に近づけば、双対ギャップ$x^Tz \rightarrow 0$となるので、$\mu \rightarrow 0$に近づく。そのため、$x_i z_i = \mu$は相補性条件$x_i z_i = 0$に近づくことがわかる。
ニュートン法を用いて、式5の非線形方程式を解き最適解を求める。
初期値として$(x_k,y_k,z_k)$を与える。初期値から$(\Delta x,\Delta y,\Delta z)$方向に微小に動かすことで、非線形方程式を満足させる。ここで、問題になるのがどれくらい微小に動かせば良いかである。
\begin{aligned}
&A^T(y_k +\Delta y)-c+(z_k +\Delta z)=0 \\
&A(x_k +\Delta x) = b \\
&(x_i+ \Delta x_i )(z_i+ \Delta z_i )=\sigma\mu_k \\
& x \geq 0, z \geq 0
\end{aligned}
\tag{7}
微小に動かしているので$\Delta x_i \Delta z_i \simeq0$とすると下式で表せられる。
\begin{aligned}
&A^T(y_k +\Delta y)-c+(z_k +\Delta z)=0 \\
&A(x_k +\Delta x) = b \\
&x_i z_i + \Delta x_i z_i + \Delta z_i x_i=\sigma\mu_k \\
& x \geq 0, z \geq 0
\end{aligned}
\tag{8}
探索方向$(\Delta x,\Delta y,\Delta z)$は下式の連立1次方程式を解くことで求められる。
\begin{aligned}
\left[
\begin{array}{lll}
O & A^T & I \\
A & O & O \\
Z_k & O & X_k
\end{array}
\right]
\left[
\begin{array}{l}
\Delta x \\
\Delta y \\
\Delta z
\end{array}
\right]
=
-
\left[
\begin{array}{r}
A^Ty_k-c+z_k &\\
Ax_k - b & \\
X_k Z_k e -\sigma\mu_k e&
\end{array}
\right]
\end{aligned}
\tag{9}
\\
X_k = diag(x_{1_k},\cdots,x_{n_k}) \\
Z_k = diag(z_{1_k},\cdots,z_{n_k}) \\
e = [1,\cdots,1]^T
より演算時間を高速化する手法もあるが本記事では素朴に解くことにする。
微小な方向に初期点を更新して、非線形方程式を満足するようにする。ただし、更新する際に非負制約$x\geq 0$、$z\geq 0$を満足するようにステップ幅$\alpha$を導入する。
\begin{aligned}
x_{k+1} \leftarrow x_k + \alpha \Delta x \\
y_{k+1} \leftarrow y_k + \alpha \Delta y \\
z_{k+1} \leftarrow z_k + \alpha \Delta z
\end{aligned}
\tag{10}
微小な方向に更新しても非線形方程式を満足しない可能性があるため、前の工程を踏襲して反復して更新し、最適解を求める。
k \leftarrow k+1
\tag{11}
しれっと式7に$\sigma$を導入しているが、$\sigma\in[0,1)$を導入することで、式12のように$\mu_k$を反復する度に指数関数的に減衰させる。つまり、$k\rightarrow \infty$としたとき、$\mu_k\rightarrow 0$となり、双対ギャップ$c^Tx-b^Ty=x^Tz=n\mu$が指数関数的に減少していく。
\begin{aligned}
\mu_{k+1} = c\mu_{k} =c^k\mu_{1}\\
c := 1 - \alpha\left(1-\sigma\right)\\
0<c<1
\end{aligned}
\tag{12}
4. モデリングツール
3章ではnumpyで実装していったが、実務ではモデリングツールを用いる。モデリングツールを用いると簡潔に記述ができる。
| 言語 | ツール | 線形計画 | 整数線形計画 | 混合整数線形計画 |
|---|---|---|---|---|
| Python | PulP | 〇 | ||
| Python | cvxpy | 〇 | 〇 | 〇 |
| Python | scipy | 〇 | ||
| Numerical Optimizer | 〇 | |||
| Matlab | Optimization toolbox | 〇 |
4.1 PulPの実装例
import pulp
# 線形計画法の定義
problem = pulp.LpProblem("test",pulp.LpMinimize)
x1 = pulp.LpVariable("x1",0,10,"Continuous")
x2 = pulp.LpVariable("x2",0,10,"Continuous")
# 目的関数をセット
problem += -5*x1 - 4*x2
# 制約条件をセット
problem += 5*x1 + 2*x2 <=30
problem += 1*x1 + 2*x2 <=14
# 解く
result = problem.solve()
# 結果を出力
print(problem)
print(x1.value())
print(x2.value())
下記に実行結果を示す。
test:
MINIMIZE
-5*x1 + -4*x2 + 0
SUBJECT TO
_C1: 5 x1 + 2 x2 <= 30
_C2: x1 + 2 x2 <= 14
VARIABLES
x1 <= 10 Continuous
x2 <= 10 Continuous
4.0
5.0
4.2 cvxpyの実装例
import cvxpy as cp
import numpy as np
A = np.array([[5,2],[1,2]])
b = np.array([30,14])
c = np.array([-5,-4])
# Define and solve the CVXPY problem.
n = len(c)
x = cp.Variable(n)
prob = cp.Problem(cp.Minimize(c.T@x),[A @ x <= b])
prob.solve()
# Print result.
print("\nThe optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("A dual solution is")
print(prob.constraints[0].dual_value)
下記に実行結果を示す。
$ The optimal value is -39.999999976234555
$ A solution x is
$ [4. 5.]
$ A dual solution is
$ [0.75000001 1.25000001]
4.3 scipyの実装例
from scipy.optimize import linprog
c = [-1, 4]
A = [[-3, 1], [1, 2]]
b = [6, 4]
x0_bounds = (None, None)
x1_bounds = (-3, None)
res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])
下記に実行結果を示す。
$ print(res)
con: array([], dtype=float64)
fun: -21.99999984082494 # may vary
message: 'Optimization terminated successfully.'
nit: 6 # may vary
slack: array([3.89999997e+01, 8.46872439e-08] # may vary
status: 0
success: True
x: array([ 9.99999989, -2.99999999]) # may vary
5. 参照文献
[1] 工学基礎 最適化とその応用、数理工学社、矢部博
[2] システム制御のための最適化理論、コロナ社、延山、瀬部
[3] 東工大 水野研究室資料
6 おわりに
実際に実装したことで理解が深まった。課題として単体法の実装方法がわからなかった。