SciPy linprog: Linear Programming in Python

linprog in SciPy solves linear programming problems: minimise a linear objective subject to linear constraints.

from scipy.optimize import linprog

result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds)

Three things about it surprise people, and all three are in the first example: it only minimises, every variable is non-negative by default, and constraints must be written as <=.

I solved each of these on SciPy 1.18.1 on Python 3.12.5, and the output is exactly what came back.

How do you maximise with scipy linprog?

A workshop makes chairs and tables, each using wood and labour, and you want the most profit. That’s a maximisation, so negate the objective:

from scipy.optimize import linprog

# A workshop makes chairs and tables.
# Profit: 40 per chair, 55 per table   -> we want the MOST profit
# Wood  : 3 units per chair, 5 per table, 240 available
# Labour: 4 hours per chair, 3 per table, 200 available

# linprog only minimises, so negate the profits to maximise them
c = [-40, -55]

A_ub = [[3, 5],        # wood
        [4, 3]]        # labour
b_ub = [240, 200]

result = linprog(c, A_ub=A_ub, b_ub=b_ub)

print("success :", result.success)
print("chairs  :", round(result.x[0], 2))
print("tables  :", round(result.x[1], 2))
print("profit  :", round(-result.fun, 2))      # negate it back

Output:

success : True
chairs  : 25.45
tables  : 32.73
profit  : 2818.18

The negation trick is the one people forget. linprog has no maximise argument, so minimising -c is how it’s done.

Remember to flip result.fun back too, or you’ll report a loss where you made a profit.

linprog bounds default to (0, None)

If you don’t pass bounds, every variable is assumed to lie between 0 and infinity. That’s right for quantities and wrong for anything that can go negative:

from scipy.optimize import linprog

c = [-1, -2]
A_ub = [[1, 1]]
b_ub = [10]

# default: every variable is assumed to be >= 0
default = linprog(c, A_ub=A_ub, b_ub=b_ub)
print("default bounds (0, None) :", default.x)

# let the first variable go negative, down to -5
negatives = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[(-5, None), (0, None)])
print("bounds [(-5,None), (0,None)]:", negatives.x)

# a per-variable range
capped = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[(2, 5), (0, 4)])
print("bounds [(2,5), (0,4)]    :", capped.x)

Output:

default bounds (0, None) : [ 0. 10.]
bounds [(-5,None), (0,None)]: [-5. 15.]
bounds [(2,5), (0,4)]    : [5. 4.]

Pass one (low, high) tuple per variable, using None for no limit. A single tuple applies to all of them.

The middle run only reaches -5 because the bound allows it. Leave the lower bound at None with nothing else holding the variable down and you’ll get status 3, unbounded, instead.

If your answer keeps coming back as zero where you expected a negative number, this is the first thing to check.

Writing a >= constraint in linprog

There’s no A_lb. Everything goes through A_ub as a <= row, so multiply a >= constraint by -1 on both sides:

from scipy.optimize import linprog

# linprog only understands <=, so a >= constraint has to be flipped by multiplying by -1.
#   2x + 3y >= 12      becomes     -2x - 3y <= -12

c = [5, 4]                       # minimise cost
A_ub = [[-2, -3]]                # the flipped constraint
b_ub = [-12]

result = linprog(c, A_ub=A_ub, b_ub=b_ub)
print("x, y :", result.x.round(3))
print("cost :", round(result.fun, 3))
print("check 2x + 3y =", round(2 * result.x[0] + 3 * result.x[1], 3), ">= 12")

Output:

x, y : [0. 4.]
cost : 16.0
check 2x + 3y = 12.0 >= 12

Equality constraints have their own pair, A_eq and b_eq, and need no flipping.

Getting the signs right is most of the work in setting these problems up, much as it is with root finding.

What linprog returns: x, fun, status and duals

The object you get back carries more than the answer:

from scipy.optimize import linprog

result = linprog([-40, -55], A_ub=[[3, 5], [4, 3]], b_ub=[240, 200])

print("x       :", result.x.round(4))
print("fun     :", round(result.fun, 4))
print("status  :", result.status, "|", result.message.split("(")[0].strip())
print("success :", result.success)
print("nit     :", result.nit)
print("slack   :", result.ineqlin.residual.round(4))     # unused capacity per constraint
print("duals   :", result.ineqlin.marginals.round(4))    # value of one more unit

