#freeze
#contents
2010/03/19からのアクセス回数 &counter;
ここでは、sage上でCVXOPTを使って最適解問題を解く方法に説明します。
このページのsageノートブックは、以下のURLで見ることができます。
[[http://www.sagenb.org/home/pub/1813]]
** 線形最適解問題 [#y285d0e5]
最初にCVXOPTのTUTORIALからSolving a linear program
をsageで解いてみます。
以下の問題がこれから解く、線形最適解問題です。
$$
\begin{array}{ll}
\mbox{minimize} & 2x_1 + x_2 \\
\mbox{subject to} & -x_1 + x_2 \leq 1 \\
& x_1 + x_2 \geq 2 \\
& x_2 \geq 0 \\
& x_1 -2x_2 \leq 4
\end{array}
$$
ここで、minimize f(w)というのは、目的関数\(f(w)\)を示し、subject to \(g(w) \le 0\)は、不等式制約\(g(w)\)を示します。
従ってこの最適化問題は、不等式制約のもとで目的関数が最も小さくなる値を求めることです。
*** 目的関数と制約条件の可視化 [#g7723967]
CVXOPTを使って解く前に、sageの可視化ツールを使って目的関数と制約条件を表示してみましょう。
黄色の領域が制約条件で囲まれた部分です。\( y = 2 x + c \) の直線で、cの値の最も小さいものが、
赤の直線です。
sageへの入力:
#pre{{
# 制約条件のプロット
rgn = [(5,6), (0.5, 1.5), (2,0), (4, 0), (5,0.5), (5,6)]
poly = polygon(rgn, rgbcolor='yellow')
plt1 = plot(lambda x :(x+1), (x, 0, 5))
plt2 = plot(lambda x :(-x+2), (x, 0, 5))
plt3 = plot(lambda x :(0.5*x-2), (x, 0, 5))
# 解の目的関数のプロット
solv = plot(lambda x :(-2*x+2.5), (x, 0, 3), rgbcolor='red')
(plt1+plt2+plt3+poly+solv).show(xmin=0, xmax=5)
}}
&ref(1.png);
*** CVXOPTに必要なインポート [#dcc4eb4c]
まずは、CVXOPTに必要なパッケージをインポートします。
sageへの入力:
#pre{{
# 線形最適解問題
# CVXOPTを使用するのに必要なインポート
# sageではcvxopt.baseからmatrixをインポートする
#
from cvxopt import solvers
from cvxopt.base import matrix as _matrix
RealNumber = float
Integer=int
}}
*** solvers.lpを使って解く [#ccc62264]
SVXOPTの本では線形最適解問題の一般問題は、次のような形式で与えられます。
不等式制約\( G x \leq h \)$ と 等式制約 \( A x = b \)の混在形式。
$$
\begin{array}{ll}
\mbox{minimize} & c^T x + d \\
\mbox{subject to} & G x \leq h \\
& A x = b
\end{array}
$$
これをSVXOPTでは、\( s \geq 0\) を導入して、 \( G x + s = h, s \geq 0 \)に分割し、
目的関数のd（定数）を除いた以下の形式で解いています。
$$
\begin{array}{ll}
\mbox{minimize} & c^T x \\
\mbox{subject to} & G x + s = h \\
& A x = b \\
& s \geq 0
\end{array}
$$
solvers.lp()を使って線形最適解問題を解いてみます。lpの使い方をみると、
以下のような引数をとり、不等号制約条件問題では、c, G, hを与えればよいこと
が分かります。
#pre{{
solvers.lp(c, G, h, A, b, solver, primalstart, dualstart)
minimize c'*x
subject to G*x + s = h
A*x = b
s >= 0
}}
従って、例題をsolvers.lpに与える引数の形式に整理すると、
$$
\begin{array}{lll}
c^T & = & \left( \begin{array}{cl} 2 & 1 \end{array} \right) \\
G & = & \left( \begin{array}{rrrr} -1 & 1 \\ -1 & -1 \\ 0 & -1 \\ 1 & -2 \end{array} \right) \\
h & = & \left( \begin{array}{r} 1 \\ -2 \\ 0 \\ 4 \end{array} \right) \end{array}
$$
sageへの入力:
#pre{{
# solvers.lpの使い方を調べます（シフト・リターンでヘルプが表示されます）
solvers.lp?
}}
sageへの入力:
#pre{{
# CVXOPTのmaxtrixを作成
# _matrixの要素の与え方がsageと異なることに注意（縦のベクトル要素を並べた形）
c = _matrix([2., 1.])
G = _matrix([[-1., -1., 0., 1.], [1., -1., -1., -2.]])
h = _matrix([1., -2., 0., 4.])
}}
*** 線形最適化例題を解く [#c7768d8c]
準備が整ったので、solvers.lpを使って例題を解いてみます。
無事、解として（0.5, 1.5）が求まりました。
sageへの入力:
#pre{{
# 最適問題を解く
sol = solvers.lp(c, G, h)
}}
#pre{{
pcost dcost gap pres dres k/t
0: 2.6471e+00 -7.0588e-01 2e+01 8e-01 2e+00 1e+00
1: 3.0726e+00 2.8437e+00 1e+00 1e-01 2e-01 3e-01
2: 2.4891e+00 2.4808e+00 1e-01 1e-02 2e-02 5e-02
3: 2.4999e+00 2.4998e+00 1e-03 1e-04 2e-04 5e-04
4: 2.5000e+00 2.5000e+00 1e-05 1e-06 2e-06 5e-06
5: 2.5000e+00 2.5000e+00 1e-07 1e-08 2e-08 5e-08
}}
sageへの入力:
#pre{{
print sol['x']
}}
#pre{{
5.0000e-01
1.5000e+00
}}
** 2次最適化問題 [#vbfbffa8]
線形の例題がうまく動いたので、今度は2次最適化問題に進みます。
*** 2次計画問題の例題 [#r7bf22d9]
CVXOPTの本に紹介されている2次計画問題の例題、次のようなものです。
$$
\begin{array}{ll}
\mbox{minimize} & 2 x_1^2 + x_2^2 + x_1 x_2 + x_1 + x_2 \\
\mbox{subject to} & x_1 \geq 0 \\
& x_1 \geq 0 \\
& x_1 + x_2 = 1
\end{array}
$$
*** 2次凸目的関数と制約条件の可視化 [#s3e91899]
線形の時と同様に2次凸目的関数と制約条件を可視化してみます。
コンター図に沿って描画した\( x + y = 1 \)の点線（青）が最も
小さな値を取るのは、(0.25, 0.75)付近であることが見て取れます。
sageへの入力:
#pre{{
# 目的関数の形を見る
var('x y')
f(x, y) = 2*x*x + y*y + x*y + x + y
cntplt = contour_plot(f, (x, 0, 1), (y, 0,1))
lst = [(x, -x + 1) for x in srange(0, 1, 0.05)]
lstplt = list_plot(lst, rgbcolor='blue')
(cntplt + lstplt).show(xmin=0, xmax=1)
}}
&ref(2.png);
*** solvers.cpを使って解く [#z6afb127]
qpのマニュアルには、2次計画問題の標準系として、次の形式で目的関数を表しています。
$$
\begin{array}{ll}
\mbox{minimize} & (1/2)x^TPx + q^Tx \\
\mbox{subject to} & G x \leq h \\
& A x = b
\end{array}
$$
例題にこの式を当てはめると、Gの不等式制約は、変数xが正の値をとるAに等式制約条件を
記述しまと、P, q, G, h, A, bは次のようになります。
$$
\begin{array}{lll}
Q & = & 2 \left( \begin{array}{cc} 2 & 1/2 \\ 1/2 & 1 \end{array} \right) \\
p^T & = & \left( \begin{array}{rr} 1 & 1 \end{array} \right) \\
G & = & \left( \begin{array}{rr} -1 & 0 \\ 0 & -1 \end{array} \right) \\
h & = & \left( \begin{array}{r} 0 \\ 0 \end{array} \right) \\
A & = & \left( \begin{array}{rr} 1 & 1 \end{array} \right) \\
b & = & \left( \begin{array}{r} 0 \end{array} \right)
\end{array}
$$
計算結果は、可視化で期待したとおり、(0.25, 0.75)が最適解となっています。
sageへの入力:
#pre{{
# 2次最適化問題
#
P = 2*_matrix([ [2, .5], [.5, 1] ])
q = _matrix([1.0, 1.0])
G = _matrix([[-1.0,0.0],[0.0,-1.0]])
h = _matrix([0.0,0.0])
A = _matrix([1.0, 1.0], (1,2))
b = _matrix(1.0)
sol=solvers.qp(P, q, G, h, A, b)
}}
#pre{{
pcost dcost gap pres dres
0: 0.0000e+00 0.0000e+00 3e+00 1e+00 0e+00
1: 1.0776e+00 1.3668e+00 6e-01 4e-01 2e-16
2: 1.8460e+00 1.8291e+00 6e-02 2e-02 4e-16
3: 1.8741e+00 1.8681e+00 8e-03 1e-03 3e-16
4: 1.8750e+00 1.8748e+00 2e-04 3e-05 5e-16
5: 1.8750e+00 1.8750e+00 3e-06 2e-07 1e-16
6: 1.8750e+00 1.8750e+00 3e-08 1e-09 6e-16
}}
sageへの入力:
#pre{{
print sol['x']
}}
#pre{{
2.5000e-01
7.5000e-01
}}
** コメント [#jebff62a]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
