Solve linear programming problem of otimization with Scipy

Question

I have the following linear prog. problem that I want to solve using Scipy.

Maximize: x0 * c + x1 * d
Such that:
          x0 * a + b * x1 >= 0
          x0 + y0 = 1
          x0, x1 belong [0,1]

I tried this:

from scipy.optimize import linprog

c = [c, d]
A = [[-a, -b], [1, 1]]
b = [0, 1]
x0_bounds = (0, 1)
x1_bounds = (0, 1)
res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])

And the answer, for a, b, c, d = 1.0, -1.0, 1.0, 5.0 should be x0 = 0.5, x1 = 0.5 but I got this:

     con: array([], dtype=float64)
     fun: 1.4471213546122847e-12
 message: 'Optimization terminated successfully.'
     nit: 5
   slack: array([-2.45517088e-13,  1.00000000e+00])
  status: 0
 success: True
       x: array([3.65893190e-14, 2.82106407e-13])

Any ideas why?

Answer 1

The linprog in scipy is incovenient sometimes because:

1. It always solves a minimization problem so if you want to maximize a objective function you have to do a workaround like in this solution to transform it into a minimization problem

2. equations which have >= need to be multiplied by -1 to become <=

3. Creating constraints together such as A_ub A_eq together, they are separates matrices, so create then individually

take a look in the docs they have a nice example as well

linprog is implementing a solver for:

so if you need to have an Inequation of the form A_lb x >= b_lb you need to convert it to a form (expression) <= (values), You can do that by multiplying both sides by -1. So -A_lb x <= -b_lb

FIRST a simple solution to solve your initial problem

from scipy.optimize import linprog
import numpy as np
####defining the minimization problem
##### those are your inputs  following scipy notation
a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
c0, c1= a, b # this is your a and b  10x0 - 3x1 is the obj func to minimize
a_ub0, a_ub1= a, b # this is your a and b in "x0*a + b*x1"  
b0_ub = 0 # this is your constraint in "x0*a + b*x1 >= 0" 
a_eq0, a_eq1= 1, 1 # this is your "1" in "1*x0 + 1*y0"
b0_eq = 1 # this is your 1 constraint in "x0 + x1 = 1" 
x0_bounds = (0, 1) # those are the bounds
x1_bounds = (0, 1)

C = [c0, c1] # you're minimizing so no need to multiply by -1 C
A_ub = [[-a_ub0, -a_ub1]] # but still need to invert the signal here
A_eq = [[a_eq0, a_eq1]]
b_ub = [b0_ub]
b_eq = [b0_eq]
res = linprog(C, A_ub=A_ub, b_ub=b_ub,
                 A_eq=A_eq, b_eq=b_eq,
                 bounds=[x0_bounds, x1_bounds])
minimization_objfunc_output = res['fun']
print(f'min value is {res["fun"]} with solutions x: {res["x"]}')


####defining your original problem
##### those are your inputs  following scipy notation
a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
c0, c1= c, d # this is your c and d 
a_ub0, a_ub1= a, b # this is your a and b in "x0*a + b*x1" 
b0_ub = 0 # this is your constraint in "x0*a + b*x1 >= 0" 
a_eq0, a_eq1= 1, 1 # this is your "1" in "1*x0 + 1*y0"
b0_eq = 1 # this is your 1 constraint in "x0 + x1 = 1" 
a_eq2, a_eq3= a, b # this is your "a" and b in "a*x0 + b*y0"
b1_eq = minimization_objfunc_output # this is your minimization constraint in "a*x0 + b* x1 = result of minimization" 
x0_bounds = (0, 1) # those are the bounds
x1_bounds = (0, 1)

C = [-c0, -c1]
A_ub = [[-a_ub0, -a_ub1]]
A_eq = [[a_eq0, a_eq1],[a_eq2, a_eq3]]
b_ub = [b0_ub]
b_eq = [b0_eq, b1_eq]
res = linprog(C, A_ub=A_ub, b_ub=b_ub,
                 A_eq=A_eq, b_eq=b_eq,
                 bounds=[x0_bounds, x1_bounds])
