scipy linear programming - formulation of problem

Question

Friends - Can someone help me formulate a LP problem using scipy in python as below, sorry for this naive ask, I am not able to get started at all with this. I could do this in excel, but finding it difficult in python (am new to this library and couldn't solve) I would be very thankful if someone could help me out please:

This is the data:

This is problem formulated

import pulp as p
import numpy as np

arr = np.array([[0.1167, 2.40,   6.95], [0.1327, 3.44, 15.1], [0.1901, 3.76, 12.7]])
arr = arr.transpose()

# create a problem
Lp_prob = p.LpProblem('Problem', p.LpMinimize)

# create variables
x1 = p.LpVariable("x1", lowBound=0, upBound=np.inf)
x2 = p.LpVariable("x2", lowBound=0, upBound=np.inf)
x3 = p.LpVariable("x3", lowBound=0, upBound=np.inf)


# define problem
Lp_prob += 6.95 * x1 + 15.1 * x2 + 12.7 * x3

# define constraints
Lp_prob += x1 * 0.1167 + x2 * .1327 + x3 * 0.1901 >= 1.95
Lp_prob += x1 * 2.4 + x2 * 3.44 + x3 * 3.76 >= 0
Lp_prob += x1 >= x2
Lp_prob += x1 >= 0
Lp_prob += x2 >= 0
Lp_prob += x3 >= 0

# see the problem created
print(Lp_prob)

status = Lp_prob.solve()

PulpSolverError: Pulp: Error while executing C:\Users\FinanceProfessional\.conda\envs\spyder-env\Lib\site-packages\pulp\apis\..\solverdir\cbc\win\64\cbc.exe

Using scipy

from scipy.optimize import linprog
arr = np.array([[0.1167, 2.40,   6.95], [0.1327, 3.44, 15.1], [0.1901, 3.76, 12.7]])
arr = arr.transpose()

c = arr[-1]
A = [arr[0], arr[1], [1,1,0]]
b = [0.09, 0, 0]
    
x0_bounds = (0, None)
x1_bounds = (0, None)
x2_bounds = (0, None)    
    
result = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds, x2_bounds], method='revised simplex')   

print(result)

con: array([], dtype=float64)
     fun: 0.0
 message: 'Optimization terminated successfully.'
     nit: 0
   slack: array([0.09, 0.  , 0.  ])
  status: 0
 success: True
       x: array([0., 0., 0.])

You may want to see if this helps you: coin-or.github.io/pulp/guides/… — jjramsey
– jjramsey, Commented Feb 23, 2022 at 14:27
In short: You cannot solve that with scipy (yet) as this is a discrete-optimization problem due to multiples of 5 (a MILP) and there is no MILP solver available (yet). — sascha
– sascha, Commented Feb 23, 2022 at 16:01
A lot of libraries can do it, including pulp. MILP solvers (or even LP solvers) are one of the most advanced / complex softwares to develop (one popular one is ~1 million lines of c++) and there is a limited amount of people working on those. Scipy suffered somewhat from it's BSD license and maybe also build-tooling which did imho not allow adding one of the very few candidates in the past. But things are changing due to a new (MIT-lic) solver-project actively developed and MILP support will come. And i bet it's much much better than what Excel provides (which has a bad rep in this regard) — sascha
– sascha, Commented Feb 23, 2022 at 17:17

lot · Accepted Answer · 2022-02-23 18:05:56Z

0

from scipy.optimize import minimize

a1, a2, a3 = 1167,1327,1907
b1,b2,b3 = 24000, 34400, 36000
c1,c2,c3 = 69500,15100,12700

x = [10000,10000,10000] 
res = minimize(
    lambda x: c1*x[0]+c2*x[1]+c3*x[2], #what we want to minimize
    x, 
    constraints = (
        {'type':'eq','fun': lambda x: x[0]*a1-x[1]*a2}, #1st subject
        {'type':'ineq','fun': lambda x: a1*x[0]+a2*x[1]+a3*x[2]-7}, #2st subject
        {'type':'ineq','fun': lambda x: b1*x[0]+b2*x[1]+b3*x[2]-0}, #3st subject
        {'type':'eq','fun': lambda x: x[0]%5+x[1]%5+x[2]%5-0}, # x1 x2 x3 are multiple of 5

                  ),
    bounds = ((0,None),(0,None),(0,None)),
    method='SLSQP',options={'disp': True,'maxiter' : 10000})

print(res)

here the output :

> Optimization terminated successfully    (Exit mode 0)
>             Current function value: 381000000.00006175
>             Iterations: 2
>             Function evaluations: 9
>             Gradient evaluations: 2
>      fun: 381000000.00006175
>      jac: array([69500., 15100., 12700.])  message: 'Optimization terminated successfully'
>     nfev: 9
>      nit: 2
>     njev: 2   status: 0  success: True
>        x: array([    0.,     0., 30000.])

I had to multiplied all value by 10000 to avoid mode 8 as explained here

I hope this is what you needed. However you should try Or-Tools, a CP library powerful and easier than scipy.

edit: answer to comment

here is a link to a google collab as the original poster cannot run this code on his side.

edited Feb 23, 2022 at 18:05

answered Feb 23, 2022 at 14:50

lot

1296 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

sascha Over a year ago

This is ignoring a lot of the solvers assumptions: twice differentiability and probably convexity too (if you are interested in a global-optimum). The original problem is a discrete-optimization problem and throwing it at a continuous NLP solver won't do.

Finance_Professional Over a year ago

thank you for your efforts. I would however, try to find alternative solutions. In reality, I have hundreds of rows and multiple columns to work as constraints. defining constraints in suggested manner would complicate it beyond a level where I would be in control of code. I try to avoid too complicated stuff. Thank you again for your contribution.

Finance_Professional Over a year ago

I did try your code, but it seems to be unsuccessful, maybe some bug at my end. Moreover, I dont think the result produced is as desired, the excel solved it with far lower quantities and cost, than the one produced with this code. thank you anyways.

Finance_Professional Over a year ago

@lot Thanks for your suggestion of Or-Tools - thats working too good for me! hopefully it will land me with desired result, just on last leg to debug some logical mistake in results. Thats incredibly simple to use as well!

Collectives™ on Stack Overflow

scipy linear programming - formulation of problem

1 Answer 1

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related