Why Does Scipy.optimize.linprog Return A Solution That Does Not Satisfy Constraints?

Am I doing something wrong or it is a bug? c = np.array([-1., 0., 0., 0., 0., 0., 0., 0., 0.]) A_ub = np.array([[ 1., -724., 911., -551., -555., -896., 478., -80., -29

Solution 1:

Yes, it looks like a bug. I usually recommend people to use alternatives to scipy's linprog because:

• bad robustness
• bad performance (especially on big sparse problems)
• does not seem to be maintained anymore; not much progress despite open issues

There is currently an interior-point-based 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).

Here is some code, which shows, that cvxpy's default solver (for these kind of problems; ECOS can solve a much broader class of problems!) ECOS and the not yet incorporated IPM-solver solve it, while linprog-simplex struggles. Keep in mind, that the code will not run on vanilla-scipy as method='interior-point' is missing:

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)

Output:

linprog:       fun:39.866871820032244message:'Optimization terminated successfully.'nit:19slack: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:0success:Truex: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.56853369114262message:'Optimization terminated successfully.'nit:9slack: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:0success:Truex: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])ECOS2.0.4-(C)embotechGmbH,ZurichSwitzerland,2012-15. Web:www.embotech.com/ECOSItpcostdcostgappresdresk/tmustepsigmaIR|BT0+2.934e+01-6.533e+02+2e+033e-015e-011e+008e+01------11-|--1+9.527e+01-1.952e+02+9e+021e-012e-018e+004e+010.61282e-01111|002+9.281e+01-1.256e+01+3e+024e-026e-025e+001e+010.68881e-01111|003+1.004e+02+6.363e+01+1e+022e-022e-022e+005e+000.73671e-01111|004+1.075e+02+9.684e+01+3e+015e-036e-039e-011e+000.83071e-01111|005+1.086e+02+1.084e+02+4e-016e-057e-051e-022e-020.98721e-04111|006+1.086e+02+1.086e+02+5e-037e-078e-071e-042e-040.98901e-04111|007+1.086e+02+1.086e+02+5e-058e-099e-091e-062e-060.98901e-04100|008+1.086e+02+1.086e+02+6e-078e-111e-102e-083e-080.98901e-04100|00OPTIMAL(withinfeastol=9.6e-11,reltol=5.2e-09,abstol=5.6e-07).Runtime:0.000424seconds.

[[-108.57]
 [  -0.  ]
 [   0.37]
 [  -0.  ]
 [   0.  ]
 [   0.37]
 [  -0.  ]
 [   0.26]
 [   0.  ]]
108.56853545568384