線形計画問題をscipy.optimizationで解こうとすると失敗する

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

min A*x
s.t.
Ceq*x-Deq=0
Cineq*x-Dineq>=0

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

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