Search code examples
pythonnumpyoptimizationscipylinear-algebra

Python linprog minimization--simplex method


I'm using scipy.optimize.linprog library to calculate the minimization using the simplex method. I'm working on this problem in my textbook and I'm hoping someone can point me in the right direction because I'm not getting the output I expect. The problem is:

Minimize          w = 10*y1 + 15*y2 + 25*y3
Subject to:       y1 + y2 + y3 >= 1000
                  y1 - 2*y2    >= 0
                            y3 >= 340
with              y1 >= 0, y2 >= 0

The code I wrote for this is:

import numpy as np
import pandas as pd
from scipy.optimize import linprog
A = np.array([
[1, 1, 1],
[1,-2, 0],
[0, 0, 1]])
b = np.array([1000,0,340])
c = np.array([-10,-15,-25])
res = linprog(c, A_ub=A, b_ub=b,
bounds=(0, None))
print('Optimal value:', res.fun, '\nX:', res.x)

Which gives the output:

Optimal value: -18400.0
X: [   0.  660.  340.]

I expect it to be:

Optimal value: -15100.0
X: [   660.  0.  340.]

I can't seem to find consistency with this function but maybe it's the way I'm using it.


Solution

  • You've set up the inputs slightly wrong; see the manual. Specifically, you have a number of sign errors.

    1. Your vector c has the wrong sign; linprog minimizes c x so c should just be the coefficients in w = c x

    2. Your vector b and matrix A have the wrong sign. Their signs should be inverted to switch from your form of constraint f(x) >= const to the desired form for the linprog method, which is a less-than-or-equal, i.e. -f(x) <= - const

    3. You are missing the final two constraints.

    4. Your proposed minimum is < 0, which is obviously impossible as w = 10*x1 + 15*x2 + 25*x3 is always positive with your constraints as x1,x2,x3>=0.

    The correct code reads:

    import numpy as np
    from scipy.optimize import linprog
    
    A = np.array([[-1, -1, -1], [-1,2, 0], [0, 0, -1], [-1, 0, 0], [0, -1, 0]])
    b = np.array([-1000, 0, -340, 0, 0])
    c = np.array([10,15,25])
    
    res = linprog(c, A_ub=A, b_ub=b,bounds=(0, None))
    
    print('Optimal value:', res.fun, '\nX:', res.x)
    # python2
    # ('Optimal value:', 15100.0, '\nX:', array([ 660.,    0.,  340.]))
    # python3
    # Optimal value: 15099.999961403426 
    # X: [6.59999996e+02 1.00009440e-07 3.40000000e+02]