「内点法で2次計画」の編集履歴（バックアップ）一覧はこちら

「内点法で2次計画」(2011/04/30 (土) 02:48:24) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

非線形計画を解く方法は主に逐次2次計画法と内点法があります。
内点法のほうが新しいアルゴリズムのようですね。
内点法ではバリアパラメータμを徐々に小さくしながら解に近づけていきます。

今回は2次計画問題を内点法を使って解いてみます。
2次計画問題は一般的に次の形で書くことができます。
object 0.5xT*D*x + cT*x
subject A*x - b > 0

参考にした文献は
"Numerical Optimization"
#amazon(0387303030)
という本で、本の中では予測子修正子法を用いています。

予測ステップでまずバリアパラメタμをどれくらいにするか計算して、
修正ステップで予測したμを用いて更新ステップを計算します。

 #!/usr/bin/python
 #coding:utf-8
 from numpy import *
 
 EPS_ZERO = 1.0e-9
 
 # 内点法(Predictor-Collector)による2次計画問題の解法
 # object  min z = 0.5*x^T*G*x+cT*x
 # subject Ax >= b
 def predictor_corrector(x, Gmat, c, Amat, b, tau=0.5, eps=1.0e-3, nloop=30):
    ndim  = Gmat.shape[0] # 次元数
    ncnst = Amat.shape[0] # 制約数
    alldim = ndim+2*ncnst
    zeros_ndim_ncnst  = zeros((ndim, ncnst))
    zeros_ncnst       = zeros((ncnst, ncnst))
    I_ncnst           = identity(ncnst)
    # yとlmdの初期値設定
    f = matrix(vstack([-Gmat*x-c, -Amat*x+b, zeros((ncnst,1))]))
    ExMat = matrix(vstack([hstack([Gmat, zeros_ndim_ncnst, -Amat.T]),
                           hstack([Amat, -I_ncnst, zeros_ncnst]),
                           hstack([zeros_ndim_ncnst.T, zeros_ncnst, zeros_ncnst])]))
    delta_xyl0 = linalg.pinv(ExMat)*f
    y   = matrix([[max([1.0, abs(delta_xyl0[ndim+i,0])])] for i in range(ncnst)])
    lmd = matrix([[max([1.0, abs(delta_xyl0[ndim+ncnst+i,0])])] for i in range(ncnst)])
    rd = Gmat*x - Amat.T*lmd + c
    rp = Amat*x - y - b
    for n in range(nloop):
        err = linalg.norm(rd)**2+linalg.norm(rp)**2+(y.T*lmd)[0,0]
        if err < eps:
            break
        # affine scaling step
        lmd_y = matrix([[lmd[i,0]*y[i,0]] for i in range(ncnst)])
        #         [ Gmat,         0, -Amat.T  ]
        # ExMat = [ Amat,        -I,       0  ]
        #         [    0, diag(lmd),  diag(y) ]
        f = matrix(vstack([-rd, -rp, -lmd_y]))
        ExMat = matrix(vstack([hstack([Gmat, zeros_ndim_ncnst, -Amat.T]),
                               hstack([Amat, -I_ncnst, zeros_ncnst]),
                               hstack([zeros_ndim_ncnst.T, diag([lmd[i,0] for i in range(ncnst)]), diag([y[i,0] for i in range(ncnst)]) ])]))
        delta_xyl_aff = linalg.pinv(ExMat)*f
        mu = (y.T*lmd/ncnst)[0,0]
        # alphaの計算
        y_lmd_cnb = matrix(vstack([y, lmd]))
        yl_dyl    = [-y_lmd_cnb[i,0]/delta_xyl_aff[ndim+i,0] for i in range(2*ncnst) if delta_xyl_aff[ndim+i,0] <= -EPS_ZERO]
        if yl_dyl == [] or  min(yl_dyl) > 1.0:
            alpha_aff = 1.0
        else:
            alpha_aff = min(yl_dyl)
        mu_aff = ((y + alpha_aff*delta_xyl_aff[ndim:(ndim+ncnst),0]).T*(lmd + alpha_aff*delta_xyl_aff[(ndim+ncnst):alldim,0])/ncnst)[0,0]
        # 中心パラメータ
        sig = (mu_aff/mu)**3
        # corrector step
        delta_lmd_y_aff = matrix([[delta_xyl_aff[ndim+i,0]*delta_xyl_aff[ndim+ncnst+i,0]] for i in range(ncnst)])
        sig_mu = matrix(sig*mu*ones((ncnst, 1)))
        f = matrix(vstack([-rd, -rp, -lmd_y-delta_lmd_y_aff+sig_mu]))
        delta_xyl = linalg.pinv(ExMat)*f
        # alpha_pri, alpha_dualの計算
        y_dy = [-tau*y[i,0]/delta_xyl[ndim+i,0] for i in range(ncnst) if delta_xyl[ndim+i,0] <= -EPS_ZERO]
        if y_dy == [] or min(y_dy) > 1.0:
            alpha_pri = 1.0
        else:
            alpha_pri = min(y_dy)
        lmd_dlmd = [-tau*lmd[i,0]/delta_xyl[ndim+ncnst+i,0] for i in range(ncnst) if delta_xyl[ndim+ncnst+i] <= -EPS_ZERO]
        if lmd_dlmd == [] or min(lmd_dlmd) > 1.0:
            alpha_dual = 1.0
        else:
            alpha_dual = min(lmd_dlmd)
        alpha_hat = min([alpha_pri, alpha_dual])
        # x, y, lmdの更新
        x += alpha_hat*delta_xyl[0:ndim,0]
        y += alpha_hat*delta_xyl[ndim:(ndim+ncnst),0]
        lmd += alpha_hat*delta_xyl[(ndim+ncnst):alldim,0]
        rp = (1.0 - alpha_pri)*rp
        rd = (1.0 - alpha_dual)*rd+(alpha_pri - alpha_dual)*Gmat*delta_xyl[0:ndim,0]
        yield x
    #return x, y, lmd
  
 c = matrix([[-2.0],
             [-4.0]])
 Dmat = matrix([[2.0, 0.0],
                [0.0, 2.0]])
 Amat = matrix([[1.0, 1.0],
                [1.0, -1.0],
                [-3.0, 1.0]])
 b = matrix([[0.5],
             [1.0],
             [-3.0]])
 # 描画
 from pylab import *
 
 ax = subplot(111, aspect='equal')
 x = arange(-3.0, 3.01, 0.15)
 y = arange(-3.0, 3.01, 0.15)
 X,Y = meshgrid(x, y)
 t = arange(-3.0, 3.01, 0.15)
 func = lambda x, y : 0.5*(Dmat[0,0]*x**2+Dmat[1,1]*y**2)+c[0,0]*x+c[1,0]*y
 const = [lambda x : -Amat[i,0]/Amat[i,1]*x + b[i,0]/Amat[i,1] for i in range(Amat.shape[0])]
 Z = func(X, Y)
 s = [const[i](t) for i in range(Amat.shape[0])]
 pcolor(X, Y, Z)
 for i in range(Amat.shape[0]):
     ax.plot(t,s[i],'k')
 for x in predictor_corrector(matrix([[0.0],[0.0]]), Dmat, c, Amat, b, tau=0.5, eps=1.0e-3, nloop=30):
     print x
     ax.plot([x[0,0]],[x[1,0]],'ro')
     axis([-3,3,-3,3])
 ax.plot([x[0,0]],[x[1,0]],'yo')
 axis([-3,3,-3,3])
 show()

