Linear programming with Scipy fails but quadratic programming succeeds in finding a solution

Question

I am trying to solve a linear programming problem with Python. Ling program fails to find a solution. But quad program works. I do not understand why, and I am not sure my formulation of the program in linprog and quad program are equivalent.

The liner programming problem, my code, and the error message from linprog are below.

Code

k = 6
n = 3
indexes = [1, 2, 5, 6, 7, 9]
V = np.zeros((1, k))
count = 0
for ID in xrange(1, 4):
    ind = count * n + ID
    p = indexes.index(ind)
    V[0, p] = 1
    count += 1

bounds = []
for i in xrange(6):
    bounds.append((0, 1))
bounds = tuple(bounds)

W2 = np.identity(k)
W2 = np.multiply(W2, -1.0)
b2 = np.transpose(np.zeros(k))


V = [[1.0, 0.0, 1.0, 0.0, 0.0, 1.0]]
W1 = [[10.,  0.,  0.,  0., 30.,  0.],
 [ 0., 10., 20.,  0.,  0.,  0.],
 [ 0.,  0.,  0., 20.,  0., 30.]]
W1 = np.array(W1)
W3 = [[1., 1., 0., 0., 0., 0.],
 [0., 0., 1., 1., 0., 0.],
 [0., 0., 0., 0., 1., 1.]]
W2 = np.array(W2)
b1 = [10., 20., 30.]
b3 = [1., 1., 1.,]
EQ = np.vstack([W1, W3]).T
Eb = np.vstack([b1, b3]).T


M = np.identity(k)
P = dot(M.T, M)
q = np.zeros(k)


def quadprog_solve_qp(P, q, W2, b2, W1, b1, W3, b3):
    qp_G = .5 * (P + P.T)   # make sure P is symmetric
    qp_a = -q
    qp_C = -np.vstack([W3, W1, W2]).T
    qp_b = -np.hstack([b3, b1, b2])
    meq = W1.shape[0]
    return quadprog.solve_qp(qp_G, qp_a, qp_C, qp_b, meq)[0]

a = quadprog_solve_qp(P=P, q=q, W2=W2, b2=b2, W1=W1, b1=b1, W3=W3, b3=b3)
print  a
V = [1] * k
res = opt.linprog(c=V, A_eq=EQ, b_eq=Eb, bounds=bounds, options={"disp": True})

Error message from linprog failing

Optimization failed. Unable to find a feasible starting point.

You can install quad program with

sudo pip install quadprog

For examples using quadprog see

https://scaron.info/blog/quadratic-programming-in-python.html

The SImplex linear solver in scipy is known to have problems. Newer versions have an interior point solver which looks more solid. — Erwin Kalvelagen, Sep 03 '18 at 15:18

Thomas Sablik · Accepted Answer · 2018-09-03T13:03:40.107

Your problem has no solution. You have the constraints [0, 0, 0.4, 0, 0, 0] * x = 0.8 and x[2] <= 1 but x[2] has to be 2 to fulfil the first constraint.

Here is an example:

'''
min [1, 1, 1, 1, 1, 1] * x

with [ 0   0   0   0   0   0   ]       [ 0    ]
     [ 0   0   0.4 0   0   0   ] * x = [ 0.4  ]
     [ 0   0   0   0.5 0   0   ]       [ 0.25 ]

and  [ 0   0   0   0   0   0   ]       [ 0     ]
     [ 0   0   0   0   0.4 0   ] * x = [ 0.04  ]
     [ 0   0   0   0   0   0.5 ]       [ 0.125 ]

and  0 <= x[i] <= 1, for 0 <= i < 6

'''

import scipy.optimize as opt
import numpy as np

k = 6
n = 3

W1 = np.zeros((n, k))
W1[1, 2] = 0.4
W1[2, 3] = 0.5
b1 = np.zeros([n, 1])
b1[1] = 0.4
b1[2] = 0.25

W3 = np.zeros((n, k))
W3[1, 4] = 0.4
W3[2, 5] = 0.5
b3 = np.zeros([n, 1])
b3[1] = 0.04
b3[2] = 0.125

c = np.array([1] * k)
A_eq = np.vstack([W1, W3])
b_eq = np.vstack([b1, b3])
bounds = np.array([(0, 1)] * k)
options = {"disp": True}
result = opt.linprog(c=c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, options=options)

print(result.x)

I fixed the dimensions of the vectors and matrices. Also I changed the entries of matrices and vectors, so that a solution can be found. Here is the output:

Optimization terminated successfully.
         Current function value: 1.850000    
         Iterations: 4
[0.   0.   1.   0.5  0.1  0.25]

I didn't try quadprog but if your quadprog function finds a solution, then I'm sure that the problems are not equivalent as there is no solution for the original problem.

Here is another example with other values:

import scipy.optimize as opt
import numpy as np

k = 6
n = 3

W1 = np.array([[0.35855621, 0, 0, 0, 0.85993925, 0 ], 
               [0, 0.35855621, 0.05538291, 0, 0, 0 ], 
               [0, 0, 0, 0.05538291, 0, 0.85993925]])
b1 = np.array([[0.35855621], 
               [0.05538291], 
               [0.85993925]])

c = np.array([1] * k)
A_eq = W1
b_eq = b1
bounds = np.array([(0, 1)] * k)
options = {"disp": True}
result = opt.linprog(c=c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, options=options)

print(result.x)

The output is:

Optimization terminated successfully.
         Current function value: 1.571416    
         Iterations: 3
[0.         0.15446089 0.         0.         0.41695528 1.        ]

I ran into a problem with your solution. It mostly work, but it gives some negative values as solutions sometimes, i.e. if you try with difference W1 and b1 for which solutions exist. — user58925, Sep 03 '18 at 12:48
Here is one example, the solution followed by W1 and b1 Optimization terminated successfully. Current function value: 3.000000 Iterations: 9 [ 2.19916938 -1.19916938 8.7635795 -7.7635795 -0.5 1.5 ] W1 [[0.35855621 0. 0. 0. 0.85993925 0. ] [0. 0.35855621 0.05538291 0. 0. 0. ] [0. 0. 0. 0.05538291 0. 0.85993925]] b1 [0.35855621 0.05538291 0.85993925] — user58925, Sep 03 '18 at 12:55
I tested your input but my output is positive. I added it to my answer. — Thomas Sablik, Sep 03 '18 at 13:04
Thanks. I've been trying with several W1, b1 combinations which have a solution. Solutions are not within the (0,1) bound sometimes. But the problem goes away when I add options = {"bland": True} to your program. — user58925, Sep 03 '18 at 13:10
Can you provide an example where the solution is wrong without `{"bland": True}`? — Thomas Sablik, Sep 04 '18 at 13:28

Linear programming with Scipy fails but quadratic programming succeeds in finding a solution

1 Answers1