Search code examples
pythonoptimizationlinear-programmingquadratic-programmingquadprog

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


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.

enter image description here

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

Solution

  • 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.        ]