今回は[[matplotlib>http://matplotlib.sourceforge.net/]]を用いて解の収束する様子を描画してみました。
下の絵は三角形の領域が実行可能領域で、青が濃いところがより最適な解であることを示しています。

#image(http://www7.atwiki.jp/whats-mind?cmd=upload&act=open&pageid=15&file=qp_test.png,width=580,height=460,title=qp_test,http://www7.atwiki.jp/whats-mind?cmd=upload&act=open&pageid=15&file=qp_test.png,center)
こんな感じで実行可能領域外からでもちゃんと解に収束しているのが分かります。
更新ベクトルを求めるのは今回は擬似逆行列で求めましたが、
行列を対称行列に変形して[[共役勾配法]]で求めたほうが速いみたいです。
内点法オープンソースの[[IPOPT>http://www.coin-or.org/Ipopt/]]なんかはさらに疎行列演算とかも使って高速化してるみたいですね。

非線形計画を解く方法は主に逐次2次計画法と内点法があります。
内点法のほうが新しいアルゴリズムのようですね。
内点法ではバリアパラメータμを徐々に小さくしながら解に近づけていきます。

今回は2次計画問題を内点法を使って解いてみます。
2次計画問題は一般的に次の形で書くことができます。
object 0.5x^T*D*x + c^T*x
subject A*x - b > 0

参考にした文献は
"Numerical Optimization"
#amazon(0387303030)
という本で、本の中では予測子修正子法を用いています。

予測ステップでまずバリアパラメタμをどれくらいにするか計算して、
修正ステップで予測したμを用いて更新ステップを計算します。

 #!/usr/bin/python
 #coding:utf-8
 from numpy import *
 
 EPS_ZERO = 1.0e-9
 
 # 内点法(Predictor-Collector)による2次計画問題の解法
 # object  min z = 0.5*x^T*G*x+c^T*x
 # subject Ax >= b
 def predictor_corrector(x, Gmat, c, Amat, b, tau=0.5, eps=1.0e-3, nloop=30):
    ndim  = Gmat.shape[0] # 次元数
    ncnst = Amat.shape[0] # 制約数
    alldim = ndim+2*ncnst
    zeros_ndim_ncnst  = zeros((ndim, ncnst))
    zeros_ncnst       = zeros((ncnst, ncnst))
    I_ncnst           = identity(ncnst)
    # yとlmdの初期値設定
    f = matrix(vstack([-Gmat*x-c, -Amat*x+b, zeros((ncnst,1))]))
    ExMat = matrix(vstack([hstack([Gmat, zeros_ndim_ncnst, -Amat.T]),
                           hstack([Amat, -I_ncnst, zeros_ncnst]),
                           hstack([zeros_ndim_ncnst.T, zeros_ncnst, zeros_ncnst])]))
    delta_xyl0 = linalg.pinv(ExMat)*f
    y   = matrix([[max([1.0, abs(delta_xyl0[ndim+i,0])])] for i in range(ncnst)])
    lmd = matrix([[max([1.0, abs(delta_xyl0[ndim+ncnst+i,0])])] for i in range(ncnst)])
    rd = Gmat*x - Amat.T*lmd + c
    rp = Amat*x - y - b
    for n in range(nloop):
        err = linalg.norm(rd)**2+linalg.norm(rp)**2+(y.T*lmd)[0,0]
        if err < eps:
            break
        # affine scaling step
        lmd_y = matrix([[lmd[i,0]*y[i,0]] for i in range(ncnst)])
        #         [ Gmat,         0, -Amat.T  ]
        # ExMat = [ Amat,        -I,       0  ]
        #         [    0, diag(lmd),  diag(y) ]
        f = matrix(vstack([-rd, -rp, -lmd_y]))
        ExMat = matrix(vstack([hstack([Gmat, zeros_ndim_ncnst, -Amat.T]),
                               hstack([Amat, -I_ncnst, zeros_ncnst]),
                               hstack([zeros_ndim_ncnst.T, diag([lmd[i,0] for i in range(ncnst)]), diag([y[i,0] for i in range(ncnst)]) ])]))
        delta_xyl_aff = linalg.pinv(ExMat)*f
        mu = (y.T*lmd/ncnst)[0,0]
        # alphaの計算
        y_lmd_cnb = matrix(vstack([y, lmd]))
        yl_dyl    = [-y_lmd_cnb[i,0]/delta_xyl_aff[ndim+i,0] for i in range(2*ncnst) if delta_xyl_aff[ndim+i,0] <= -EPS_ZERO]
        if yl_dyl == [] or  min(yl_dyl) > 1.0:
            alpha_aff = 1.0
        else:
            alpha_aff = min(yl_dyl)
        mu_aff = ((y + alpha_aff*delta_xyl_aff[ndim:(ndim+ncnst),0]).T*(lmd + alpha_aff*delta_xyl_aff[(ndim+ncnst):alldim,0])/ncnst)[0,0]
        # 中心パラメータ
        sig = (mu_aff/mu)**3
        # corrector step
        delta_lmd_y_aff = matrix([[delta_xyl_aff[ndim+i,0]*delta_xyl_aff[ndim+ncnst+i,0]] for i in range(ncnst)])
        sig_mu = matrix(sig*mu*ones((ncnst, 1)))
        f = matrix(vstack([-rd, -rp, -lmd_y-delta_lmd_y_aff+sig_mu]))
        delta_xyl = linalg.pinv(ExMat)*f
        # alpha_pri, alpha_dualの計算
        y_dy = [-tau*y[i,0]/delta_xyl[ndim+i,0] for i in range(ncnst) if delta_xyl[ndim+i,0] <= -EPS_ZERO]
        if y_dy == [] or min(y_dy) > 1.0:
            alpha_pri = 1.0
        else:
            alpha_pri = min(y_dy)
        lmd_dlmd = [-tau*lmd[i,0]/delta_xyl[ndim+ncnst+i,0] for i in range(ncnst) if delta_xyl[ndim+ncnst+i] <= -EPS_ZERO]
        if lmd_dlmd == [] or min(lmd_dlmd) > 1.0:
            alpha_dual = 1.0
        else:
            alpha_dual = min(lmd_dlmd)
        alpha_hat = min([alpha_pri, alpha_dual])
        # x, y, lmdの更新
        x += alpha_hat*delta_xyl[0:ndim,0]
        y += alpha_hat*delta_xyl[ndim:(ndim+ncnst),0]
        lmd += alpha_hat*delta_xyl[(ndim+ncnst):alldim,0]
        rp = (1.0 - alpha_pri)*rp
        rd = (1.0 - alpha_dual)*rd+(alpha_pri - alpha_dual)*Gmat*delta_xyl[0:ndim,0]
        yield x
    #return x, y, lmd
  
 c = matrix([[-2.0],
             [-4.0]])
 Dmat = matrix([[2.0, 0.0],
                [0.0, 2.0]])
 Amat = matrix([[1.0, 1.0],
                [1.0, -1.0],
                [-3.0, 1.0]])
 b = matrix([[0.5],
             [1.0],
             [-3.0]])
 # 描画
 from pylab import *
 
 ax = subplot(111, aspect='equal')
 x = arange(-3.0, 3.01, 0.15)
 y = arange(-3.0, 3.01, 0.15)
 X,Y = meshgrid(x, y)
 t = arange(-3.0, 3.01, 0.15)
 func = lambda x, y : 0.5*(Dmat[0,0]*x**2+Dmat[1,1]*y**2)+c[0,0]*x+c[1,0]*y
 const = [lambda x : -Amat[i,0]/Amat[i,1]*x + b[i,0]/Amat[i,1] for i in range(Amat.shape[0])]
 Z = func(X, Y)
 s = [const[i](t) for i in range(Amat.shape[0])]
 pcolor(X, Y, Z)
 for i in range(Amat.shape[0]):
     ax.plot(t,s[i],'k')
 for x in predictor_corrector(matrix([[0.0],[0.0]]), Dmat, c, Amat, b, tau=0.5, eps=1.0e-3, nloop=30):
     print x
     ax.plot([x[0,0]],[x[1,0]],'ro')
     axis([-3,3,-3,3])
 ax.plot([x[0,0]],[x[1,0]],'yo')
 axis([-3,3,-3,3])
 show()

今回は[[matplotlib>http://matplotlib.sourceforge.net/]]を用いて解の収束する様子を描画してみました。
下の絵は三角形の領域が実行可能領域で、青が濃いところがより最適な解であることを示しています。

#image(http://www7.atwiki.jp/whats-mind?cmd=upload&act=open&pageid=15&file=qp_test.png,width=580,height=460,title=qp_test,http://www7.atwiki.jp/whats-mind?cmd=upload&act=open&pageid=15&file=qp_test.png,center)
こんな感じで実行可能領域外からでもちゃんと解に収束しているのが分かります。
更新ベクトルを求めるのは今回は擬似逆行列で求めましたが、
行列を対称行列に変形して[[共役勾配法]]で求めたほうが速いみたいです。
内点法オープンソースの[[IPOPT>http://www.coin-or.org/Ipopt/]]なんかはさらに疎行列演算とかも使って高速化してるみたいですね。