線形計画問題に対する内点法 in Python
内点法にもいろいろある。ここでやるのはもちろん一番簡単なやつ。
・主アフィンスケーリング法 - 経営工学専攻 - 東京工業大学
http://www.me.titech.ac.jp/~mizu_lab/text/PDF-IP/IP2A-affine.pdf
を見てPythonでナイーブに実装した。
・行列演算はライブラリ任せ(これは速度面でも好ましい)
・本体はコメントを除けば20行足らず
・例外処理なし
・バグ放置
こんな感じで簡単に線形計画問題の最適解が求まります。
例題
以前の記事で商用ソルバを使って解いた問題をまた解く。
GurobiとPythonで数理最適化 - matsulibの日記
とりあえず前処理として標準形に変形する。
加えて、手続きの出発地点となる初期内点x0を与える必要がある。
下のコードでは x0=[20, 10, 2200, 1400, 670] としている。
ソースコード
#! /usr/bin/env python # -*- coding:utf-8 -*- from pylab import * def affine_scaling(c, A, b, x0, alpha=0.5, eps=10e-4): ''' 次のような標準形の線形計画問題を解く minimize (c^T)*x subject to A*x = b x_i >= 0 (i=1,...,n) Parameters ---------- c : matrix, shape (n, 1); 決定変数の係数ベクトル A : matrix, shape (m, n); 制約式の係数行列 b : matrix, shape (m, 1); 制約式の右辺 x0 : matrix, shape (n, 1); 決定変数の初期値 alpha : float in (0, 1); ステップサイズ(方向ベクトルの重み),デフォルト値:0.5 eps : float > 0; 収束の精度,デフォルト値:10e-4 Returns ------- xk : matrix, shape (n, 1); 最適解 ''' xk = x0 delta_x = [inf] while norm(delta_x) > eps: # Xk : 対角成分が xk の各要素で残りは0の行列 Xk = diagflat(xk.T) # 逆行列の計算を避けるため # yk = (A*Xk*Xk*A.T).I * (A*Xk*Xk*c) を計算する代わりに # (A*Xk*Xk*A.T)*yk = A*Xk*Xk*c を解く yk = solve(A*Xk*Xk*A.T, A*Xk*Xk*c) zk = c - A.T*yk if all(zk==0): break delta_x = -(Xk*Xk*zk)/norm(Xk*zk) xk = xk + alpha*delta_x return xk if __name__ == '__main__': # 問題 c = matrix([[-400], [-300], [0], [0], [0]]) A = matrix([[60, 40, 1, 0, 0], [20, 30, 0, 1, 0], [20, 10, 0, 0, 1]]) b = matrix([[3800], [2100], [1200]]) # 初期値 x1, x2 = 20, 10 x0 = matrix([[x1], [x2], [b[0,0] - A[0,0]*x1 - A[0,1]*x2], [b[1,0] - A[1,0]*x1 - A[1,1]*x2], [b[2,0] - A[2,0]*x1 - A[2,1]*x2]]) x_opt = affine_scaling(c, A, b, x0) for i,xi in enumerate(x_opt): print "x{} : {}".format(i+1, round(xi)) print 'Obj :', round(c.T * x_opt) # Output # x1 : 30.0 # x2 : 50.0 # x3 : 0.0 # x4 : 0.0 # x5 : 100.0 # Obj : -27000.0
matrix*arrayは内積を計算してmatrix型の値を返す。
ちなみにarray*arrayは要素ごとの積となる。
arrayでの内積はdot(array,array)またはarray.dot(array)。
結果
初期値を変えてみたときのイテレーションごとの解を図示する。matplotlibで。
シンプレックス法が隣接する頂点を移動するのに対し、これは確かに内側から最適解に近づいてることがわかる。まさに内点法。
既知のバグ
・パラメータの選び方が難しい。
・epsをやたら小さく設定すると変な答えが出る。逆に大すぎる場合は当然、解が十分収束する前に終了する。
・初期値で自明な実行可能解(x1,x2)=(0,0)を使うと失敗する。