Output:

x       : [25.4545 32.7273]
fun     : -2818.1818
status  : 0 | Optimization terminated successfully.
success : True
nit     : 2
slack   : [0. 0.]
duals   : [-9.0909 -3.1818]

Attribute	What it holds
x	The optimal values, one per variable
fun	The objective at that point
status	0 solved, 2 infeasible, 3 unbounded
success	True only when status is 0
ineqlin.residual	Slack: how much of each constraint is unused
ineqlin.marginals	Shadow prices: the value of relaxing a constraint

Those marginals are the part most tutorials skip. They tell you which constraint is actually holding you back, which is usually the question behind the question.

Integer variables with linprog integrality

You can’t make 26.7 chairs. Pass integrality=1 per variable and SciPy solves a mixed-integer problem instead:

from scipy.optimize import linprog

# you cannot ship 4.7 crates, so force whole numbers
c = [-40, -55]
A_ub = [[3, 5], [4, 3]]
b_ub = [240, 200]

continuous = linprog(c, A_ub=A_ub, b_ub=b_ub)
integers = linprog(c, A_ub=A_ub, b_ub=b_ub, integrality=[1, 1])

print("continuous:", continuous.x.round(3), "profit", round(-continuous.fun, 2))
print("integer   :", integers.x.round(3), "profit", round(-integers.fun, 2))
print("mip gap   :", integers.mip_gap)

Output:

continuous: [25.455 32.727] profit 2818.18
integer   : [25. 33.] profit 2815.0
mip gap   : 0.0

Rounding the continuous answer yourself is not the same thing and can land outside the constraints. Let the solver do it.

linprog status codes: infeasible and unbounded

Two failures come up constantly, and the status code tells you which you have:

from scipy.optimize import linprog

# infeasible: x <= 1 and x >= 5 at the same time
infeasible = linprog([1], A_ub=[[1], [-1]], b_ub=[1, -5])
print("status", infeasible.status, "->", infeasible.message.split("(")[0].strip())

# unbounded: maximise x with nothing holding it back
unbounded = linprog([-1], A_ub=[[-1]], b_ub=[0], bounds=[(None, None)])
print("status", unbounded.status, "->", unbounded.message.split("(")[0].strip())

print()
print("0 = solved, 1 = iteration limit, 2 = infeasible, 3 = unbounded, 4 = numerical trouble")

Output:

status 2 -> The problem is infeasible.
status 3 -> The problem is unbounded.

0 = solved, 1 = iteration limit, 2 = infeasible, 3 = unbounded, 4 = numerical trouble

Infeasible means your constraints contradict each other. Unbounded almost always means a missing constraint or a bound you meant to set.

Check result.status before using result.x. On a failed solve it is still populated, and it is meaningless.

linprog methods: HiGHS and the deprecated names

HiGHS is the default and has been since SciPy 1.9. You rarely need to choose:

import warnings
from scipy.optimize import linprog

problem = dict(c=[-1, -2], A_ub=[[1, 1], [2, 1]], b_ub=[4, 6])

for method in ("highs", "highs-ds", "highs-ipm", "simplex", "interior-point"):
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        result = linprog(method=method, **problem)
        note = "deprecated" if caught else ""
    print(f"{method:<15} x={result.x.round(3)} status={result.status} {note}")

Output:

highs           x=[0. 4.] status=0 
highs-ds        x=[0. 4.] status=0 
highs-ipm       x=[0. 4.] status=0 
simplex         x=[0. 4.] status=0 deprecated
interior-point  x=[0. 4.] status=0 deprecated

simplex, interior-point and revised simplex still run, but each emits a DeprecationWarning and returns the same answer.

Their warning text has been stale for a while, announcing a removal in SciPy 1.11 that still hasn’t happened. Treat them as legacy names and write highs in new code.

More optimisation and array guides:

• SciPy root finding with brentq and root
• Create a 2D array with NumPy
• NumPy data types
• np.add.at for repeated indices
• Find the maximum value in an array

Frequently asked questions