

以下の制約付き線形計画問題を解こうとしています。

[code]

    min A*x
    s.t.
    Ceq*x-Deq=0
    Cineq*x-Dineq>=0
    
[/code]

今回、罰金法を使って制約なし問題に変換して解いています。  
しかし以下のコードを実行すると、"Linear search failed"というエラーが出て最適化に失敗します。L-BFGS-
BとTNCを試しましたが、両方ダメでした。線形計画問題なのでソルバで解けるはずで、何か初歩的なミスをしていると思うのですが、何がいけないのでしょうか？どうしたら最適化できますか？  
将来非線形問題に拡張しようと思っているので、できれば非線形のソルバを使いたいです。  
よろしくお願いします。

[code]

    import numpy as np
    from scipy.optimize import minimize
    
    # objective function: A*x
    A = np.mat([0., 0., 0., 0., 0., 0., 0., 0., 0., -0.5, -0.5, 0., 0., 0., 0., 0., 0., 1.0, 1.0, 0.]);
    
    # matrices for constraints
    Crg = np.mat([ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]);
    Drg = np.mat([0.]).T;
    
    Chm = np.mat([ 1.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]);
    Dhm = np.mat([ 1.]).T;
    
    Cex = np.mat([ 1.,1.,1.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]);
    Dex = np.mat([ 1.]).T;
    
    Cg = np.mat([[ 1.,1.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,0.],
     [ 0.,1.,1.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,0.]]);
    Dg = np.mat([ 0.,1.]).T;
    
    Cdm = np.mat([[ 1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,0., 0.,0.],
     [ 0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,-1., 0.,0.],
     [ 0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,0.],
     [ 0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,-1.]]);
    Ddm = np.mat([ 0.,0.,0.,0.]).T;
    
    Cdr = np.mat([[ 0.,0.,0.,0.,0.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0.,1.,0., 0.,0.],
     [ 0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0.,1., 0.,0.],
     [ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0., 1.,0.],
     [ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0., 0.,1.]]);
    Ddr = np.mat([ 0.,0.,0.,0.]).T;
    
    Cx0 = np.mat([[ 0.,0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,0.],
     [ 0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,0.]]);
    Dx0 = np.mat([1.,0.]).T;
    
    Cz0 = np.mat([[ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0., 0.,0.],
     [ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0., 0.,0.]]);
    Dz0 = np.mat([1.,0.]).T;
    
    Cnc = np.mat([[ 0.,0.,0.,0.,1.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,0.],
     [ 0.,0.,0.,0.,0.,0.,1.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 0.,0.]]);
    Dnc = np.mat([1.,1.]).T;
    
    Cnp = np.mat([[ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,1.,0.,0.,0.,0., 0.,0.],
     [ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,1.,0.,0., 0.,0.]]);
    Dnp = np.mat([1.,1.]).T;
    
    Cfc = np.mat([[ 0.,0.,0.,0.,1.,0.,-1.,0.,0.,-1.,1.,0.,0.,0.,0.,0.,0.,0., 0.,0.],
     [ 0.,0.,0.,0.,0.,1.,0.,-1.,0.,1.,-1.,0.,0.,0.,0.,0.,0.,0., 0.,0.]]);
    Dfc = np.mat([0.,0.]).T;
    
    Cfp = np.mat([[ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,-1.,0.,0.,-1., 1.,0.],
     [ 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,-1.,0.,1.,-1.,0.]]);
    Dfp = np.mat([0.,0.]).T;
    
    # equality constraints Cx-D=0
    Ceq = [Crg,Chm,Cex,Cdm,Cx0,Cz0,Cnc,Cnp,Cfc,Cfp];
    Deq = [Drg,Dhm,Dex,Ddm,Dx0,Dz0,Dnc,Dnp,Dfc,Dfp];
    
    # inequality constraints Cx-D=>0
    Cineq = [Cg,Cdr];
    Dineq = [Dg,Ddr];
    
    bnd = ((0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1));
    
    x0 = [ 0.,1.,0.,0.,1.,0.,1.,0.,0.,0.,0.,0.,1.,0.,0.,1.,0.,1.,0.,0.];
    #xopt = [ 0.,1.,0.,0.,1.,0.,0.,1.,0.,1.,0.,0.,1.,0.,0.,1.,0.,1.,0.,0.];
    
    wpenalty = 100;
    
    def l(x):
      ret = (A*np.mat(x).T)[0,0];
    
      # penalties: lambda*|Cx-D|
      for i,C in enumerate(Ceq):
        c=np.abs(C*np.mat(x).T-Deq[i]);
        lmd = np.mat(np.ones(len(Deq[i])))*wpenalty;
        ret += (lmd*c)[0,0];
      for i,C in enumerate(Cineq):
        c=C*np.mat(x).T-Dineq[i];
        j,k = np.where(c>0);
        c[j,k] = 0.0;
        c = np.abs(c);
        lmd = np.mat(np.ones(len(Dineq[i])))*wpenalty;
        ret += (lmd*c)[0,0];
      print ret;
      return ret;
    
    def dl(x):
      ret = A.A.flatten();
    
      # differetials of penalties: lambda*C*sign(x)
      for i,C in enumerate(Ceq):
        lmd = np.mat(np.ones(len(Deq[i])))*wpenalty;
        ret += (lmd*C).A.flatten()*np.sign(x);
      for i,C in enumerate(Cineq):
        lmd = np.ones(len(Dineq[i]))*wpenalty;
        c=(C*np.mat(x).T-Dineq[i]).T.A.flatten();
        lmd[np.where(c>0)] = 0.0;
        ret += (np.mat(lmd)*C).A.flatten()*np.sign(x);
      return ret;
    
    #res = minimize(l, x0, method='L-BFGS-B',jac=dl,bounds=bnd,options={'disp': True});
    res = minimize(l, x0, method='TNC',jac=dl,bounds=bnd,options={'disp': True});
    print res.fun; # should be 0.5
    print res.x; # should be equal to xopt
    
[/code]


CC BY-SA 3.0
2015-09-27T03:02:54.873
16987
2019-07-01T15:36:09.900
12462
question
python
線形計画問題をscipy.optimizationで解こうとすると失敗する


自己解決しました。  
matlabの最適化関数を使ったら全く同じ入力なのにうまくいきました。  
したがってscipy.optimizeを使うのはおすすめしません。


CC BY-SA 3.0
2015-09-29T12:58:00.407
17066
2015-09-29T12:58:00.407
12493
16987
answer