matsulibの日記

Ingredients as Code

線形計画問題に対する内点法 in Python

内点法にもいろいろある。ここでやるのはもちろん一番簡単なやつ。

・主アフィンスケーリング法 - 経営工学専攻 - 東京工業大学
http://www.me.titech.ac.jp/~mizu_lab/text/PDF-IP/IP2A-affine.pdf
を見てPythonでナイーブに実装した。

・行列演算はライブラリ任せ(これは速度面でも好ましい)
・本体はコメントを除けば20行足らず
・例外処理なし
・バグ放置
こんな感じで簡単に線形計画問題の最適解が求まります。

例題

以前の記事で商用ソルバを使って解いた問題をまた解く。
GurobiとPythonで数理最適化 - matsulibの日記
とりあえず前処理として標準形に変形する。


\begin{eqnarray} 
\mbox{minimize}&~& -400x_1 - 300x_2
\end{eqnarray}

\begin{eqnarray}
\mbox{subject to}
&~& 60x_1 + 40x_2 + x_3 = 3800 \\
&~& 20x_1 + 30x_2 + x_4 = 2100 \\
&~& 20x_1 + 10x_2 + x_5 = 1200 \\
&~& x_i \geq 0, \:\:\: i=1,2,3,4,5
\end{eqnarray}

加えて、手続きの出発地点となる初期内点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で。
f:id:matsulib:20140208225101p:plain:w500
シンプレックス法が隣接する頂点を移動するのに対し、これは確かに内側から最適解に近づいてることがわかる。まさに内点法

既知のバグ

・パラメータの選び方が難しい。
・epsをやたら小さく設定すると変な答えが出る。逆に大すぎる場合は当然、解が十分収束する前に終了する。
・初期値で自明な実行可能解(x1,x2)=(0,0)を使うと失敗する。

もっと良い方法

今回解いたのは線形計画問題だが、内点法はもっと一般的な二次計画問題も解ける。
こちらに素晴らしいPythonのコードがある。
whats-mind @ ウィキ - 内点法で2次計画

二次計画問題線形計画問題も含むので今回の例題をこれで解いてみる。
f:id:matsulib:20140306165632p:plain:w500
等高線図になってて青が濃いほど良い値をとる。最適解に一直線に向かってて良い感じ。