1813-CVXOPTを使って最適解問題を解く.sagews

()

CVXOPTを使って最適解問題を解く

ここでは、sage上でCVXOPTを使って最適解問題を解く方法に説明します。

線形最適解問題

最初に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(x)$を示し、subject to g(w) $\le$ 0は、不等式制約$g(w)を示します。

従ってこの最適化問題は、不等式制約のもとで目的関数が最も小さくなる値を求めることです。

目的関数と制約条件の可視化

CVXOPTを使って解く前に、sageの可視化ツールを使って目的関数と制約条件を表示してみましょう。

黄色の領域が制約条件で囲まれた部分です。$y = 2 x + c$ の直線で、cの値の最も小さいものが、 赤の直線です。

# 制約条件のプロット
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)

CVXOPTに必要なインポート

まずは、CVXOPTに必要なパッケージをインポートします。

# 線形最適解問題
# CVXOPTを使用するのに必要なインポート
# sageではcvxopt.baseからmatrixをインポートする
#
from cvxopt import solvers
from cvxopt.base import matrix as _matrix
RealNumber = float 
Integer=int

solvers.lpを使って解く

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を与えればよいこと が分かります。

    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} $$

# solvers.lpの使い方を調べます（シフト・リターンでヘルプが表示されます）
solvers.lp?

<function lp at 0x9898d30>

# CVXOPTのmaxtrixを作成
# _matrixの要素の与え方がsageと異なることに注意（縦のベクトル要素を並べた形）
c = _matrix([2., 1.])
G = _matrix([[-1., -1., 0., 1.], [1., -1., -1., -2.]])
h = _matrix([1., -2., 0., 4.])

線形最適化例題を解く

準備が整ったので、solvers.lpを使って例題を解いてみます。

無事、解として（0.5, 1.5）が求まりました。

# 最適問題を解く
sol = solvers.lp(c, G, h)

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

print sol['x']

5.0000e-01 1.5000e+00

2次最適化問題

線形の例題がうまく動いたので、今度は2次最適化問題に進みます。

2次計画問題の例題

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次凸目的関数と制約条件の可視化

線形の時と同様に2次凸目的関数と制約条件を可視化してみます。

コンター図に沿って描画した$x + y = 1$の点線（青）が最も 小さな値を取るのは、(0.25, 0.75)付近であることが見て取れます。

# 目的関数の形を見る
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)

solvers.cpを使って解く

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)が最適解となっています。

# 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)

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

print sol['x']

2.5000e-01 7.5000e-01