Solving with Python

Why it matters

Real models have hundreds or millions of variables. Nobody solves those by hand. An OR analyst's day-to-day work is writing a model in code, handing it to a solver, checking the solver's status, and reading the answer. Python has two popular, free routes: scipy.optimize.linprog for models already in matrix form, and PuLP for writing models with named variables and constraints. This lesson shows both, then has you build a small simplex solver in numpy so the solver is no longer a black box.

The running example: a furniture shop

A shop makes desks, tables and chairs, with profits of 60, 30 and 20.

maximize   60d + 30t + 20c
subject to 8d +   6t +    c <= 48
           4d +   2t + 1.5c <= 20
           2d + 1.5t + 0.5c <= 8
           d, t, c >= 0

Route 1: scipy's linprog

linprog always minimizes, and it takes constraints as A_ub @ x <= b_ub (plus optional A_eq @ x == b_eq). To maximize profit, minimize its negative.

import numpy as np
from scipy.optimize import linprog

c = np.array([-60, -30, -20])          # negate: linprog minimizes
A_ub = np.array([[8, 6, 1],
                 [4, 2, 1.5],
                 [2, 1.5, 0.5]])
b_ub = np.array([48, 20, 8])

res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[(0, None)] * 3, method="highs")
print(res.status, res.message)
print(res.x)          # [2. 0. 8.]
print(-res.fun)       # 280.0  (undo the negation)

Always check res.status (0 means an optimal solution was found; other codes mean infeasible, unbounded or a stopped run) before you use res.x. With the HiGHS methods, res.ineqlin.marginals holds the dual values of the <= constraints, which module 4 explains as shadow prices. For a minimization they are zero or negative; because you negated the objective, flip their sign to read them as "extra profit per extra unit of resource".

Route 2: PuLP

PuLP lets you write the model the way you would on paper. It ships with the open-source CBC solver.

import pulp

prob = pulp.LpProblem("furniture", pulp.LpMaximize)
d = pulp.LpVariable("desks", lowBound=0)
t = pulp.LpVariable("tables", lowBound=0)
c = pulp.LpVariable("chairs", lowBound=0)

prob += 60 * d + 30 * t + 20 * c, "profit"
prob += 8 * d + 6 * t + c <= 48, "lumber"
prob += 4 * d + 2 * t + 1.5 * c <= 20, "finishing"
prob += 2 * d + 1.5 * t + 0.5 * c <= 8, "carpentry"

prob.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[prob.status])                     # Optimal
print(d.value(), t.value(), c.value())                # 2.0 0.0 8.0
print(pulp.value(prob.objective))                     # 280.0
print(prob.constraints["finishing"].pi)               # shadow price

Named constraints make reports readable, and changing lowBound=0 to cat="Integer" or cat="Binary" turns the model into an integer program (next module). Both scipy and PuLP are installed with pip install scipy pulp. In this course's browser exercises only numpy is available, so you will write the solver yourself.

Under the hood: the tableau simplex

The numpy version follows the steps from the last lesson.

1. Build the tableau [A | I | b] with the objective row [-c | 0 | 0].
2. Pick the entering column with the most negative objective-row entry. If none is negative, stop: the current corner is optimal.
3. Ratio test: for rows where the column entry is positive, divide the right-hand side by it; the smallest ratio picks the pivot row. If no entry is positive, the LP is unbounded.
4. Pivot: divide the pivot row by the pivot element, then subtract multiples of it from every other row so the column becomes a unit vector.
5. Read the answer: each basic variable equals its row's right-hand side; the rest are zero. The bottom-right cell is the optimal objective.

This simple version assumes every constraint is <= with a non-negative right-hand side, so the slack variables give a starting corner. Real solvers handle >= and = constraints (with a "phase 1"), degeneracy, numerical tolerances and sparse matrices, which is why you use them in practice.

Common mistakes

• Forgetting that linprog minimizes. You get the worst plan, or a sign error on the objective.
• Using linprog's default bounds by accident. They are (0, None), which is right for most models, but free variables need (None, None).
• Ignoring the status. An "answer" from an infeasible or unbounded run is meaningless.
• Mixing units across rows. The solver will not notice that one row is in hours and another in minutes.
• Hard-coding numbers inside the model. Keep data (costs, capacities) in tables or files so the model can be rerun.

Quick recap

• linprog(c, A_ub, b_ub, ...) minimizes; negate c to maximize and negate res.fun back.
• PuLP writes models with named variables and constraints and solves with CBC.
• Check solver status before reading results.
• The tableau simplex is pivots chosen by the most negative cost and the minimum ratio test.

Try it

In the Do tab, put the furniture model in linprog's form and finish the ratio test of a numpy simplex solver.