为什么scipy.optimize.linprog return 是不满足约束的解呢？
Why does scipy.optimize.linprog return a solution that does not satisfy constraints?
我是做错了什么还是一个错误？
c = np.array([-1., 0., 0., 0., 0., 0., 0., 0., 0.])
A_ub = np.array([[ 1., -724., 911., -551., -555., -896., 478., -80., -293.],
[ 1., 566., 42., 937., 233., 883., 392., -909., 57.],
[ 1., -208., -894., 539., 321., 532., -924., 942., 55.],
[ 1., 857., -859., 83., 462., -265., -971., 826., 482.],
[ 1., 314., -424., 245., -424., 194., -443., -104., -429.],
[ 1., 540., 679., 361., 149., -827., 876., 633., 302.],
[ 0., -1., -0., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -1., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -1., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -1., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -1., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -1., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -1., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -0., -1.],
[ 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
b_ub = np.array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 1., 1., 1., 1., 1., 1., 1.])
A_eq = np.array([[ 0., 1., 1., 1., 1., 1., 1., 1., 1.]])
b_eq = np.array([[ 1.]])
bounds = [(None, None), (None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None)]
soln = scipy.optimize.linprog(c, A_ub, b_ub, A_eq, b_eq, bounds)
print(soln)
print(A_ub.dot(soln['x']) <= b_ub)
print(A_ub.dot(soln['x']) - b_ub)
print(abs(A_eq.dot(soln['x']) - b_eq))](url)
产出
fun: 39.866871820032244
message: 'Optimization terminated successfully.'
nit: 19
slack: array([ 0. , 0. , 0. , 0. , 365.161, 0. ,
0. , 0. , 0. , 0. , 0.167, 0.61 ,
0.115, 0.518, 1. , 1. , 1. , 1. ,
0.833, 0.39 , 0.885, 0.482])
status: 0
success: True
x: array([-39.812, -0.281, -0.625, 0. , -0.031, -0.125, 0.625,
0.312, 0.406])
[ True True False False True False False False True False False True
True True True True True True True True True True]
[-121.5 -358.812 240.125 121.781 -357.781 350.656 0.281 0.625
0. 0.031 0.125 -0.625 -0.312 -0.406 -1.281 -1.625
-1. -1.031 -1.125 -0.375 -0.688 -0.594]
[[ 0.719]]
从A_ub的中间部分可以看出解一定要满足soln['x'][1:] > 0，但是显然不满足！
我是不是做错了什么？
Python 3.6
scipy==0.18.1
是的，它看起来像一个错误。我通常建议人们使用 scipy 的 linprog 的替代品，因为：
- 鲁棒性差
- 性能差（尤其是在大的稀疏问题上）
- 好像不维护了；尽管存在未解决的问题，但进展不大
目前有一个基于内部点的solver in development for scipy, which solves this problem correctly. Of course you can also use the much better alternatives available (e.g. through cvxpy)。
这是一些代码，它显示了 cvxpy 的默认求解器（对于这类问题；ECOS 可以解决更广泛的 class 问题！）ECOS 和尚未合并的 IPM -solver 解决它，而 linprog-simplex 挣扎。请记住，代码不会 运行 在 vanilla-scipy 上，因为缺少 method='interior-point'：
import numpy as np
from scipy.optimize import linprog
c = np.array([-1., 0., 0., 0., 0., 0., 0., 0., 0.])
A_ub = np.array([[ 1., -724., 911., -551., -555., -896., 478., -80., -293.],
[ 1., 566., 42., 937., 233., 883., 392., -909., 57.],
[ 1., -208., -894., 539., 321., 532., -924., 942., 55.],
[ 1., 857., -859., 83., 462., -265., -971., 826., 482.],
[ 1., 314., -424., 245., -424., 194., -443., -104., -429.],
[ 1., 540., 679., 361., 149., -827., 876., 633., 302.],
[ 0., -1., -0., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -1., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -1., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -1., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -1., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -1., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -1., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -0., -1.],
[ 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
b_ub = np.array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 1., 1., 1., 1., 1., 1., 1.])
A_eq = np.array([[ 0., 1., 1., 1., 1., 1., 1., 1., 1.]])
b_eq = np.array([[ 1.]])
bounds = [(None, None), (None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None)]
soln = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method='simplex')
print('linprog: ', soln)
soln = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method='interior-point')
print('linprog IPM: ', soln)
from cvxpy import *
x = Variable(A_eq.shape[1])
constraints = []
constraints.append(A_ub*x <= b_ub)
constraints.append(A_eq*x == b_eq)
objective = Minimize(c*x)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
print(np.round(x.value, 2))
print(problem.value)
输出：
linprog: fun: 39.866871820032244
message: 'Optimization terminated successfully.'
nit: 19
slack: array([ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 3.65161351e+02, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 1.66922982e-01, 6.10104031e-01,
1.14953670e-01, 5.18298274e-01, 1.00000000e+00,
1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
8.33077018e-01, 3.89895969e-01, 8.85046330e-01,
4.81701726e-01])
status: 0
success: True
x: array([ -3.98125000e+01, -2.81250000e-01, -6.25000000e-01,
0.00000000e+00, -3.12500000e-02, -1.25000000e-01,
6.25000000e-01, 3.12500000e-01, 4.06250000e-01])
linprog IPM: con: array([ -3.93646007e-08])
fun: 108.56853369114262
message: 'Optimization terminated successfully.'
nit: 9
slack: array([ 1.23442489e+02, -1.13909794e-05, -7.46066462e-07,
3.12539380e+02, 2.20799190e+02, -1.41898576e-05,
2.48742446e-09, 3.71114180e-01, 1.07937459e-09,
1.33093066e-08, 3.70892271e-01, 2.03637625e-09,
2.57993546e-01, 2.36290348e-08, 9.99999998e-01,
6.28885820e-01, 9.99999999e-01, 9.99999987e-01,
6.29107729e-01, 9.99999998e-01, 7.42006454e-01,
9.99999976e-01])
status: 0
success: True
x: array([ -1.08568534e+02, 2.48742446e-09, 3.71114180e-01,
1.07937459e-09, 1.33093066e-08, 3.70892271e-01,
2.03637625e-09, 2.57993546e-01, 2.36290348e-08])
ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +2.934e+01 -6.533e+02 +2e+03 3e-01 5e-01 1e+00 8e+01 --- --- 1 1 - | - -
1 +9.527e+01 -1.952e+02 +9e+02 1e-01 2e-01 8e+00 4e+01 0.6128 2e-01 1 1 1 | 0 0
2 +9.281e+01 -1.256e+01 +3e+02 4e-02 6e-02 5e+00 1e+01 0.6888 1e-01 1 1 1 | 0 0
3 +1.004e+02 +6.363e+01 +1e+02 2e-02 2e-02 2e+00 5e+00 0.7367 1e-01 1 1 1 | 0 0
4 +1.075e+02 +9.684e+01 +3e+01 5e-03 6e-03 9e-01 1e+00 0.8307 1e-01 1 1 1 | 0 0
5 +1.086e+02 +1.084e+02 +4e-01 6e-05 7e-05 1e-02 2e-02 0.9872 1e-04 1 1 1 | 0 0
6 +1.086e+02 +1.086e+02 +5e-03 7e-07 8e-07 1e-04 2e-04 0.9890 1e-04 1 1 1 | 0 0
7 +1.086e+02 +1.086e+02 +5e-05 8e-09 9e-09 1e-06 2e-06 0.9890 1e-04 1 0 0 | 0 0
8 +1.086e+02 +1.086e+02 +6e-07 8e-11 1e-10 2e-08 3e-08 0.9890 1e-04 1 0 0 | 0 0
OPTIMAL (within feastol=9.6e-11, reltol=5.2e-09, abstol=5.6e-07).
Runtime: 0.000424 seconds.
[[-108.57]
[ -0. ]
[ 0.37]
[ -0. ]
[ 0. ]
[ 0.37]
[ -0. ]
[ 0.26]
[ 0. ]]
108.56853545568384
我是做错了什么还是一个错误？
c = np.array([-1., 0., 0., 0., 0., 0., 0., 0., 0.])
A_ub = np.array([[ 1., -724., 911., -551., -555., -896., 478., -80., -293.],
[ 1., 566., 42., 937., 233., 883., 392., -909., 57.],
[ 1., -208., -894., 539., 321., 532., -924., 942., 55.],
[ 1., 857., -859., 83., 462., -265., -971., 826., 482.],
[ 1., 314., -424., 245., -424., 194., -443., -104., -429.],
[ 1., 540., 679., 361., 149., -827., 876., 633., 302.],
[ 0., -1., -0., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -1., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -1., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -1., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -1., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -1., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -1., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -0., -1.],
[ 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
b_ub = np.array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 1., 1., 1., 1., 1., 1., 1.])
A_eq = np.array([[ 0., 1., 1., 1., 1., 1., 1., 1., 1.]])
b_eq = np.array([[ 1.]])
bounds = [(None, None), (None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None)]
soln = scipy.optimize.linprog(c, A_ub, b_ub, A_eq, b_eq, bounds)
print(soln)
print(A_ub.dot(soln['x']) <= b_ub)
print(A_ub.dot(soln['x']) - b_ub)
print(abs(A_eq.dot(soln['x']) - b_eq))](url)
产出
fun: 39.866871820032244
message: 'Optimization terminated successfully.'
nit: 19
slack: array([ 0. , 0. , 0. , 0. , 365.161, 0. ,
0. , 0. , 0. , 0. , 0.167, 0.61 ,
0.115, 0.518, 1. , 1. , 1. , 1. ,
0.833, 0.39 , 0.885, 0.482])
status: 0
success: True
x: array([-39.812, -0.281, -0.625, 0. , -0.031, -0.125, 0.625,
0.312, 0.406])
[ True True False False True False False False True False False True
True True True True True True True True True True]
[-121.5 -358.812 240.125 121.781 -357.781 350.656 0.281 0.625
0. 0.031 0.125 -0.625 -0.312 -0.406 -1.281 -1.625
-1. -1.031 -1.125 -0.375 -0.688 -0.594]
[[ 0.719]]
从A_ub的中间部分可以看出解一定要满足soln['x'][1:] > 0，但是显然不满足！
我是不是做错了什么？
Python 3.6 scipy==0.18.1
是的，它看起来像一个错误。我通常建议人们使用 scipy 的 linprog 的替代品，因为：
- 鲁棒性差
- 性能差（尤其是在大的稀疏问题上）
- 好像不维护了；尽管存在未解决的问题，但进展不大
目前有一个基于内部点的solver in development for scipy, which solves this problem correctly. Of course you can also use the much better alternatives available (e.g. through cvxpy)。
这是一些代码，它显示了 cvxpy 的默认求解器（对于这类问题；ECOS 可以解决更广泛的 class 问题！）ECOS 和尚未合并的 IPM -solver 解决它，而 linprog-simplex 挣扎。请记住，代码不会 运行 在 vanilla-scipy 上，因为缺少 method='interior-point'：
import numpy as np
from scipy.optimize import linprog
c = np.array([-1., 0., 0., 0., 0., 0., 0., 0., 0.])
A_ub = np.array([[ 1., -724., 911., -551., -555., -896., 478., -80., -293.],
[ 1., 566., 42., 937., 233., 883., 392., -909., 57.],
[ 1., -208., -894., 539., 321., 532., -924., 942., 55.],
[ 1., 857., -859., 83., 462., -265., -971., 826., 482.],
[ 1., 314., -424., 245., -424., 194., -443., -104., -429.],
[ 1., 540., 679., 361., 149., -827., 876., 633., 302.],
[ 0., -1., -0., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -1., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -1., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -1., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -1., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -1., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -1., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -0., -1.],
[ 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
b_ub = np.array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 1., 1., 1., 1., 1., 1., 1.])
A_eq = np.array([[ 0., 1., 1., 1., 1., 1., 1., 1., 1.]])
b_eq = np.array([[ 1.]])
bounds = [(None, None), (None, None), (None, None), (None, None), (None, None),
(None, None), (None, None), (None, None), (None, None)]
soln = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method='simplex')
print('linprog: ', soln)
soln = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method='interior-point')
print('linprog IPM: ', soln)
from cvxpy import *
x = Variable(A_eq.shape[1])
constraints = []
constraints.append(A_ub*x <= b_ub)
constraints.append(A_eq*x == b_eq)
objective = Minimize(c*x)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
print(np.round(x.value, 2))
print(problem.value)
输出：
linprog: fun: 39.866871820032244
message: 'Optimization terminated successfully.'
nit: 19
slack: array([ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 3.65161351e+02, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 1.66922982e-01, 6.10104031e-01,
1.14953670e-01, 5.18298274e-01, 1.00000000e+00,
1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
8.33077018e-01, 3.89895969e-01, 8.85046330e-01,
4.81701726e-01])
status: 0
success: True
x: array([ -3.98125000e+01, -2.81250000e-01, -6.25000000e-01,
0.00000000e+00, -3.12500000e-02, -1.25000000e-01,
6.25000000e-01, 3.12500000e-01, 4.06250000e-01])
linprog IPM: con: array([ -3.93646007e-08])
fun: 108.56853369114262
message: 'Optimization terminated successfully.'
nit: 9
slack: array([ 1.23442489e+02, -1.13909794e-05, -7.46066462e-07,
3.12539380e+02, 2.20799190e+02, -1.41898576e-05,
2.48742446e-09, 3.71114180e-01, 1.07937459e-09,
1.33093066e-08, 3.70892271e-01, 2.03637625e-09,
2.57993546e-01, 2.36290348e-08, 9.99999998e-01,
6.28885820e-01, 9.99999999e-01, 9.99999987e-01,
6.29107729e-01, 9.99999998e-01, 7.42006454e-01,
9.99999976e-01])
status: 0
success: True
x: array([ -1.08568534e+02, 2.48742446e-09, 3.71114180e-01,
1.07937459e-09, 1.33093066e-08, 3.70892271e-01,
2.03637625e-09, 2.57993546e-01, 2.36290348e-08])
ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +2.934e+01 -6.533e+02 +2e+03 3e-01 5e-01 1e+00 8e+01 --- --- 1 1 - | - -
1 +9.527e+01 -1.952e+02 +9e+02 1e-01 2e-01 8e+00 4e+01 0.6128 2e-01 1 1 1 | 0 0
2 +9.281e+01 -1.256e+01 +3e+02 4e-02 6e-02 5e+00 1e+01 0.6888 1e-01 1 1 1 | 0 0
3 +1.004e+02 +6.363e+01 +1e+02 2e-02 2e-02 2e+00 5e+00 0.7367 1e-01 1 1 1 | 0 0
4 +1.075e+02 +9.684e+01 +3e+01 5e-03 6e-03 9e-01 1e+00 0.8307 1e-01 1 1 1 | 0 0
5 +1.086e+02 +1.084e+02 +4e-01 6e-05 7e-05 1e-02 2e-02 0.9872 1e-04 1 1 1 | 0 0
6 +1.086e+02 +1.086e+02 +5e-03 7e-07 8e-07 1e-04 2e-04 0.9890 1e-04 1 1 1 | 0 0
7 +1.086e+02 +1.086e+02 +5e-05 8e-09 9e-09 1e-06 2e-06 0.9890 1e-04 1 0 0 | 0 0
8 +1.086e+02 +1.086e+02 +6e-07 8e-11 1e-10 2e-08 3e-08 0.9890 1e-04 1 0 0 | 0 0
OPTIMAL (within feastol=9.6e-11, reltol=5.2e-09, abstol=5.6e-07).
Runtime: 0.000424 seconds.
[[-108.57]
[ -0. ]
[ 0.37]
[ -0. ]
[ 0. ]
[ 0.37]
[ -0. ]
[ 0.26]
[ 0. ]]
108.56853545568384