4 Ways to Solve Linear Programming in Python
SciPy | PuLP | Pyomo | Google OR-Tools
In my previous article, I have demonstrated how to solve linear programming problems using the graphical method. This article will show you how to solve linear programming problems in Python using four different open-source libraries — Scipy, PuLP, Pyomo, and Google OR-Tools.
Note: These libraries do not solve optimization problems directly. Instead, they allow us to solve optimization problems by interfacing with open-source solvers like Cbc, Glpk, HiGHS and commercial solvers like Gurobi.
The model formulation
Let us solve the following linear programming problem:
SciPy
The first option is SciPy’s optimize.linprog. It is quite easy to use, considering many Python users are familiar with the SciPy library. A plus point is that it interfaces with HiGHS, a high-performance linear programming solver. However, SciPy’s linprog only solves minimization problems. Therefore, we need to provide our objective function in minimization form. To do that, we need to “flip” the sign of the coefficients of our objective function — the equivalent of maximizing 10x1 + 5x2 is minimizing -10x1 - 5x2. Note that the linprog function also does not support lower bounded constraints. Any lower bounded constraints (Ax >= b) have to be reformulated (-Ax <= -b). Let’s look at the code example below:
from scipy.optimize import linprog
# declare the decision variable bounds
x1_bounds = (0, None)
x2_bounds = (0, None)
# declare coefficients of the objective function
c = [-10, -5]
# declare the inequality constraint matrix
A = [[1, 1],
[10, 0],
[0, 5]]
# declare the inequality constraint vector
b = [24, 100, 100]
# solve
results = linprog(c=c, A_ub=A, b_ub=b, bounds=[x1_bounds, x2_bounds], method='highs-ds')
# print results
if results.status == 0: print(f'The solution is optimal.')
print(f'Objective value: z* = {results.fun}')
print(f'Solution: x1* = {results.x[0]}, x2* = {results.x[1]}')The solution is optimal.
Objective value: z* = -170.0
Solution: x1* = 10.0, x2* = 14.0
Hooray! The optimization terminated successfully. Hold on — our initial model is a maximization problem! Don’t forget to “flip” the sign of the objective value, -170, to 170.
PuLP
PuLP is one of my favourites for solving simple linear programming problems. Unlike SciPy, PulP does not require us to provide inputs in the form of matrices or vectors. Instead, we can declare the objective function and constraints explicitly. Besides that, PuLP is more flexible than SciPy in model formulation. It supports both minimization and maximization objective functions. The constraints can also be either lower bounded or upper bounded.
from pulp import *
model = pulp.LpProblem('linear_programming', LpMaximize)
# get solver
solver = getSolver('PULP_CBC_CMD')
# declare decision variables
x1 = LpVariable('x1', lowBound = 0, cat = 'continuous')
x2 = LpVariable('x2', lowBound = 0, cat = 'continuous')
# declare objective
model += 10*x1 + 5*x2
# declare constraints
model += x1 + x2 <= 24
model += 10*x1 <= 100
model += 5*x2 <= 100
# solve
results = model.solve(solver=solver)
# print results
if LpStatus[results] == 'Optimal': print('The solution is optimal.')
print(f'Objective value: z* = {value(model.objective)}')
print(f'Solution: x1* = {value(x1)}, x2* = {value(x2)}')The solution is optimal.
Objective value: z* = 170.0
Solution: x1* = 10.0, x2* = 14.0
Pyomo
Pyomo is a good choice for modelling complex optimization problems. It interfaces a wide range of optimization solvers, not just limited to linear solvers. For simple LP problems, we can declare expressions explicitly like the example below. We can also generate rule-based expressions for the objective function and constraints using Python functions (not covered in this article). One disadvantage of Pyomo is that we need to install the solvers (such as Cbc for linear programming) before running the optimization.
from pyomo.environ import *
model = ConcreteModel('linear_programming')
# declare decision variables
model.x1 = Var(domain = NonNegativeReals)
model.x2 = Var(domain = NonNegativeReals)
# declare objective
model.obj = Objective(expr = 10*model.x1 + 5*model.x2, sense = maximize)
# declare constraints
model.c1 = Constraint(expr = model.x1 + model.x2 <= 24)
model.c2 = Constraint(expr = 10*model.x1 <= 100)
model.c3 = Constraint(expr = 5*model.x2 <= 100)
# solve
results = SolverFactory('cbc').solve(model)
# print results
if results.solver.termination_condition == TerminationCondition.optimal: print('The solution is optimal.')
print(f'Objective value: z* = {value(model.obj)}')
print(f'Solution: x1* = {value(model.x1)}, x2* = {value(model.x2)}')The solution is optimal.
Objective value: z* = 170.0
Solution: x1* = 10.0, x2* = 14.0
Google OR-Tools
Google OR-Tools is an exciting one — besides Python, it also supports other languages like Java and C++. In addition to interfacing with external solvers, Google has also developed in-house solvers such as Glop. Have a look at the example below:
from ortools.linear_solver import pywraplp
# get solver
solver = pywraplp.Solver.CreateSolver('GLOP')
# declare decision variables
x1 = solver.NumVar(0.0, solver.infinity(), 'x1')
x2 = solver.NumVar(0.0, solver.infinity(), 'x2')
# declare objective
solver.Maximize(10*x1 + 5*x2)
# declare constraints
solver.Add(x1 + x2 <= 24)
solver.Add(10*x1 <= 100)
solver.Add(5*x2 <= 100)
# solve
results = solver.Solve()
# print results
if results == pywraplp.Solver.OPTIMAL: print(f'The solution is optimal.')
print(f'Objective value: z* = {solver.Objective().Value()}')
print(f'Solution: x1* = {x1.solution_value()}, x2* = {x2.solution_value()}')The solution is optimal.
Objective value: z* = 170.0
Solution: x1* = 10.0, x2* = 14.0
Recap
This article has given examples of solving linear programming problems in Python using SciPy, PuLP, Pyomo, and Google OR-Tools. In short, SciPy’s linprog and PuLP are meant to solve linear programming problems and are relatively easy to use. Pyomo and Google OR-Tools support a wide range of optimization problems, not limited to LP problems. If you are not a fan of Python — you might prefer Google OR-Tools!