print(f'max value is {res["fun"]} with solutions x: {res["x"]}')

OUTPUT:

min value is 5.210942788380635e-12 with solutions x: [0.23076923 0.76923077]
max value is -5.897435897450446 with solutions x: [0.23076923 0.76923077]

SECOND a more generic solution to adapt any problem definition to LinProg format

from scipy.optimize import linprog
import numpy as np
##### Ideal generic inputs
def convertInputs_toLinProg(objective, 
                            c, 
                            A_GtE, b_GtE,
                            A_LtE, b_LtE,
                            A_Eq, b_Eq,
                            list_of_bounds):
    A_GtE = -1 * A_GtE
    b_GtE = -1 * b_GtE
    if not (A_GtE is None) and not (A_LtE is None): 
        A_ub = np.vstack([A_GtE,A_GtE])
    elif not (A_GtE is None) and (A_LtE is None):
        A_ub = A_GtE
    elif (A_GtE is None) and not (A_LtE is None):
        A_ub = A_GtE

    if objective == "maximize":
        c = -1*c
    return {
            "c":c, 
            "A_ub": A_ub,
            "A_eq": A_Eq,
            "b_ub": b_GtE,
            "b_eq": b_Eq,
            "bounds": list_of_bounds 
           }

a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
####defining the minimization problem here  
##### those are your inputs  following scipy notation
# objective function:
#    10*x0 - 3*x1 
c = [a,     b]
c = np.array(c)
#                           10*x0 -3*x1 >= 0
A_greaterEq, b_greaterEq =[[a,   b]],   [0]
A_greaterEq, b_greaterEq =np.array(A_greaterEq), np.array(b_greaterEq)
#                      1*x0 + 1*x1 = 1
A_equals, b_equals = [[1,     1]],   [1]
A_equals, b_equals = np.array(A_equals), np.array(b_equals) 
# bounds are exactly how you define them
x0_bounds, x1_bounds = (0, 1), (0, 1)
inputs = convertInputs_toLinProg(objective="minimize",
                                 c=c, 
                                 A_GtE= A_greaterEq, b_GtE= b_greaterEq,
                                 A_LtE= None,       b_LtE= None,
                                 A_Eq=  A_equals,   b_Eq= b_equals,
                                 list_of_bounds=[x0_bounds, x1_bounds])
res = linprog(**inputs)
minimization_objfunc_output = res['fun']
print(f'min value is {res["fun"]} with solutions x: {res["x"]}')
##### generic definition of your problem
####defining your original problem
##### those are your inputs  following scipy notation
# objective function:
#   1*x0 + 5*x1 
a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
c = [c,d]
c = np.array(c)
#                           1*x0 -1*x1 >= 0
A_greaterEq, b_greaterEq =[[a,   b]],   [0]
A_greaterEq, b_greaterEq =np.array(A_greaterEq), np.array(b_greaterEq)
#                      1*x0 + 1*x1 = 1
A_equals, b_equals = [[1,     1],[a,     b]],   [1,minimization_objfunc_output]
A_equals, b_equals = np.array(A_equals), np.array(b_equals) 
# bounds are exactly how you define them
x0_bounds, x1_bounds = (0, 1), (0, 1)

inputs = convertInputs_toLinProg(objective="maximize",
                                 c=c, 
                                 A_GtE= A_greaterEq, b_GtE= b_greaterEq,
                                 A_LtE= None,       b_LtE= None,
                                 A_Eq=  A_equals,   b_Eq= b_equals,
                                 list_of_bounds=[x0_bounds, x1_bounds])

res = linprog(**inputs)
print(f'max value is {res["fun"]} with solutions x: {res["x"]}')

OUTPUT 2

min value is 5.210942788380635e-12 with solutions x: [0.23076923 0.76923077]
max value is 5.897435897450446 with solutions x: [0.23076923 0.76923077]