2. Linear Programming

Linear programming (LP) optimizes a linear objective subject to linear equality and inequality constraints. It is the workhorse of Operations Research: it models production, blending, logistics and finance directly, it is the engine inside every integer programming solver, and its duality theory explains what scarce resources are worth. This chapter goes from formulating LPs, through their geometry and the simplex method, to duality, sensitivity analysis and the algorithms used for very large models.

<a id="21-lp-modeling"></a>

2.1 LP Modeling

This subchapter shows how to turn word problems into linear programs and how to rewrite any LP into the standard forms that algorithms expect.

Every LP in this chapter has the general shape

$$ \max \text{ or } \min \ c^\top x \quad \text{s.t.} \quad a_i^\top x \ (\le, =, \ge)\ b_i \ \ (i = 1, \dots, m), \qquad x_j \ge 0,\ x_j \le 0 \text{ or free}. $$

The bakery model from Chapter 1 ($\max 3x_1 + 5x_2$ subject to $x_1 \le 4$, $2x_2 \le 12$, $3x_1 + 2x_2 \le 18$, $x \ge 0$, optimum $(2, 6)$ with profit 36) is the running example.

Standard and Canonical Forms

Why it matters. Users write LPs in whatever shape suits the problem, with maximization or minimization, $\le$, $\ge$ and $=$ constraints, and signed or free variables. Algorithms need one fixed shape. The simplex method works on equalities with nonnegative variables, and duality theory is easiest to state for inequalities. Converting between forms is a mechanical skill you will use constantly.

Theory. Two forms are used throughout this guide:

Form	Objective	Constraints	Variables	Used for
Standard form	$\min c^\top x$	$Ax = b$ (usually $b \ge 0$)	$x \ge 0$	Simplex method, basic solutions
Canonical (inequality) form	$\max c^\top x$	$Ax \le b$	$x \ge 0$	Graphical method, duality

Terminology varies between textbooks. Some call the equality form "canonical" and some use "standard form" with $\max$. Some also say a system is "in canonical form with respect to a basis" when the basic columns form an identity matrix, which is the shape of a simplex tableau. Always check which convention a source uses. The mathematics is the same.

Any LP can be converted into either form without changing its set of optimal solutions:

Original element	Conversion
$\max c^\top x$	$\min (-c)^\top x$, then negate the optimal value
$a^\top x \le b$	$a^\top x + s = b$, slack $s \ge 0$
$a^\top x \ge b$	$a^\top x - e = b$, surplus $e \ge 0$
$a^\top x = b$ (inequality form needed)	$a^\top x \le b$ and $-a^\top x \le -b$
$b_i < 0$ in an equality	Multiply the row by $-1$
$x_j \le 0$	Substitute $x_j = -x_j'$, $x_j' \ge 0$
$x_j$ free	Substitute $x_j = x_j^+ - x_j^-$, both $\ge 0$

Think of the forms as file formats: the content (the feasible region and the optimum) is the same, only the encoding changes.

Example: the bakery in standard form. Add one slack per constraint and flip the objective:

$$ \begin{aligned} \min \quad & -3x_1 - 5x_2 + 0s_1 + 0s_2 + 0s_3 \ \text{s.t.} \quad & x_1 + s_1 = 4 \ & 2x_2 + s_2 = 12 \ & 3x_1 + 2x_2 + s_3 = 18 \ & x_1, x_2, s_1, s_2, s_3 \ge 0 \end{aligned} \qquad A = \begin{bmatrix} 1 & 0 & 1 & 0 & 0 \ 0 & 2 & 0 & 1 & 0 \ 3 & 2 & 0 & 0 & 1 \end{bmatrix} $$

The last three columns form an identity matrix, which gives the simplex method an obvious starting point (Section 2.3).

import numpy as np
from scipy.optimize import linprog

# Canonical form: max c^T x  s.t.  A x <= b, x >= 0
c = np.array([3, 5])
A = np.array([[1, 0], [0, 2], [3, 2]])
b = np.array([4, 12, 18])

def to_standard_form(c, A, b):
    """max c^T x, Ax <= b, x >= 0   ->   min c_std^T z, A_std z = b, z >= 0."""
    m, n = A.shape
    A_std = np.hstack([A, np.eye(m)])          # one slack column per row
    c_std = np.concatenate([-c, np.zeros(m)])  # max -> min, slacks cost nothing
    return c_std, A_std, b

c_std, A_std, b_std = to_standard_form(c, A, b)
res = linprog(c_std, A_eq=A_std, b_eq=b_std, bounds=[(0, None)] * len(c_std))

print("z =", res.x)            # [2. 6. 2. 0. 0.] -> x1, x2, s1, s2, s3
print("max value =", -res.fun) # 36.0, after undoing the sign flip

The standard form solution also reports the slacks: 2 hours of bread oven are left over, and the cake oven and mixer are fully used.

Key takeaways

• Standard form: $\min c^\top x$, $Ax = b$, $x \ge 0$. Canonical form: $\max c^\top x$, $Ax \le b$, $x \ge 0$.
• Every LP can be rewritten in either form by sign flips, slack/surplus variables and variable substitutions.
• Conversions change the number of variables and constraints, but not the optimal decisions.
• Check each textbook's or solver's convention before reading its formulas or its output.

Practice

1. (Beginner) Convert to standard form: $\max 2x_1 - x_2$ s.t. $x_1 + x_2 \ge 3$, $x_1 - 2x_2 \le 4$, $x_1 \ge 0$, $x_2 \le 0$.
2. (Intermediate) Extend to_standard_form to accept a list of senses ("<=", ">=", "=") and add slack or surplus columns only where needed. Test it on Exercise 1.
3. (Interview) "Why do simplex implementations insist on equality constraints and nonnegative variables, while users rarely write models that way?" Hint: think about what a basic solution is and how easy it is to compute one when all constraints are equalities.

Slack and Surplus Variables

Why it matters. Slack and surplus variables do two jobs. Technically, they turn inequalities into equalities so that linear algebra applies. In the business sense, their values tell you how much of each resource is unused, or by how much a requirement is exceeded, which is often the first thing a manager asks about a plan.

Theory.

• For a constraint $a_i^\top x \le b_i$, the slack is $s_i = b_i - a_i^\top x \ge 0$: the unused amount of resource $i$.
• For a constraint $a_i^\top x \ge b_i$, the surplus (or excess) is $e_i = a_i^\top x - b_i \ge 0$: how much the requirement is over-satisfied.
• A constraint is binding when its slack or surplus is $0$.

Analogy: with a monthly budget of 1000 and spending of 850, the slack is 150. If you must walk at least 8000 steps and walk 9200, the surplus is 1200.

Slack and surplus variables have zero cost in the objective, and they have one structural difference:

	Slack ($\le$)	Surplus ($\ge$)
Equation	$a^\top x + s = b$	$a^\top x - e = b$
Coefficient in $A$	$+1$	$-1$
Gives a ready-made starting basis at $x = 0$?	Yes: $s = b \ge 0$	No: $e = -b < 0$ is infeasible
Fix for the starting basis	Not needed	Artificial variables (two-phase or Big-M, Section 2.3)

 Bakery at the optimum x = (2, 6):

 resource       used      available   slack   status
 bread oven     1*2 =  2      4         2     non-binding (idle capacity)
 cake oven      2*6 = 12     12         0     binding     (bottleneck)
 mixer      3*2+2*6 = 18     18         0     binding     (bottleneck)

Example: a diet problem with surplus. Four foods with costs (in EUR) and nutrient contents, and minimum daily requirements:

Food	Cost	Calories	Chocolate (oz)	Sugar (oz)	Fat (oz)
Brownie	0.50	400	3	2	2
Ice cream (scoop)	0.20	200	2	2	4
Cola (bottle)	0.30	150	0	4	1
Cheesecake (slice)	0.80	500	0	4	5
Minimum		500	6	10	8

import numpy as np
from scipy.optimize import linprog

cost = np.array([0.5, 0.2, 0.3, 0.8])
N = np.array([[400, 200, 150, 500],     # calories
              [  3,   2,   0,   0],     # chocolate
              [  2,   2,   4,   4],     # sugar
              [  2,   4,   1,   5]])    # fat
req = np.array([500, 6, 10, 8])

# linprog only accepts <=, so  N x >= req  becomes  -N x <= -req
res = linprog(cost, A_ub=-N, b_ub=-req)

print("servings:", res.x)               # [0. 3. 1. 0.] -> 3 ice creams, 1 cola
print("cost    :", round(res.fun, 2))   # 0.9
# For -N x <= -req, linprog's 'slack' is -req - (-N x) = N x - req = the SURPLUS
for name, e in zip(["calories", "chocolate", "sugar", "fat"], res.slack):
    print(f"{name:<10} surplus = {e:g}")
# calories   surplus = 250    <- 750 calories eaten, 500 required
# chocolate  surplus = 0      <- binding
# sugar      surplus = 0      <- binding
# fat        surplus = 5

Chocolate and sugar are binding, so they drive the cost. The calorie requirement could be raised to 750 without changing the plan.

Key takeaways

• Slack = unused resource in a $\le$ constraint, surplus = excess over a $\ge$ requirement. Both are nonnegative and cost nothing.
• A zero slack or surplus marks a binding constraint, which is a bottleneck of the current plan.
• Slack columns form an identity matrix, which gives the simplex method a free starting basis. Surplus columns do not, so artificial variables are needed.
• Solvers often report slacks directly. Check the sign convention when $\ge$ rows were negated.

Practice

1. (Beginner) For the bakery plan $x = (3, 4)$, compute all three slacks and say which constraints are binding.
2. (Intermediate) In the diet example, add an upper limit of 1000 calories. Write it with a slack variable, re-solve, and check whether the plan changes.
3. (Interview) "A solver reports a slack of $-10^{-9}$ on a constraint. Is the solution infeasible?" Hint: think about floating-point arithmetic and feasibility tolerances.

Product Mix Problems

Why it matters. "Which products should we make, and how many, given limited resources?" is the oldest LP application and still one of the most common. Its structure (maximize profit, $\le$ resource constraints) is the template for capacity planning in factories, data centers and service operations.

Theory. With $n$ products and $m$ resources:

• $x_j$ = units of product $j$ to make, $p_j$ = unit profit (price minus variable cost)
• $a_{ij}$ = amount of resource $i$ used per unit of product $j$, $b_i$ = availability of resource $i$
• optional market limits $u_j$ (maximum sales) and contracts $l_j$ (minimum delivery)

$$ \max \sum_{j=1}^n p_j x_j \quad \text{s.t.} \quad \sum_{j=1}^n a_{ij} x_j \le b_i \ \ (i = 1..m), \qquad l_j \le x_j \le u_j . $$

Modeling checklist:

1. Use contribution margins. Fixed costs (rent, salaries) do not depend on $x$ and do not belong in the objective of an LP.
2. Express every resource in one unit (hours, kg) and match $a_{ij}$ and $b_i$.
3. Put market limits in variable bounds rather than extra constraint rows, because solvers handle bounds more efficiently.

Example: furniture workshop. A workshop makes desks, tables and chairs:

Resource	Desk	Table	Chair	Available
Lumber (board ft)	8	6	1	48
Finishing (h)	4	2	1.5	20
Carpentry (h)	2	1.5	0.5	8
Profit (EUR)	60	30	20

At most 5 tables can be sold.

import numpy as np
from scipy.optimize import linprog

products  = ["desk", "table", "chair"]
profit    = np.array([60, 30, 20])
resources = ["lumber", "finishing", "carpentry"]
A = np.array([[8, 6,   1  ],
              [4, 2,   1.5],
              [2, 1.5, 0.5]])
b = np.array([48, 20, 8])
bounds = [(0, None), (0, 5), (0, None)]   # market limit as a bound, not a row

res = linprog(-profit, A_ub=A, b_ub=b, bounds=bounds)

for name, q in zip(products, res.x):
    print(f"{name:<6}: {q:g}")
print("profit:", -res.fun)
for name, s in zip(resources, res.slack):
    print(f"{name:<9} slack: {s:g}")
# desk  : 2
# table : 0          <- tables are not made, despite the market for 5
# chair : 8
# profit: 280.0
# lumber    slack: 24 <- half the lumber is unused
# finishing slack: 0
# carpentry slack: 0

Tables are not produced. Per finishing and carpentry hour, the two bottlenecks, they earn less than desks and chairs. Section 2.5 shows how to compute exactly how much a table's price would have to rise (by 5) before it becomes worth making.

Key takeaways

• Product mix: maximize contribution margin subject to $\le$ resource limits and market bounds.
• Use variable (marginal) profits only. Fixed costs do not affect the optimal mix.
• Products that look profitable in isolation may be excluded because they use bottleneck resources inefficiently.
• Slack values identify idle resources, and binding constraints identify where extra capacity would pay off.

Practice

1. (Beginner) Write the algebraic model (objective and constraints) for the furniture workshop by hand.
2. (Intermediate) Add a minimum contract of 1 table and re-solve. How much profit does the contract cost?
3. (Interview) "Our product with the highest unit margin is never produced by the optimizer. Is the model broken?" Hint: compare margin per unit of the binding resources, not margin per product unit.

Diet and Blending Problems

Why it matters. Diet and blending models choose a cheapest mix of ingredients that meets quality specifications. They are used for animal feed, fuels, alloys, fertilizers and investment portfolios. They also teach an important modeling trick: turning ratio (percentage) specifications into linear constraints.

Theory.

Diet problem (Stigler, 1945, one of the first problems solved with simplex): choose amounts $x_j$ of foods with cost $c_j$ and nutrient content $a_{ij}$ so that every nutrient $i$ meets its requirement $r_i$:

$$ \min \sum_j c_j x_j \quad \text{s.t.} \quad \sum_j a_{ij} x_j \ge r_i, \qquad x \ge 0 . $$

The previous topic solved a small diet model.

Blending problem: raw materials $j$ with quality $q_j$ (octane, protein %, sulfur %) are mixed into a product whose average quality must reach a specification. Specifications are ratios:

$$ \frac{\sum_j q_j x_j}{\sum_j x_j} \ge Q_{\min}. $$

A ratio is not linear, but since $\sum_j x_j > 0$ you can multiply through by the denominator:

$$ \sum_j q_j x_j \ge Q_{\min} \sum_j x_j \quad \Longleftrightarrow \quad \sum_j (q_j - Q_{\min}), x_j \ge 0 . $$

Analogy: a blend is a weighted average. Ingredients above the target "pay" for ingredients below it.

	Diet	Blending
Variables	Amount of each food	Amount of each raw material (often per product)
Constraints	Nutrient minima/maxima	Quality ratios, supply limits, demand
Typical trick	$\ge$ rows with surplus	Ratio constraint linearized as $\sum (q_j - Q)x_j \ge 0$

Example: gasoline blending. A refinery must produce 5000 barrels of gasoline with average octane at least 89 from three crude streams:

Stream	Octane	Cost / barrel	Available
1	82	40	3000
2	90	55	2000
3	98	70	2000

import numpy as np
from scipy.optimize import linprog

octane = np.array([82, 90, 98])
cost   = np.array([40, 55, 70])
avail  = np.array([3000, 2000, 2000])
target, demand = 89, 5000

# Ratio spec  sum(q x) / sum(x) >= 89   ->   sum((q - 89) x) >= 0   ->   -(q - 89) x <= 0
A_ub = [-(octane - target)]
b_ub = [0]
A_eq = [np.ones(3)]                       # produce exactly 5000 barrels
b_eq = [demand]
bounds = [(0, a) for a in avail]          # supply limits as bounds

res = linprog(cost, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print("barrels:", res.x)                             # [2625.  375. 2000.]
print("cost   :", res.fun)                           # 265625.0
print("octane :", octane @ res.x / res.x.sum())      # 89.0 -> spec is binding

The cheapest blend uses all of the expensive high-octane stream 3 to "lift" as much cheap stream 1 as possible, and the octane specification is met exactly.

Key takeaways

• Diet: minimize cost subject to $\ge$ nutrient requirements. Blending: minimize cost subject to quality specifications on the mix.
• Ratio constraints $\frac{\sum q_j x_j}{\sum x_j} \ge Q$ become linear as $\sum (q_j - Q) x_j \ge 0$.
• Quality specifications are usually binding at the optimum, because overshooting quality costs money.
• Pure-cost diet models produce unpalatable answers. Add variety and upper-bound constraints to make them realistic.

Practice

1. (Beginner) Linearize the constraint "at most 30% of the blend may come from stream 1".
2. (Intermediate) Add a sulfur specification to the gasoline model (sulfur 0.3%, 0.2%, 0.1% for streams 1-3, blend at most 0.18%). Re-solve and identify which specification is binding.
3. (Interview) "A feed company wants its cost-minimizing blend to stay stable when ingredient prices fluctuate slightly. What would you look at?" Hint: objective coefficient ranging (Section 2.5) tells you how far each price can move before the recipe changes.

Production Planning Models

Why it matters. Demand varies over time, but capacity does not. Production planning decides when to produce: build inventory ahead of peaks, use overtime, or subcontract. Multi-period LPs are the backbone of Sales and Operations Planning (S&OP) and Material Requirements Planning (MRP) systems.

Theory. For periods $t = 1, \dots, T$ with demand $d_t$:

• $x_t$ = regular production (capacity $C_t$, cost $c_t$), $y_t$ = overtime production (cost $o_t > c_t$)
• $I_t$ = inventory at the end of period $t$ (holding cost $h$ per unit per period), $I_0$ given

The key element is the inventory balance (flow conservation) constraint:

$$ I_{t-1} + x_t + y_t - d_t = I_t \qquad t = 1, \dots, T $$

$$ \min \sum_{t=1}^{T} \left( c_t x_t + o_t y_t + h I_t \right) \quad \text{s.t. balance}, \quad 0 \le x_t \le C_t, \quad y_t, I_t \ge 0 . $$

Requiring $I_t \ge 0$ means no backorders. Allowing backorders replaces $I_t$ by $I_t^+ - I_t^-$ with a penalty on $I_t^-$, using the free-variable trick from the next topic.

 Inventory flows between periods (a network!):

   x1+y1        x2+y2        x3+y3        x4+y4
     |            |            |            |
     v     I1     v     I2     v     I3     v
 I0->[ Q1 ]------>[ Q2 ]------>[ Q3 ]------>[ Q4 ]--> I4
     |            |            |            |
     v            v            v            v
    d1           d2           d3           d4

Example: a sailboat builder. Demand for the next four quarters is 40, 60, 75, 25 boats, starting inventory is 10, regular capacity is 40 boats per quarter at 400 each, overtime is unlimited at 450 each, and holding costs 20 per boat per quarter.

import numpy as np
from scipy.optimize import linprog

d = [40, 60, 75, 25]
T, I0 = len(d), 10
# Variable layout: [x_1..x_T | y_1..y_T | I_1..I_T]
c = np.r_[np.full(T, 400), np.full(T, 450), np.full(T, 20)]

A_eq = np.zeros((T, 3 * T))
b_eq = np.zeros(T)
for t in range(T):
    A_eq[t, t] = 1            # + x_t
    A_eq[t, T + t] = 1        # + y_t
    A_eq[t, 2 * T + t] = -1   # - I_t
    if t > 0:
        A_eq[t, 2 * T + t - 1] = 1   # + I_{t-1}
    b_eq[t] = d[t] - (I0 if t == 0 else 0)  # known I_0 moves to the right side

bounds = [(0, 40)] * T + [(0, None)] * (2 * T)
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds)

x, y, I = res.x.reshape(3, T)
print("regular  :", x)        # [40. 40. 40. 25.]
print("overtime :", y)        # [ 0. 10. 35.  0.]
print("inventory:", I)        # [10.  0.  0.  0.]
print("cost     :", res.fun)  # 78450.0

Regular capacity is fully used in the first three quarters, the 10 boats in stock are carried into quarter 2, and overtime is used only in the peak quarters 2 and 3. Building a boat early on regular time and holding it (400 + 20 per quarter) would beat overtime (450), but there is no spare regular capacity to do it. Building early on overtime only adds holding cost, so the model never does that.

Key takeaways

• Multi-period models link periods through inventory balance: start + production - demand = end.
• Regular and overtime production are separate variables with separate costs and capacities.
• Building inventory early is worthwhile only when production cost + holding cost < later production cost.
• Production planning LPs are network flow problems over time (Chapter 3).

Practice

1. (Beginner) Write the four inventory balance equations of the sailboat example explicitly.
2. (Intermediate) Require an ending inventory of at least 10 boats ($I_4 \ge 10$) and allow backorders at a penalty of 100 per boat per quarter. Re-solve.
3. (Interview) "How would you extend this model to 50 products sharing the same assembly line?" Hint: index variables by product and period. The capacity constraint then sums over products, while each product keeps its own balance constraint.

Converting Free Variables and Inequalities

Why it matters. Many natural quantities can be negative: net cash flow, deviation from a target, a regression coefficient, a change in inventory. Absolute values and two-sided limits also show up everywhere. Knowing how to fit these into $x \ge 0$ and linear constraints greatly widens the range of problems LP can model.

Theory.

Situation	Transformation	Why it is valid
Free variable $x$	$x = x^+ - x^-$, $x^+, x^- \ge 0$	Any real number is a difference of two nonnegatives
$x \le 0$	$x = -x'$, $x' \ge 0$	Sign flip
$x \ge l$ ($l \ne 0$)	$x = x' + l$, $x' \ge 0$, shift $b$ accordingly	Translation
$x \le u$	Row $x + s = u$, or let the solver handle the bound	Bounded simplex avoids the extra row
$a^\top x = b$ to inequalities	$a^\top x \le b$ and $a^\top x \ge b$	Both directions together force equality
$a^\top x \ge b$ to $\le$	$-a^\top x \le -b$	Multiply by $-1$
$\lvert a^\top x - b\rvert \le t$	$a^\top x - b \le t$ and $-(a^\top x - b) \le t$	Absolute value is the max of two linear functions
$\min \lvert z \rvert$	$z = z^+ - z^-$, minimize $z^+ + z^-$	At optimum, at most one of $z^+, z^-$ is positive
Free variable, alternative	Solve an equality for $x$ and substitute	Removes the variable entirely

Why the split $x = x^+ - x^-$ works. The columns of $x^+$ and $x^-$ are negatives of each other, so they are linearly dependent, and a basic solution (Section 2.2) can never have both basic. Therefore at a vertex at most one of them is nonzero, and the decomposition is "clean". With an objective that minimizes $x^+ + x^-$, the optimum also never wastes value on both.

Example: least absolute deviations (L1) regression. Fit $y \approx a + bx$ minimizing $\sum_i |y_i - a - bx_i|$. The parameters $a, b$ are free, and each residual is split into a positive part $e_i^+$ and a negative part $e_i^-$:

$$ \min \sum_i (e_i^+ + e_i^-) \quad \text{s.t.} \quad a + b x_i + e_i^+ - e_i^- = y_i, \quad a, b \text{ free}, \quad e^\pm \ge 0 . $$

import numpy as np
from scipy.optimize import linprog

x = np.array([1, 2, 3, 4, 5, 6.0])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1, 30.0])   # last point is an outlier
n = len(x)

# Variable layout: [a, b, e+_1..e+_n, e-_1..e-_n]
c = np.r_[0, 0, np.ones(2 * n)]                  # minimize sum of |residuals|
A_eq = np.c_[np.ones(n), x, np.eye(n), -np.eye(n)]
bounds = [(None, None)] * 2 + [(0, None)] * (2 * n)   # a, b free: no sign restriction

res = linprog(c, A_eq=A_eq, b_eq=y, bounds=bounds)
a, b = res.x[:2]
print(f"L1 fit : y = {a:.2f} + {b:.2f} x")       # y = 0.10 + 2.00 x

slope, intercept = np.polyfit(x, y, 1)          # ordinary least squares, for comparison
print(f"L2 fit : y = {intercept:.2f} + {slope:.2f} x")   # y = -5.95 + 4.56 x

The L1 fit follows the five regular points and ignores the outlier, while least squares is pulled far off. Passing bounds=(None, None) asks the solver to do the free-variable conversion internally. Writing $a = a^+ - a^-$ by hand gives the same result.

Key takeaways

• Free variables are split as $x^+ - x^-$, or eliminated by substitution. Nonpositive variables are sign-flipped.
• Equalities become two inequalities, and $\ge$ becomes $\le$ by multiplying by $-1$.
• Absolute values in constraints or in a minimized objective can be linearized. Maximizing $|z|$ cannot, because that problem is non-convex.
• Modern solvers accept free variables and bounds directly, so convert by hand only when an algorithm or textbook requires it.

Practice

1. (Beginner) Rewrite $\min |x_1 - 3| + |x_2 + 1|$ s.t. $x_1 + x_2 = 4$ (with $x_1, x_2$ free) as an LP with only nonnegative variables.
2. (Intermediate) Change the regression code to minimize the maximum absolute residual (Chebyshev fit) instead of the sum. How does the outlier affect this fit?
3. (Interview) "Why can we linearize $\min |z|$ but not $\max |z|$?" Hint: draw the feasible set ${(z, t) : t \ge |z|}$ and the set ${(z, t) : t \le |z|}$, and check which one is convex.

<a id="22-geometry-of-linear-programming"></a>

2.2 Geometry of Linear Programming

This subchapter explains why LP optima sit at corners of the feasible region, and links those corners to the algebraic objects (basic feasible solutions) that the simplex method manipulates.

Graphical Solution Method

Why it matters. With two decision variables you can solve an LP by drawing it. Real models have thousands of variables, but the picture shows almost everything that matters about LP: the feasible region is a polygon, the objective is a family of parallel lines, and the optimum is the last point those lines touch. Every concept later in this chapter (vertices, degeneracy, unboundedness, sensitivity) has a picture here first.

Theory. For $\max c_1x_1 + c_2x_2$:

1. Draw each constraint as its boundary line $a_{i1}x_1 + a_{i2}x_2 = b_i$, and keep the side that satisfies the inequality (test the origin when it is not on the line).
2. Intersect all halfspaces, including $x_1, x_2 \ge 0$, to obtain the feasible region.
3. Draw an iso-profit line $c_1x_1 + c_2x_2 = z$ for some value $z$. Every point on it has the same objective value.
4. Slide the line parallel to itself in the improving direction, which is the gradient $c = (c_1, c_2)$, until it is about to leave the region.
5. Read off the last point(s) touched. That point is optimal. Solve the two binding constraint equations to get exact coordinates.

Analogy: the iso-profit line is a ruler laid on the map. You push it across the map in the direction of $c$, and the last corner of the region the ruler touches is the optimum.

 Bakery: max 3 x1 + 5 x2

  x2
  9 +\
    |  \  z = 45 (infeasible: line misses the region)
  7 +    \.
    |  .   \
  6 *=========* (2,6)   <- last point touched: z* = 3*2 + 5*6 = 36
    |#########|\ .
    |##########\  \ .     iso-profit lines 3x1 + 5x2 = z
  3 +###########* (4,3)    have slope -3/5; c = (3,5) points up-right
    |###########|
    |###########|   z = 15
    *-----------*--------> x1
  (0,0)       (4,0)

The direction of $c$ is perpendicular to the iso-profit lines. Changing the objective coefficients rotates the lines and can move the optimum to a different corner, which is the topic of sensitivity analysis (Section 2.5).

import numpy as np
import matplotlib
matplotlib.use("Agg")                      # render to a file, no window needed
import matplotlib.pyplot as plt

x1 = np.linspace(0, 7, 400)
fig, ax = plt.subplots(figsize=(6, 5))

# Constraint boundaries
ax.axvline(4, color="tab:blue", label="x1 <= 4 (bread oven)")
ax.axhline(6, color="tab:orange", label="2 x2 <= 12 (cake oven)")
ax.plot(x1, (18 - 3 * x1) / 2, color="tab:green", label="3x1 + 2x2 <= 18 (mixer)")

# Feasible region = polygon through the vertices, listed in order
verts = np.array([[0, 0], [4, 0], [4, 3], [2, 6], [0, 6]])
ax.fill(verts[:, 0], verts[:, 1], alpha=0.3, color="gray", label="feasible region")

# Iso-profit lines 3 x1 + 5 x2 = z for increasing z
for z in (15, 30, 36):
    ax.plot(x1, (z - 3 * x1) / 5, "k--", lw=0.8)
    ax.annotate(f"z={z}", xy=(0.1, (z - 0.3) / 5))

ax.plot(2, 6, "ro")
ax.annotate("optimum (2, 6)", xy=(2.1, 6.2))
ax.set(xlim=(0, 7), ylim=(0, 9), xlabel="x1 (bread)", ylabel="x2 (cake)")
ax.legend(loc="upper right", fontsize=8)
fig.savefig("bakery_graphical.png", dpi=120)
print("saved bakery_graphical.png")

Possible outcome	Picture
Unique optimum	The line last touches a single corner
Multiple optima	The line is parallel to an edge and last touches the whole edge
Unbounded	The region extends forever in the direction of $c$
Infeasible	The halfspaces have no common point, so there is nothing to draw

Key takeaways

• Plot the constraints, shade the feasible region, slide an iso-profit line in the direction of $c$.
• The optimum is where the line last touches the region, always at a corner and possibly along a whole edge.
• The objective gradient $c$ is perpendicular to the iso-profit lines and points toward improvement.
• The method only works in 2 (or, with effort, 3) dimensions, but its geometry carries over to any dimension.

Practice

1. (Beginner) Solve graphically: $\max x_1 + x_2$ s.t. $x_1 + 2x_2 \le 8$, $3x_1 + x_2 \le 9$, $x \ge 0$.
2. (Intermediate) For the bakery, find an objective $c$ for which every point on the edge between $(2, 6)$ and $(4, 3)$ is optimal. Verify with the plotting code by drawing one iso-profit line through both points.
3. (Interview) "Explain to a non-technical manager why the optimal plan is always at a 'corner'." Hint: if a point is in the middle of an edge, moving along the edge in one of the two directions does not make the objective worse, because the objective is linear.

Extreme Points and Vertices

Why it matters. The fact that "optima are at corners" only helps if we can define, recognize and enumerate corners precisely in $n$ dimensions, where no picture exists. Three different definitions, one geometric, one optimization-based and one algebraic, turn out to describe the same points. That equivalence is what lets the simplex method work with matrices while thinking in geometry.

Theory. Let $P = {x \in \mathbb{R}^n : Ax \le b}$ be a polyhedron (Chapter 1.2).

Definition	Statement	Intuition
Extreme point	$x \in P$ is not $\lambda y + (1-\lambda) z$ for distinct $y, z \in P$, $0 < \lambda < 1$	Not in the middle of any segment inside $P$
Vertex	Some $c$ makes $x$ the unique maximizer of $c^\top x$ over $P$	Some objective picks out exactly this point
Basic feasible solution	$x \in P$, and $n$ linearly independent constraints are active at $x$	Exactly determined by its tight constraints

Theorem. For a nonempty polyhedron, the three notions coincide: extreme point $\iff$ vertex $\iff$ basic feasible solution.

Two further facts:

• Adjacency. Two vertices are adjacent if they share $n - 1$ linearly independent active constraints. The segment between them is an edge of $P$. The simplex method moves along edges between adjacent vertices.
• Counting. A polyhedron with $m$ constraints in $\mathbb{R}^n$ has at most $\binom{m}{n}$ vertices, and the count can be exponential. The unit cube $[0,1]^n$ has only $2n$ constraints but $2^n$ vertices.
• Existence. A nonempty polyhedron has at least one vertex if and only if it contains no full line. Every LP in standard form ($x \ge 0$) satisfies this.

  A point and its active constraints (n = 2):

        vertex: 2 independent       edge point: 1 active   interior: 0 active
        constraints active
            \  /                         |                     . . .
             \/  <- (2,6)                * <- (3, 4.5)         . * .  <- (1,1)
             /\                          |                     . . .
     unique optimum of some c     midpoint of 2 vertices    midpoint of many

import numpy as np

# Bakery region as A x <= b, including -x <= 0 for nonnegativity
A = np.array([[1, 0], [0, 2], [3, 2], [-1, 0], [0, -1]], dtype=float)
b = np.array([4, 12, 18, 0, 0], dtype=float)
names = ["bread", "cake", "mixer", "x1>=0", "x2>=0"]

def classify(x, tol=1e-9):
    x = np.asarray(x, dtype=float)
    if np.any(A @ x > b + tol):
        return "infeasible", []
    active = np.where(np.abs(A @ x - b) <= tol)[0]           # tight constraints
    rank = np.linalg.matrix_rank(A[active]) if len(active) else 0
    kind = "vertex (BFS)" if rank == len(x) else f"not a vertex (rank {rank})"
    return kind, [names[i] for i in active]

for p in [(2, 6), (4, 3), (3, 4.5), (1, 1), (5, 5)]:
    kind, act = classify(p)
    print(f"{str(p):<10} {kind:<22} active: {act}")
# (2, 6)     vertex (BFS)           active: ['cake', 'mixer']
# (4, 3)     vertex (BFS)           active: ['bread', 'mixer']
# (3, 4.5)   not a vertex (rank 1)  active: ['mixer']        <- on the edge
# (1, 1)     not a vertex (rank 0)  active: []               <- interior
# (5, 5)     infeasible             active: []

The vertices $(2, 6)$ and $(4, 3)$ share one active constraint (the mixer), which is $n - 1 = 1$, so they are adjacent, and the point $(3, 4.5)$ lies on the edge between them.

Key takeaways

• Extreme point (geometric), vertex (unique optimizer) and basic feasible solution (algebraic) are the same thing.
• In $\mathbb{R}^n$, a vertex has $n$ linearly independent active constraints. Fewer than $n$ means the point lies on an edge, face or in the interior.
• Adjacent vertices share $n - 1$ independent active constraints and are joined by an edge.
• The number of vertices can grow exponentially, so explicit enumeration does not scale.

Practice

1. (Beginner) Use classify to check the points $(0, 6)$, $(4, 0)$ and $(2, 3)$.
2. (Intermediate) Write a function that, given a vertex, returns all its adjacent vertices by dropping one active constraint at a time and moving along the resulting edge until a new constraint becomes tight.
3. (Interview) "How many vertices can a 3D polytope defined by 6 inequalities have, and can you give an example that reaches the maximum?" Hint: the cube has 6 faces. Compare its vertex count with $\binom{6}{3}$ and explain why most triples of planes do not meet at a feasible point.

Basic Feasible Solutions

Why it matters. Computers cannot "look at corners". They solve linear systems. Basic feasible solutions are the algebraic version of vertices: pick a set of columns, solve a square system, check signs. Every simplex iteration is exactly this computation, so understanding bases is the key to reading simplex output, sensitivity reports and solver logs.

Theory. Take an LP in standard form, $Ax = b$, $x \ge 0$, with $A \in \mathbb{R}^{m \times n}$ of full row rank $m$ ($m < n$).

1. Choose $m$ columns of $A$ that are linearly independent. Their indices form a basis $\mathcal{B}$, and their matrix $B$ is invertible.
2. Set the other $n - m$ variables, the nonbasic ones $x_N$, to zero.
3. Solve $Bx_B = b$, so $x_B = B^{-1}b$. The result is a basic solution.
4. If additionally $x_B \ge 0$, it is a basic feasible solution (BFS).

$$ Ax = Bx_B + Nx_N = b, \qquad x_N = 0 \ \Rightarrow \ x_B = B^{-1} b . $$

Connection to geometry: setting $n - m$ variables to zero makes $n - m$ of the bounds $x_j \ge 0$ active, which together with the $m$ equalities gives $n$ active constraints, exactly the vertex condition from the previous topic.

Term	Meaning
Basic variables $x_B$	The $m$ variables solved for. They may be positive or zero
Nonbasic variables $x_N$	Fixed at zero
Basic solution	Any $x$ obtained this way (may violate $x \ge 0$)
BFS	Basic solution with $x \ge 0$, which corresponds to a vertex
Degenerate BFS	A BFS with some basic variable equal to 0

There are at most $\binom{n}{m}$ bases, so there are finitely many BFSs. Different bases can, however, give the same BFS (degeneracy).

Example: all bases of the bakery. In standard form there are $n = 5$ variables ($x_1, x_2, s_1, s_2, s_3$) and $m = 3$ equations, so at most $\binom{5}{3} = 10$ bases:

import numpy as np
from itertools import combinations

A = np.array([[1, 0, 1, 0, 0],
              [0, 2, 0, 1, 0],
              [3, 2, 0, 0, 1]], dtype=float)
b = np.array([4, 12, 18], dtype=float)
profit = np.array([3, 5, 0, 0, 0])
names = ["x1", "x2", "s1", "s2", "s3"]

for basis in combinations(range(5), 3):
    B = A[:, basis]
    label = "{" + ",".join(names[j] for j in basis) + "}"
    if abs(np.linalg.det(B)) < 1e-12:            # columns dependent: not a basis
        print(f"{label:<12} singular")
        continue
    x = np.zeros(5)
    x[list(basis)] = np.linalg.solve(B, b)       # x_B = B^-1 b, x_N = 0
    status = "BFS" if np.all(x >= -1e-12) else "infeasible"
    print(f"{label:<12} (x1,x2)=({x[0]:g},{x[1]:g})  {status:<10} z={profit @ x:g}")
# {x1,x2,s1}   (x1,x2)=(2,6)  BFS        z=36   <- optimal vertex
# {x1,x2,s2}   (x1,x2)=(4,3)  BFS        z=27
# {x1,x2,s3}   (x1,x2)=(4,6)  infeasible z=42   <- s3 = -6: outside the region
# {x1,s1,s2}   (x1,x2)=(6,0)  infeasible z=18
# {x1,s1,s3}   singular
# {x1,s2,s3}   (x1,x2)=(4,0)  BFS        z=12
# {x2,s1,s2}   (x1,x2)=(0,9)  infeasible z=45
# {x2,s1,s3}   (x1,x2)=(0,6)  BFS        z=30
# {x2,s2,s3}   singular
# {s1,s2,s3}   (x1,x2)=(0,0)  BFS        z=0    <- the all-slack starting basis

The 5 BFSs are exactly the 5 vertices drawn in Chapter 1.2. The infeasible basic solutions are intersections of constraint lines outside the region. Some of them have a higher objective than 36, which is why feasibility must be checked.

Key takeaways

• A basis is a set of $m$ linearly independent columns. Nonbasic variables are 0 and $x_B = B^{-1}b$.
• A basic solution is feasible (a BFS) iff $x_B \ge 0$, and BFSs correspond one-to-one with vertices (up to degeneracy).
• There are at most $\binom{n}{m}$ bases, which is finite but astronomically large for real models.
• Infeasible basic solutions can have better objective values than the optimum, so feasibility is not optional.

Practice

1. (Beginner) For the basis ${x_2, s_1, s_3}$, solve $Bx_B = b$ by hand and confirm the vertex $(0, 6)$.
2. (Intermediate) Explain geometrically why the basis ${x_1, s_1, s_3}$ is singular. Which two constraint lines would it try to intersect?
3. (Interview) "A model has 1000 constraints and 5000 variables. Why can't we just enumerate all bases?" Hint: estimate $\binom{5000}{1000}$ using logarithms, then compare with the number of atoms in the universe ($\approx 10^{80}$).

Fundamental Theorem of LP

Why it matters. This theorem turns an optimization over infinitely many points into a search over finitely many basic feasible solutions. It is the reason the simplex method can stop after a finite number of steps, and the reason solvers return corner solutions.

Theory. For an LP in standard form $\min{c^\top x : Ax = b,\ x \ge 0}$ with $\operatorname{rank}(A) = m$:

1. If there is a feasible solution, there is a basic feasible solution.
2. If there is an optimal solution, there is an optimal basic feasible solution.

Combined with the three-way classification from Chapter 1.3, every LP is exactly one of: infeasible, unbounded, or optimal with an optimal BFS.

Proof idea ("purification"). Take any feasible (or optimal) $x$. Let $S = {j : x_j > 0}$ be its support.

• If the columns $A_S$ are linearly independent, $x$ is already basic.
• Otherwise there is $d \ne 0$ with $d_j = 0$ for $j \notin S$ and $Ad = 0$. Then $x + td$ still satisfies $Ax = b$ for any $t$. Choose the sign of $d$ so that $c^\top d \le 0$, so the objective does not get worse. Increase $t$ until the first component $x_j + t d_j$ hits $0$. (If no component decreases and $c^\top d < 0$, the LP is unbounded.)
• The new point is still feasible, at least as good, and has a strictly smaller support. Repeat. Since the support shrinks every time, this ends at a BFS after at most $n$ steps.

  Purification in the bakery (maximize, so the objective must not decrease):

  x2
  6 *  <------------ step 2: slide up the x2 axis until cake oven is tight -> (0,6) BFS, z=30
    |
    |
  2.65 *  <--------- step 1: slide along a null-space direction until x1 = 0, z=13.3
    |   \
  1 +    * (1,1) start, z=8, interior: support has 5 positive variables
    +----+-----------> x1

import numpy as np
from scipy.linalg import null_space

# Bakery in standard form, written as a minimization of -profit
A = np.array([[1, 0, 1, 0, 0],
              [0, 2, 0, 1, 0],
              [3, 2, 0, 0, 1]], dtype=float)
c = np.array([-3, -5, 0, 0, 0], dtype=float)

x = np.array([1, 1, 3, 10, 13], dtype=float)    # feasible interior point (x1, x2) = (1, 1)
assert np.allclose(A @ x, [4, 12, 18])

while True:
    S = np.where(x > 1e-9)[0]                   # support: positive variables
    N = null_space(A[:, S])                     # directions with A_S d_S = 0
    if N.shape[1] == 0:                         # columns independent -> x is a BFS
        break
    d = np.zeros(len(x))
    d[S] = N[:, 0]
    if c @ d > 0:                               # never move uphill in cost
        d = -d
    if not np.any(d < -1e-12):                  # nothing ever hits zero along d
        d = -d                                  # (c @ d == 0 here, else unbounded)
    shrinking = d < -1e-12
    t = np.min(-x[shrinking] / d[shrinking])    # step until first variable hits 0
    x = x + t * d
    x[np.abs(x) < 1e-9] = 0
    print(f"(x1, x2) = ({x[0]:.3f}, {x[1]:.3f})   profit = {-(c @ x):.2f}   "
          f"support = {np.count_nonzero(x)}")
# (x1, x2) = (0.000, 2.653)   profit = 13.26   support = 4
# (x1, x2) = (0.000, 6.000)   profit = 30.00   support = 3   <- a BFS, never worse than the start

The vertex reached depends on the directions chosen (the null space basis returned by your SciPy version may differ). The guarantee is only that the BFS is at least as good as the starting point. Starting from an optimal point, purification returns an optimal BFS.

Key takeaways

• If an LP is feasible, it has a BFS. If it has an optimum, some BFS is optimal.
• Every LP is infeasible, unbounded, or has an optimal BFS. There is no fourth case.
• The proof moves along null-space directions that keep $Ax = b$ and do not worsen the objective, shrinking the support each time.
• The search space reduces from infinitely many points to at most $\binom{n}{m}$ BFSs. The simplex method searches them intelligently.

Practice

1. (Beginner) Give an LP whose optimum is attained at a non-basic point as well as at a BFS (hint: multiple optima), and identify both points.
2. (Intermediate) Modify the purification code to start from the optimal non-vertex point of your Exercise 1 and verify it ends at an optimal BFS.
3. (Interview) "Does the fundamental theorem hold for $\max{x_1 : x_1 + x_2 \le 1}$ with $x_1, x_2$ free?" Hint: the feasible region contains a line. Does it have any vertices at all?

Unbounded and Infeasible Problems

Why it matters. In practice, "infeasible" and "unbounded" almost always mean a modeling or data error: a missing constraint, a wrong sign, a typo in a capacity. Recognizing them quickly, and knowing that each comes with a certificate that proves it, saves hours of debugging.

Theory.

Infeasible: the feasible region is empty. In 2D, the halfspaces have no common point. A certificate of infeasibility is a nonnegative combination of the constraints that produces a contradiction such as $0 \le -1$ (Farkas' lemma, Section 2.4).

Unbounded: the region is nonempty and contains a ray $x_0 + t d$, $t \ge 0$, along which the objective improves without limit. For $\max{c^\top x : Ax \le b}$, the certificate is a feasible $x_0$ and a direction $d$ with $Ad \le 0$ and $c^\top d > 0$.

An unbounded region does not imply an unbounded problem:

	Bounded objective	Unbounded objective
Bounded region (polytope)	Always	Impossible
Unbounded region	Possible: $c$ points "into" the region	Possible: $c$ points along a ray

  Unbounded region, unbounded problem:      Unbounded region, finite optimum:
  max x1 + x2,  -x1 + x2 <= 1, x >= 0        min x1 + x2, same region

  x2   /####                                 x2   /####
     /#######  ->  c = (1,1): push            /#######
   /##########       up-right forever       /##########
  *############ ...                         *############   <- optimum (0,0):
  +-------------> x1                       +-------------> x1  c points back into
                                                                the region

In the simplex method:

• Unbounded: an entering variable has no positive entry in its column, so nothing limits how far it can increase (Section 2.3, ratio test).
• Infeasible: Phase I of the two-phase method ends with a positive sum of artificial variables.

from scipy.optimize import linprog

cases = {
    # max x1 + x2  s.t.  -x1 + x2 <= 1, x >= 0           -> unbounded
    "unbounded":  dict(c=[-1, -1], A_ub=[[-1, 1]], b_ub=[1]),
    # min x1 + x2 on the same region                      -> finite optimum at (0,0)
    "finite":     dict(c=[1, 1],   A_ub=[[-1, 1]], b_ub=[1]),
    # x1 + x2 <= 2  and  x1 + x2 >= 3                     -> infeasible
    "infeasible": dict(c=[1, 1],   A_ub=[[1, 1], [-1, -1]], b_ub=[2, -3]),
}
for name, kw in cases.items():
    res = linprog(**kw)                     # default bounds: x >= 0
    print(f"{name:<10} status={res.status}  x={res.x}")
# unbounded  status=3  x=None
# finite     status=0  x=[0. 0.]
# infeasible status=2  x=None

linprog status	Meaning
0	Optimal solution found
1	Iteration limit reached
2	Problem appears infeasible
3	Problem appears unbounded
4	Numerical difficulties

Always check res.status before reading res.x. With status 2 or 3, res.x is not a meaningful solution.

Key takeaways

• Infeasible means an empty region. Unbounded means the objective improves forever along a feasible ray.
• An unbounded region can still have a finite optimum. It depends on the direction of $c$.
• Both outcomes have short certificates: a contradicting constraint combination (infeasible) or an improving ray (unbounded).
• In real models these statuses usually signal modeling errors: missing constraints, wrong signs, or bad data.

Practice

1. (Beginner) For each LP, decide graphically whether it is infeasible, unbounded or optimal: (a) $\max x_1$ s.t. $x_1 - x_2 \le 1$, $x \ge 0$; (b) $\max x_2$ s.t. $x_1 - x_2 \le 1$, $x \ge 0$; (c) $\min x_1$ s.t. $x_1 + x_2 \ge 5$, $x_1 + x_2 \le 3$, $x \ge 0$.
2. (Intermediate) For the unbounded case in the code, find a feasible point $x_0$ and a direction $d$ with $Ad \le 0$, $d \ge 0$ and $c^\top d > 0$, and verify both conditions numerically.
3. (Interview) "A production model returns 'unbounded' with infinite profit. What do you check first?" Hint: which resource or market limit is missing? Look at which variables grow along the reported ray.

Degeneracy

Why it matters. Degeneracy is common in practice: assignment problems, network flows and models with many redundant constraints are full of it. It can make the simplex method perform pivots that do not move at all (stalling) and, in rare cases, loop forever (cycling). It also makes sensitivity reports harder to interpret, because shadow prices stop being unique.

Theory.

• Geometrically, a vertex in $\mathbb{R}^n$ is degenerate when more than $n$ constraints are active there (the vertex is "over-determined").
• Algebraically, a BFS is degenerate when at least one basic variable equals 0.
• The two are linked: if more than $n$ constraints are tight, several different bases describe the same point, and in each of them some basic variable is 0.

  Non-degenerate vertex (n = 2):       Degenerate vertex:
  exactly 2 lines meet                 3 lines meet at one point

        \   /                                \  |  /
         \ /                                  \ | /
          *                                    \|/
         / \                                    *
                                               /|\

Consequences for the simplex method:

Effect	What happens
Degenerate pivot	The basis changes, but the ratio test gives step 0, so the point and objective stay the same
Stalling	Many consecutive degenerate pivots, so there is no progress for a long time
Cycling	A sequence of degenerate pivots returns to an earlier basis, so the method loops forever
Non-unique duals	Different optimal bases give different shadow prices (Section 2.5)

Remedies include anti-cycling pivoting rules (Bland's rule, lexicographic rule, Section 2.3) and the small random perturbations of $b$ that commercial solvers apply.

Example: a degenerate bakery. Add a labor constraint $x_1 + x_2 \le 8$. It passes exactly through the optimum $(2, 6)$, so three constraints (cake oven, mixer, labor) are now tight there:

import numpy as np
from itertools import combinations

# Bakery + labor row x1 + x2 <= 8, standard form with slacks s1..s4
A = np.array([[1, 0, 1, 0, 0, 0],
              [0, 2, 0, 1, 0, 0],
              [3, 2, 0, 0, 1, 0],
              [1, 1, 0, 0, 0, 1]], dtype=float)
b = np.array([4, 12, 18, 8], dtype=float)
names = ["x1", "x2", "s1", "s2", "s3", "s4"]

for basis in combinations(range(6), 4):
    B = A[:, basis]
    if abs(np.linalg.det(B)) < 1e-12:
        continue
    x = np.zeros(6)
    x[list(basis)] = np.linalg.solve(B, b)
    if np.all(x >= -1e-12) and np.allclose(x[:2], [2, 6]):   # BFSs at the vertex (2,6)
        zeros = [names[j] for j in basis if abs(x[j]) < 1e-12]
        print("basis", [names[j] for j in basis], "-> basic variables at 0:", zeros)
# basis ['x1', 'x2', 's1', 's2'] -> basic variables at 0: ['s2']
# basis ['x1', 'x2', 's1', 's3'] -> basic variables at 0: ['s3']
# basis ['x1', 'x2', 's1', 's4'] -> basic variables at 0: ['s4']

Three different bases describe the same point, and each has a basic slack at zero. A simplex pivot between them changes the basis without moving, which is a degenerate pivot.

Key takeaways

• Degenerate vertex: more than $n$ active constraints. Degenerate BFS: a basic variable equals 0.
• Several bases can represent one degenerate vertex, and pivoting among them makes no progress.
• Degeneracy can cause stalling and, rarely, cycling. Anti-cycling rules and perturbation prevent it.
• Under degeneracy, dual values (shadow prices) may not be unique.

Practice

1. (Beginner) Is the origin of the bakery a degenerate vertex? Count the active constraints.
2. (Intermediate) Change the labor constraint to $x_1 + x_2 \le 8.5$ and rerun the code. How many bases give the new optimal vertex now, and why?
3. (Interview) "Why are assignment problems highly degenerate when solved as LPs?" Hint: an $n \times n$ assignment LP has $2n - 1$ independent equality constraints, but a feasible assignment has only $n$ positive variables.

<a id="23-the-simplex-method"></a>

2.3 The Simplex Method

This subchapter develops the simplex method (Dantzig, 1947): start at a vertex, move along an improving edge to an adjacent vertex, and stop when no edge improves. It covers the mechanics (tableau, pivots), the decisions inside each iteration (entering and leaving variables), how to find a starting vertex, and how to avoid infinite loops.

  The simplex loop

  +---------------------------+
  | start at a BFS            |  (all-slack basis, or Phase I / Big-M)
  +-------------+-------------+
                v
  +---------------------------+    no improving variable
  | pricing: any improving    | -------------------------> OPTIMAL
  | nonbasic variable?        |
  +-------------+-------------+
                v yes: choose entering variable (pivoting rule)
  +---------------------------+    no positive entry
  | ratio test: which basic   | -------------------------> UNBOUNDED
  | variable hits 0 first?    |
  +-------------+-------------+
                v
  +---------------------------+
  | pivot: update the basis   | ---> back to pricing
  +---------------------------+

Simplex Tableau

Why it matters. The tableau is a compact table that holds everything the simplex method needs: the current basic solution, its objective value and the reduced costs that tell you whether and how to improve. Working through tableaux by hand is the best way to understand what solvers do, and the final tableau is also the source of all sensitivity information (Section 2.5).

Theory. For a maximization problem in standard form with slacks, write the objective as an equation $z - c^\top x = 0$ and stack it on top of the constraints:

          | x_1 ... x_n (all variables)  | rhs
  --------+------------------------------+------
     z    |  z-row coefficients          | current objective value
  --------+------------------------------+------
   x_B1   |                              | value of x_B1
    ...   |  constraint rows (B^-1 A)    | ...
   x_Bm   |                              | value of x_Bm

Invariants of a valid tableau:

• Each basic variable has a unit column: 1 in its own row, 0 in every other row, including the z-row.
• The rhs column holds the values of the basic variables ($B^{-1}b$) and, in the z-row, the objective value.
• A z-row entry under a nonbasic variable is the negative of its reduced cost, the rate at which $z$ would change if that variable were increased from 0.

For maximization in this convention: a negative z-row entry means increasing that variable improves $z$. When the z-row has no negative entries, the tableau is optimal.

A pivot on element $(r, q)$ (row $r$, column $q$) makes $x_q$ basic in place of the current basic variable of row $r$. It is one step of Gauss-Jordan elimination: divide row $r$ by the pivot element, then subtract multiples of row $r$ from every other row so that column $q$ becomes a unit column.

Example: the bakery, step by step. Exact fractions keep the tableaux readable:

from fractions import Fraction as F

def show(T, basis, names):
    """Print a tableau. Row 0 is the z-row, the last column is the right-hand side."""
    print("basis |" + "".join(f"{n:>7}" for n in names) + " |    rhs")
    labels = ["z"] + [names[j] for j in basis]
    for label, row in zip(labels, T):
        print(f"{label:>5} |" + "".join(f"{str(v):>7}" for v in row[:-1]) + f" | {str(row[-1]):>6}")

def pivot(T, basis, r, q):
    """Gauss-Jordan pivot: variable q enters, the basic variable of row r (r >= 1) leaves."""
    T[r] = [v / T[r][q] for v in T[r]]                 # scale pivot row so the pivot is 1
    for i in range(len(T)):
        if i != r and T[i][q] != 0:                    # clear column q in every other row
            f = T[i][q]
            T[i] = [a - f * b for a, b in zip(T[i], T[r])]
    basis[r - 1] = q                                   # row r now belongs to x_q

names = ["x1", "x2", "s1", "s2", "s3"]
T = [[F(v) for v in row] for row in [
    [-3, -5, 0, 0, 0,  0],     # z - 3x1 - 5x2 = 0
    [ 1,  0, 1, 0, 0,  4],     # bread oven
    [ 0,  2, 0, 1, 0, 12],     # cake oven
    [ 3,  2, 0, 0, 1, 18]]]    # mixer
basis = [2, 3, 4]              # start: all slacks basic -> vertex (0, 0)

show(T, basis, names)
pivot(T, basis, r=2, q=1)      # x2 enters (most negative: -5), s2 leaves (ratio 12/2 = 6)
show(T, basis, names)
pivot(T, basis, r=3, q=0)      # x1 enters (-3), s3 leaves (ratio 6/3 = 2)
show(T, basis, names)

basis |     x1     x2     s1     s2     s3 |    rhs        vertex (0, 0), z = 0
    z |     -3     -5      0      0      0 |      0
   s1 |      1      0      1      0      0 |      4
   s2 |      0      2      0      1      0 |     12
   s3 |      3      2      0      0      1 |     18

basis |     x1     x2     s1     s2     s3 |    rhs        vertex (0, 6), z = 30
    z |     -3      0      0    5/2      0 |     30
   s1 |      1      0      1      0      0 |      4
   x2 |      0      1      0    1/2      0 |      6
   s3 |      3      0      0     -1      1 |      6

basis |     x1     x2     s1     s2     s3 |    rhs        vertex (2, 6), z = 36
    z |      0      0      0    3/2      1 |     36        <- no negative entry: OPTIMAL
   s1 |      0      0      1    1/3   -1/3 |      2
   x2 |      0      1      0    1/2      0 |      6
   x1 |      1      0      0   -1/3    1/3 |      2

Reading the final tableau: $x_1 = 2$, $x_2 = 6$, $s_1 = 2$ (2 idle bread-oven hours), $z = 36$. The z-row entries under the slacks, $3/2$ and $1$, are the shadow prices of the cake oven and the mixer (Section 2.4).

Key takeaways

• The tableau stores the z-row and the constraint rows in the current basis. Basic columns are unit vectors.
• The rhs column gives the basic variable values and the objective value. Nonbasic variables are 0.
• For maximization, a negative z-row entry identifies an improving variable. No negative entries means optimal.
• A pivot is one Gauss-Jordan elimination step that swaps one variable into the basis and one out.

Practice

1. (Beginner) Solve $\max 2x_1 + 3x_2$ s.t. $x_1 + x_2 \le 4$, $x_1 + 3x_2 \le 6$, $x \ge 0$ with tableaux by hand, then check with the code.
2. (Intermediate) In the bakery, start by pivoting $x_1$ in first (instead of $x_2$). How many pivots are needed, and which vertices are visited?
3. (Interview) "What exactly does the number in the z-row under a slack variable mean in the final tableau?" Hint: increase that slack by one unit, which is the same as removing one unit of the resource, and see how $z$ changes.

Pivoting Rules

Why it matters. At each iteration several nonbasic variables may be improving. The rule used to choose one (the pivoting rule) does not affect correctness, but it decides how many iterations are needed, how much each iteration costs, and whether the method can loop forever. Modern solvers spend much of their engineering effort on this choice.

Theory. Geometrically, every improving nonbasic variable corresponds to an improving edge leaving the current vertex. A pivoting rule chooses which edge to walk along.

Rule	Entering variable	Strengths	Weaknesses
Dantzig (largest coefficient)	Most negative z-row entry	Cheap, simple	Depends on scaling. Exponential on Klee-Minty cubes. Can cycle
Greatest improvement	Largest total gain $\lvert \bar c_j \rvert \cdot \theta_j$	Big steps	Needs a ratio test per candidate. Still exponential in the worst case
Steepest edge	Largest $\lvert \bar c_j \rvert / \lVert \eta_j \rVert$ (gain per unit of distance)	Fewest iterations in practice	Needs edge norms $\lVert \eta_j \rVert$, updated each iteration
Devex	Approximate steepest edge	Most of the benefit at a fraction of the cost	Approximation
Bland	Lowest-index improving variable	Provably never cycles	Very slow in practice
Partial pricing	Best among a subset of columns	Cheap for $n \gg m$	More iterations

The leaving variable is fixed by the ratio test (a later topic). Only ties among leaving candidates leave room for a rule.

Example: the same LP, three rules. This builds on show and pivot from the previous topic:

from fractions import Fraction as F

# Reuses pivot() from the Simplex Tableau topic above.

def step(T, j):
    """Length of the step along column j (minimum ratio), or infinity."""
    ratios = [T[i][-1] / T[i][j] for i in range(1, len(T)) if T[i][j] > 0]
    return min(ratios) if ratios else float("inf")

def simplex(T, basis, rule="dantzig", max_iter=100, trace=None):
    """Maximize using a z-row tableau. Returns (status, number of pivots)."""
    for it in range(max_iter):
        z = T[0][:-1]
        improving = [j for j, v in enumerate(z) if v < 0]
        if not improving:
            return "optimal", it
        if rule == "bland":
            q = min(improving)                                # lowest index
        elif rule == "dantzig":
            q = min(improving, key=lambda j: z[j])            # most negative entry
        else:                                                 # "greatest"
            q = max(improving, key=lambda j: -z[j] * step(T, j))
        rows = [i for i in range(1, len(T)) if T[i][q] > 0]
        if not rows:
            return "unbounded", it
        # Ratio test. Ties go to the smallest basic-variable index.
        r = min(rows, key=lambda i: (T[i][-1] / T[i][q], basis[i - 1]))
        pivot(T, basis, r, q)
        if trace:
            trace(T, basis)
    return "iteration limit", max_iter

def bakery():
    rows = [[-3, -5, 0, 0, 0, 0], [1, 0, 1, 0, 0, 4], [0, 2, 0, 1, 0, 12], [3, 2, 0, 0, 1, 18]]
    return [[F(v) for v in r] for r in rows], [2, 3, 4]

def vertex(T, basis):
    vals = {j: T[i + 1][-1] for i, j in enumerate(basis)}
    return f"({vals.get(0, 0)}, {vals.get(1, 0)}) z={T[0][-1]}"

for rule in ["dantzig", "greatest", "bland"]:
    T, basis = bakery()
    path = ["(0, 0) z=0"]
    status, its = simplex(T, basis, rule, trace=lambda T, b: path.append(vertex(T, b)))
    print(f"{rule:<9} {its} pivots: " + " -> ".join(path))
# dantzig   2 pivots: (0, 0) z=0 -> (0, 6) z=30 -> (2, 6) z=36
# greatest  2 pivots: (0, 0) z=0 -> (0, 6) z=30 -> (2, 6) z=36
# bland     3 pivots: (0, 0) z=0 -> (4, 0) z=12 -> (4, 3) z=27 -> (2, 6) z=36

  Two routes around the bakery polygon:

  x2
  6 *<---------* (2,6)      Dantzig / greatest:  (0,0) -> (0,6) -> (2,6)
    ^ Dantzig   ^           Bland:               (0,0) -> (4,0) -> (4,3) -> (2,6)
    |            \
  3 |            * (4,3)
    |            ^ Bland
    *----------->*-----> x1
  (0,0)        (4,0)

Both routes reach the same optimum, because pivoting rules change the path, not the answer. Bland's rule picks $x_1$ only because it has the lower index, and goes the long way around.

Key takeaways

• A pivoting rule chooses the entering variable, which means choosing an improving edge out of the current vertex.
• Dantzig's rule is simple but scale-dependent. Steepest edge (and its approximation Devex) is the practical standard in modern solvers.
• Bland's rule guarantees termination but is slow, so it is used as a fallback against cycling.
• For every known deterministic rule, there are LPs that force exponentially many pivots (Section 2.6).

Practice

1. (Beginner) Multiply the bakery's cake-oven row by 10 ($20x_2 \le 120$). Does Dantzig's rule choose a different first entering variable? Does the solution change?
2. (Intermediate) Add a "steepest" rule to simplex that divides $\lvert \bar c_j \rvert$ by the Euclidean norm of the edge direction (the column $-B^{-1}A_j$ for basic variables, plus a 1 for $x_j$ itself). Compare pivot counts on a few random LPs.
3. (Interview) "Why don't solvers simply use the rule that gives the biggest improvement per iteration?" Hint: compare the cost of one iteration (all ratio tests) with the number of iterations saved.

Optimality and Entering Variable Selection

Why it matters. The simplex method needs a test that proves a vertex is optimal without comparing it to all other vertices. Reduced costs provide exactly that: a local check of the edges leaving the current vertex is enough to certify global optimality, thanks to convexity.

Theory. Split the variables into basic $x_B$ and nonbasic $x_N$. From $Bx_B + Nx_N = b$, we get $x_B = B^{-1}b - B^{-1}Nx_N$. Substitute into the objective $c^\top x = c_B^\top x_B + c_N^\top x_N$:

$$ z = \underbrace{c_B^\top B^{-1} b}{\text{current value}} + \sum{j \in N} \underbrace{\left(c_j - c_B^\top B^{-1} A_j\right)}_{\bar c_j \text{ (reduced cost)}} x_j . $$

• $y^\top = c_B^\top B^{-1}$ are the simplex multipliers, which are also the dual variables (Section 2.4).
• $\bar c_j = c_j - y^\top A_j$ is the change in $z$ per unit increase of nonbasic $x_j$, when the basic variables adjust to keep $Ax = b$.
• For basic variables, $\bar c_j = 0$ always.

Optimality test (maximization): if $\bar c_j \le 0$ for all nonbasic $j$, the current BFS is optimal. The reason is that any feasible $x$ has $x_N \ge 0$, so $z \le c_B^\top B^{-1}b$. For minimization, flip the sign: optimal when all $\bar c_j \ge 0$.

Alternative optima: if the BFS is optimal and some nonbasic $\bar c_j = 0$, pivoting $x_j$ in reaches another optimal solution. The LP has multiple optima (the iso-profit line is parallel to an edge).

Entering variable selection: any $j$ with $\bar c_j > 0$ (maximization) improves the objective. The pivoting rule picks one.

import numpy as np

A = np.array([[1, 0, 1, 0, 0],
              [0, 2, 0, 1, 0],
              [3, 2, 0, 0, 1]], dtype=float)
names = ["x1", "x2", "s1", "s2", "s3"]

def reduced_costs(c, basis):
    B = A[:, basis]
    y = np.linalg.solve(B.T, c[basis])            # y^T = c_B^T B^-1  (simplex multipliers)
    cbar = c - y @ A                              # c_j - y^T A_j for every column
    return {names[j]: float(round(cbar[j], 3)) for j in range(len(c)) if j not in basis}

c = np.array([3, 5, 0, 0, 0], dtype=float)
print("start   {s1,s2,s3}:", reduced_costs(c, [2, 3, 4]))   # {'x1': 3.0, 'x2': 5.0}  -> improve
print("middle  {s1,x2,s3}:", reduced_costs(c, [2, 1, 4]))   # {'x1': 3.0, 's2': -2.5} -> x1 improves
print("optimal {s1,x2,x1}:", reduced_costs(c, [2, 1, 0]))   # {'s2': -1.5, 's3': -1.0} -> optimal

c_alt = np.array([7.5, 5, 0, 0, 0])                          # bread profit raised to 7.5
print("c1=7.5  {s1,x2,x1}:", reduced_costs(c_alt, [2, 1, 0]))  # {'s2': 0.0, 's3': -2.5}
# s2 has reduced cost 0 -> pivoting s2 in gives another optimal vertex, (4, 3)

Key takeaways

• Reduced cost $\bar c_j = c_j - c_B^\top B^{-1}A_j$: the rate of change of $z$ when nonbasic $x_j$ increases.
• Maximization is optimal when all reduced costs are $\le 0$. Minimization is optimal when all are $\ge 0$.
• A zero reduced cost on a nonbasic variable at the optimum signals alternative optimal solutions.
• The check is local (edges at one vertex) but certifies a global optimum because LPs are convex.

Practice

1. (Beginner) Using the code, compute the reduced costs at the basis ${x_1, s_2, s_3}$ (vertex $(4, 0)$). Which variables could enter?
2. (Intermediate) Prove that $\bar c_j = 0$ for every basic variable, directly from the formula.
3. (Interview) "The simplex method only looks at neighboring vertices. Why can it claim global optimality?" Hint: the formula for $z$ above holds for every feasible point, not only for neighbors.

Ratio Test and Leaving Variable

Why it matters. Choosing the entering variable says which edge to walk along. The ratio test says how far to walk: as far as possible while staying feasible. It also detects unboundedness (the edge never ends) and degeneracy (the step has length zero).

Theory. Let $x_q$ enter. Keeping the other nonbasic variables at 0 and increasing $x_q = \theta \ge 0$, the basic variables change as

$$ x_B(\theta) = \bar b - \theta , \bar a_q, \qquad \bar b = B^{-1} b, \quad \bar a_q = B^{-1} A_q . $$

• A row with $\bar a_{iq} \le 0$ never limits $\theta$, because the basic variable stays the same or grows.
• A row with $\bar a_{iq} > 0$ reaches zero at $\theta = \bar b_i / \bar a_{iq}$.

The minimum ratio test picks the first basic variable to hit zero:

$$ \theta^* = \min_{i ,:, \bar a_{iq} > 0} \frac{\bar b_i}{\bar a_{iq}}, \qquad \text{leaving variable} = x_{B_r} \text{ where } r \text{ attains the minimum.} $$

Situation	Meaning
No $\bar a_{iq} > 0$	The edge is an infinite ray: unbounded
$\theta^* > 0$	Move to a new vertex. The objective improves by $\lvert \bar c_q \rvert \theta^*$
$\theta^* = 0$	Degenerate pivot: basis changes, point does not move
Ties in the minimum	Several variables hit 0 together, so the next BFS is degenerate

Analogy: you are filling several tanks from one pipe. Each tank (row) fills at its own rate. The first tank to overflow decides when you must stop, and tanks that drain or stay level never stop you.

  Bakery, first pivot: x2 enters (moving up the x2 axis from (0,0))

  x2
  9 +  <- mixer would stop us at x2 = 18/2 = 9
    |
  6 +  <- cake oven stops us first at x2 = 12/2 = 6   => theta* = 6, s2 leaves
    |
    |     bread oven never stops us (coefficient 0 in column x2)
    *-------------> x1

from fractions import Fraction as F

def ratio_test(column, rhs, basic_names):
    """Return the index of the leaving row, or None if the entering column is unbounded."""
    best = None
    for i, (a, b) in enumerate(zip(column, rhs)):
        if a > 0:
            ratio = F(b) / F(a)
            print(f"  {basic_names[i]}: {b}/{a} = {ratio}")
            if best is None or ratio < best[0]:
                best = (ratio, i)
        else:
            print(f"  {basic_names[i]}: coefficient {a} <= 0, no limit")
    if best is None:
        print("  -> no positive entry: UNBOUNDED")
        return None
    print(f"  -> theta* = {best[0]}, leaving variable {basic_names[best[1]]}")
    return best[1]

print("Bakery, x2 entering at (0,0):")
ratio_test(column=[0, 2, 2], rhs=[4, 12, 18], basic_names=["s1", "s2", "s3"])
print("Bakery, x1 entering at (0,6):")
ratio_test(column=[1, 0, 3], rhs=[4, 6, 6], basic_names=["s1", "x2", "s3"])
print("max x1 + x2 s.t. -x1 + x2 <= 1, x1 entering:")
ratio_test(column=[-1], rhs=[1], basic_names=["s1"])
# Bakery, x2 entering at (0,0):
#   s1: coefficient 0 <= 0, no limit
#   s2: 12/2 = 6
#   s3: 18/2 = 9
#   -> theta* = 6, leaving variable s2
# Bakery, x1 entering at (0,6):
#   s1: 4/1 = 4
#   x2: coefficient 0 <= 0, no limit
#   s3: 6/3 = 2
#   -> theta* = 2, leaving variable s3
# max x1 + x2 s.t. -x1 + x2 <= 1, x1 entering:
#   s1: coefficient -1 <= 0, no limit
#   -> no positive entry: UNBOUNDED

Key takeaways

• The ratio test considers only rows with a positive entry in the entering column and picks the smallest ratio $\bar b_i / \bar a_{iq}$.
• It keeps the next solution feasible by stopping exactly when the first basic variable reaches 0.
• No positive entry means unbounded. A zero ratio means a degenerate pivot. Ties create degenerate vertices.
• Including rows with zero or negative entries, or picking a non-minimal ratio, produces an infeasible basis.

Practice

1. (Beginner) Explain why a row with $\bar a_{iq} < 0$ is skipped, using the formula $x_{B_i}(\theta) = \bar b_i - \theta \bar a_{iq}$.
2. (Intermediate) Construct a tableau row set where the ratio test has a tie, pivot on one of the tied rows, and show that the next BFS is degenerate.
3. (Interview) "Solvers use tolerances in the ratio test, for example ignoring pivot elements smaller than $10^{-9}$. Why?" Hint: think about dividing by a tiny, possibly round-off-polluted number. Look up the Harris ratio test.

Two-Phase Simplex Method

Why it matters. The simplex method must start at a BFS. With only $\le$ constraints and $b \ge 0$, the slacks provide one for free. Many real models contain $\ge$ and $=$ constraints (demands, balances, blends), where no starting vertex is obvious, and the model may even be infeasible. The two-phase method solves both problems in a clean, numerically safe way.

Theory.

1. Write the LP in standard form with $b \ge 0$. For every row that has no slack variable with coefficient $+1$ (all $\ge$ and $=$ rows), add an artificial variable $a_i \ge 0$. The artificial columns plus the slacks now form an identity matrix, so there is an obvious starting basis.
2. Phase I: minimize $w = \sum_i a_i$ subject to the augmented constraints.
   • If $w^* > 0$, no solution has all artificials at zero, so the original LP is infeasible.
   • If $w^* = 0$, the final basis is a BFS of the original LP. Pivot out any artificial that is still basic at value 0 (a degenerate case), then delete the artificial columns.
3. Phase II: restore the original objective, price out the basic columns (make their z-row entries zero by subtracting multiples of their rows), and run the ordinary simplex method.

Artificial variables have no meaning in the model. They are scaffolding: you build it to get started and remove it before Phase II.

Example.

$$ \min\ 4x_1 + x_2 \quad \text{s.t.} \quad 3x_1 + x_2 = 3, \quad 4x_1 + 3x_2 \ge 6, \quad x_1 + 2x_2 \le 4, \quad x \ge 0 . $$

Standard form: $3x_1 + x_2 + a_1 = 3$, $4x_1 + 3x_2 - e_2 + a_2 = 6$, $x_1 + 2x_2 + s_3 = 4$. The starting basis is ${a_1, a_2, s_3}$.

from fractions import Fraction as F

# Reuses show(), pivot() from the Simplex Tableau topic and simplex() from Pivoting Rules.

names = ["x1", "x2", "e2", "s3", "a1", "a2"]
rows = [[3, 1,  0, 0, 1, 0, 3],     # 3x1 +  x2           + a1      = 3
        [4, 3, -1, 0, 0, 1, 6],     # 4x1 + 3x2 - e2           + a2 = 6
        [1, 2,  0, 1, 0, 0, 4]]     #  x1 + 2x2      + s3           = 4
T = [[F(0)] * 7] + [[F(v) for v in r] for r in rows]
basis = [4, 5, 3]

# Phase I: min a1 + a2  ==  max -(a1 + a2)  ->  z-row: z + a1 + a2 = 0
T[0][4] = T[0][5] = F(1)
for i in (1, 2):                                   # price out basic a1, a2
    T[0] = [z - v for z, v in zip(T[0], T[i])]
simplex(T, basis)
print("Phase I optimum w* =", -T[0][-1])           # 0 -> the original LP is feasible
show(T, basis, names)

# Phase II: drop artificial columns, restore min 4x1 + x2  ==  max -4x1 - x2
T = [row[:4] + [row[-1]] for row in T]
T[0] = [F(4), F(1), F(0), F(0), F(0)]              # z + 4x1 + x2 = 0
for i, j in enumerate(basis, start=1):             # price out the current basic columns
    if T[0][j] != 0:
        f = T[0][j]
        T[0] = [z - f * v for z, v in zip(T[0], T[i])]
simplex(T, basis)
show(T, basis, names[:4])
print("min cost =", -T[0][-1])                     # 17/5

Phase I optimum w* = 0
basis |     x1     x2     e2     s3     a1     a2 |    rhs
    z |      0      0      0      0      1      1 |      0
   x1 |      1      0    1/5      0    3/5   -1/5 |    3/5
   x2 |      0      1   -3/5      0   -4/5    3/5 |    6/5
   s3 |      0      0      1      1      1     -1 |      1
basis |     x1     x2     e2     s3 |    rhs
    z |      0      0      0    1/5 |  -17/5
   x1 |      1      0      0   -1/5 |    2/5
   x2 |      0      1      0    3/5 |    9/5
   e2 |      0      0      1      1 |      1
min cost = 17/5

Phase I found the feasible vertex $(3/5, 6/5)$. Phase II improved it to the optimum $(2/5, 9/5)$ with cost $17/5$.

Key takeaways

• Artificial variables create an identity starting basis for $\ge$ and $=$ rows. They must be driven to zero.
• Phase I minimizes the sum of artificials. A positive optimum proves the LP is infeasible.
• Phase II starts from Phase I's final basis with the true objective, after pricing out basic columns.
• Basic artificials at value 0 after Phase I must be pivoted out (or their row is redundant) before Phase II.

Practice

1. (Beginner) Which rows of the example need an artificial variable, and why does row 3 not need one?
2. (Intermediate) Change the first constraint to $3x_1 + x_2 = 30$ and the third to $x_1 + 2x_2 \le 4$. Run Phase I and interpret $w^*$.
3. (Interview) "How does a solver know an LP is infeasible, and what can it tell you about why?" Hint: Phase I's final dual values give a Farkas certificate (Section 2.4), a weighted combination of constraints that is contradictory.

Big-M Method

Why it matters. The Big-M method is the single-phase alternative to two-phase simplex. It is common in textbooks and spreadsheets, and the same "large penalty" idea appears everywhere in integer programming (Chapter 4). Understanding why it is numerically risky explains why solvers prefer two phases.

Theory. Add artificial variables as in two-phase, but solve one LP with a huge penalty $M$ on them:

$$ \min\ c^\top x + M \sum_i a_i \qquad (\text{or } \max\ c^\top x - M \sum_i a_i). $$

If $M$ is large enough, any solution with an artificial variable at a positive value is so expensive that the simplex method drives all artificials to 0 if it can. If an artificial remains positive at the optimum, the original LP is infeasible (provided $M$ was large enough).

	Two-phase	Big-M
Number of LPs	Two	One
Extra parameter	None	$M$, which must be "large enough"
Numerics	Safe	Mixing $M \approx 10^6$ with data $\approx 1$ causes round-off
Infeasibility detection	Clean: $w^* > 0$	Ambiguous if $M$ is too small
Used in	Commercial solvers	Hand calculations, teaching, quick spreadsheet models

By hand, $M$ is treated symbolically: z-row entries look like $4 - 7M$, and you compare them by the coefficient of $M$ first.

Example. The same LP as in the two-phase topic, with numeric values of $M$:

import numpy as np
from scipy.optimize import linprog

# Columns: x1, x2, e2, s3, a1, a2
A_eq = np.array([[3, 1,  0, 0, 1, 0],
                 [4, 3, -1, 0, 0, 1],
                 [1, 2,  0, 1, 0, 0]], dtype=float)
b_eq = np.array([3, 6, 4], dtype=float)

for M in [1e6, 10, 1]:
    c = [4, 1, 0, 0, M, M]                       # artificials penalized by M
    res = linprog(c, A_eq=A_eq, b_eq=b_eq)
    x1, x2, *_, a1, a2 = res.x
    print(f"M = {M:>9g}: x = ({x1:.2f}, {x2:.2f}), a = ({a1:g}, {a2:g}), objective {res.fun:.2f}")
# M =     1e+06: x = (0.40, 1.80), a = (0, 0), objective 3.40   <- correct
# M =        10: x = (0.40, 1.80), a = (0, 0), objective 3.40   <- still large enough
# M =         1: x = (0.00, 2.00), a = (1, 0), objective 3.00   <- WRONG: a1 = 1 > 0

With $M = 1$, violating the first constraint by one unit costs only 1, which is cheaper than satisfying it. The "optimum" is infeasible for the original problem, and nothing in the output warns you except the positive artificial.

Key takeaways

• Big-M adds artificials with a huge cost $M$ and solves a single LP.
• If an artificial is positive at the optimum (with $M$ large enough), the original LP is infeasible.
• Too small an $M$ gives wrong answers. Too large an $M$ causes numerical errors.
• Prefer two-phase in code. Use Big-M for hand calculations with symbolic $M$.

Practice

1. (Beginner) Write the initial Big-M tableau (with symbolic $M$) for the example and identify the first entering variable.
2. (Intermediate) Find, by bisection with the code above, the smallest $M$ for which the method returns the correct solution. Relate it to the dual value of the first constraint.
3. (Interview) "Your colleague's model uses $M = 10^{12}$ and the solver reports slightly infeasible solutions. What is going on?" Hint: floating-point numbers have about 16 significant digits, and the solver's feasibility tolerance is around $10^{-6}$ relative to the coefficients.

Cycling and Bland's Rule

Why it matters. Each nondegenerate pivot strictly improves the objective, so the simplex method can never revisit a basis, and it terminates. Degenerate pivots do not improve anything, so a sequence of them can return to an earlier basis and repeat forever. Bland's rule is the simplest guarantee that this never happens.

Theory.

• Cycling: a sequence of degenerate pivots $B_1 \to B_2 \to \dots \to B_k = B_1$. The objective stays constant and the algorithm never stops.
• Cycling needs degeneracy, and in practice it is rare. Stalling (very long but finite degenerate sequences) is much more common in large degenerate models.

Bland's rule (smallest-subscript rule, 1977):

1. Entering: among all improving variables, choose the one with the smallest index.
2. Leaving: among rows tied in the ratio test, choose the basic variable with the smallest index.

Theorem. With Bland's rule, the simplex method never cycles and therefore terminates in finitely many iterations.

Anti-cycling technique	Idea	Used in practice?
Bland's rule	Smallest-index entering and leaving	As a fallback (slow in general)
Lexicographic rule	Break ratio-test ties lexicographically on rows of $B^{-1}$	In theory and some implementations
Perturbation	Add tiny random amounts to $b$ to destroy degeneracy	Yes, the standard in commercial solvers

Example: Beale's cycling LP.

$$ \begin{aligned} \max\ & \tfrac34 x_4 - 20 x_5 + \tfrac12 x_6 - 6 x_7 \ \text{s.t.}\ & x_1 + \tfrac14 x_4 - 8x_5 - x_6 + 9x_7 = 0 \ & x_2 + \tfrac12 x_4 - 12x_5 - \tfrac12 x_6 + 3x_7 = 0 \ & x_3 + x_6 = 1, \qquad x \ge 0 \end{aligned} $$

from fractions import Fraction as F

# Reuses pivot() from the Simplex Tableau topic and simplex() from Pivoting Rules.

def beale():
    rows = [[0, 0, 0, F(-3, 4), 20, F(-1, 2), 6, 0],      # z-row
            [1, 0, 0, F(1, 4), -8, -1, 9, 0],
            [0, 1, 0, F(1, 2), -12, F(-1, 2), 3, 0],
            [0, 0, 1, 0, 0, 1, 0, 1]]
    return [[F(v) for v in r] for r in rows], [0, 1, 2]    # basis {x1, x2, x3}

def tracer(start):
    seen = [sorted(j + 1 for j in start)]                  # the starting basis counts as seen
    def trace(T, basis):
        key = sorted(j + 1 for j in basis)
        tag = "  <- seen before: CYCLING" if key in seen else ""
        seen.append(key)
        print(f"  basis {key}  z = {T[0][-1]}{tag}")
    return trace

print("Dantzig's rule:")
T, basis = beale()
print(simplex(T, basis, "dantzig", max_iter=6, trace=tracer(basis)))

print("Bland's rule:")
T, basis = beale()
print(simplex(T, basis, "bland", trace=tracer(basis)))

Dantzig's rule:
  basis [2, 3, 4]  z = 0
  basis [3, 4, 5]  z = 0
  basis [3, 5, 6]  z = 0
  basis [3, 6, 7]  z = 0
  basis [1, 3, 7]  z = 0
  basis [1, 2, 3]  z = 0  <- seen before: CYCLING
('iteration limit', 6)
Bland's rule:
  basis [2, 3, 4]  z = 0
  basis [3, 4, 5]  z = 0
  basis [3, 5, 6]  z = 0
  basis [1, 3, 6]  z = 0
  basis [1, 2, 6]  z = 1/2
  basis [1, 4, 6]  z = 5/4
('optimal', 6)

Dantzig's rule returns to the starting basis after 6 degenerate pivots and would repeat forever. Bland's rule deviates at the fourth pivot (it picks $x_1$ instead of $x_7$), escapes the degenerate vertex and reaches the optimum $z = 5/4$.

Key takeaways

• Cycling is an infinite loop of degenerate pivots. It requires degeneracy and is rare, but possible.
• Bland's rule (smallest index for entering and leaving) provably prevents cycling.
• Bland's rule is slow in general, so solvers use perturbation or switch to Bland's rule only when stalling is detected.
• Stalling (long finite degenerate sequences) is the practical problem, and it is also addressed by perturbation.

Practice

1. (Beginner) In Beale's LP, why is every pivot of the Dantzig run degenerate? Look at the right-hand sides of the first two rows.
2. (Intermediate) Implement the perturbation method: add $\varepsilon^i$ (e.g. $\varepsilon = 10^{-3}$) to the right-hand side of row $i$ and run Dantzig's rule on Beale's LP. Does it still cycle?
3. (Interview) "If cycling is so rare, why do textbooks spend so much time on it?" Hint: termination proofs. Without an anti-cycling rule, the simplex method is not even guaranteed to be an algorithm.

Revised Simplex Method

Why it matters. The tableau updates all $m \times n$ entries at every pivot, even though each iteration only needs the z-row, one column and the rhs. Real models have millions of mostly-zero columns. The revised simplex method computes only what it needs from the original sparse data and a factorization of $B$, and it is the form implemented in every serious LP solver.

Theory. Each iteration needs only:

Step	Computation	Tableau equivalent
1. Basic solution	Solve $B x_B = b$	rhs column
2. Pricing	Solve $B^\top y = c_B$, then $\bar c_j = c_j - y^\top A_j$	z-row
3. Entering column	Solve $B u = A_q$	Column $q$
4. Ratio test	$\min {x_{B_i} / u_i : u_i > 0}$	Same
5. Update	Replace column $r$ of $B$ by $A_q$	Full pivot

The algorithm never forms $B^{-1}$ explicitly. It keeps an LU factorization of $B$ and updates it after each basis change (product form of the inverse, Forrest-Tomlin update). It refactorizes from scratch every ~100 iterations for numerical stability.

	Tableau simplex	Revised simplex
Stores	Dense $m \times n$ tableau	Original sparse $A$ + factorization of $B$ ($m \times m$)
Work per iteration	$O(mn)$ always	$O(m^2)$ + pricing, much less with sparsity
Exploits sparsity	No (fill-in destroys it)	Yes
Numerical control	Errors accumulate	Periodic refactorization
Best for	Hand calculation, teaching	All practical implementations

When $n \gg m$, pricing (step 2) dominates. Partial pricing scans only part of the columns each iteration, and column generation (Section 2.6) takes this to its limit.

import numpy as np

def revised_simplex(c, A, b, basis, names, max_iter=100):
    """min c^T x  s.t.  A x = b, x >= 0, from a feasible starting basis."""
    m, n = A.shape
    basis = list(basis)
    for it in range(1, max_iter + 1):
        B = A[:, basis]                                # solvers keep an LU factorization instead
        x_B = np.linalg.solve(B, b)                    # 1. basic solution
        y = np.linalg.solve(B.T, c[basis])             # 2. multipliers: B^T y = c_B
        nonbasic = [j for j in range(n) if j not in basis]
        cbar = {j: c[j] - y @ A[:, j] for j in nonbasic}   # pricing
        q = min(cbar, key=cbar.get)                    # most negative (minimization)
        if cbar[q] >= -1e-12:
            x = np.zeros(n)
            x[basis] = x_B
            return x, y
        u = np.linalg.solve(B, A[:, q])                # 3. entering column B^-1 A_q
        if np.all(u <= 1e-12):
            raise ValueError("LP is unbounded")
        ratios = np.where(u > 1e-12, x_B / np.where(u > 1e-12, u, 1), np.inf)
        r = int(np.argmin(ratios))                     # 4. ratio test
        print(f"iter {it}: {names[q]} enters, {names[basis[r]]} leaves, "
              f"step {ratios[r]:g}, objective {c[basis] @ x_B:g} -> "
              f"{c[basis] @ x_B + cbar[q] * ratios[r]:g}")
        basis[r] = q                                   # 5. basis update
    raise RuntimeError("iteration limit")

A = np.array([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=float)
b = np.array([4, 12, 18], dtype=float)
c = np.array([-3, -5, 0, 0, 0], dtype=float)           # min -profit
x, y = revised_simplex(c, A, b, basis=[2, 3, 4], names=["x1", "x2", "s1", "s2", "s3"])
print("x =", x, " y =", y)
# iter 1: x2 enters, s2 leaves, step 6, objective 0 -> -30
# iter 2: x1 enters, s3 leaves, step 2, objective -30 -> -36
# x = [2. 6. 2. 0. 0.]  y = [ 0.  -1.5 -1. ]

The multipliers $y$ computed for pricing are the dual solution (with signs from the minimization form): $0$, $1.5$ and $1$ per hour of bread oven, cake oven and mixer.

Key takeaways

• Revised simplex performs the same pivots as the tableau method, but computes only the needed pieces from $A$ and a factorization of $B$.
• Each iteration solves a few $m \times m$ systems ($Bx_B = b$, $B^\top y = c_B$, $Bu = A_q$) and prices nonbasic columns.
• It exploits sparsity, controls round-off by refactorizing, and is the basis of all production simplex codes.
• The pricing vector $y = B^{-\top}c_B$ is the dual solution at optimality.

Practice

1. (Beginner) For the final basis ${s_1, x_2, x_1}$ of the bakery, write $B$ and verify $x_B = B^{-1}b$ by hand.
2. (Intermediate) Replace np.linalg.solve with a scipy.linalg.lu_factor computed once per iteration and reused for all three solves (use trans=1 for $B^\top$).
3. (Interview) "An LP has 1000 rows and 10 million columns. Which part of the simplex iteration dominates, and how would you speed it up?" Hint: pricing touches every column. Consider partial pricing, parallel pricing, or not storing the columns at all (column generation).

<a id="24-duality-theory"></a>

2.4 Duality Theory

Every LP has a twin, its dual, built from the same data. This subchapter shows how to construct it, proves that the two have equal optimal values, and explains why dual variables are prices: the value of one more unit of each resource.

The bakery primal (P) and its dual (D) are used throughout:

$$ \text{(P)}\ \max\ 3x_1 + 5x_2 \ \text{ s.t. } \begin{cases} x_1 \le 4 & (y_1) \ 2x_2 \le 12 & (y_2) \ 3x_1 + 2x_2 \le 18 & (y_3) \end{cases} x \ge 0 \qquad \text{(D)}\ \min\ 4y_1 + 12y_2 + 18y_3 \ \text{ s.t. } \begin{cases} y_1 + 3y_3 \ge 3 & (x_1) \ 2y_2 + 2y_3 \ge 5 & (x_2) \end{cases} y \ge 0 $$

Constructing the Dual Problem

Why it matters. The dual is not an arbitrary transformation. It answers a natural question: "what is the best upper bound on the primal optimum that I can prove by combining the constraints?" Being able to write the dual quickly and correctly is needed for optimality certificates, sensitivity analysis, decomposition methods and many proofs in later chapters.

Theory. Deriving the dual from bounds. Multiply each bakery constraint by a nonnegative weight $y_i$ and add them up:

$$ y_1(x_1) + y_2(2x_2) + y_3(3x_1 + 2x_2) \le 4y_1 + 12y_2 + 18y_3 ;\Longleftrightarrow; (y_1 + 3y_3),x_1 + (2y_2 + 2y_3),x_2 \le 4y_1 + 12y_2 + 18y_3 . $$

If the coefficient of each $x_j$ is at least the objective coefficient ($y_1 + 3y_3 \ge 3$ and $2y_2 + 2y_3 \ge 5$), then because $x \ge 0$,

$$ 3x_1 + 5x_2 \le (y_1 + 3y_3)x_1 + (2y_2 + 2y_3)x_2 \le 4y_1 + 12y_2 + 18y_3 . $$

So every such $y$ gives an upper bound on the profit. The dual asks for the best (smallest) such bound, which is exactly problem (D). For example, $y = (0, 1.5, 1)$ gives the bound $18 + 18 = 36$, which equals the profit of $(2, 6)$, so that plan is proven optimal.

General recipe (primal max, dual min):

Primal (max)		Dual (min)
Objective coefficients $c$	$\leftrightarrow$	Right-hand sides $c$
Right-hand sides $b$	$\leftrightarrow$	Objective coefficients $b$
Constraint matrix $A$	$\leftrightarrow$	Transpose $A^\top$
Constraint $i$: $\le$	$\leftrightarrow$	$y_i \ge 0$
Constraint $i$: $\ge$	$\leftrightarrow$	$y_i \le 0$
Constraint $i$: $=$	$\leftrightarrow$	$y_i$ free
Variable $x_j \ge 0$	$\leftrightarrow$	Constraint $j$: $\ge$
Variable $x_j \le 0$	$\leftrightarrow$	Constraint $j$: $\le$
Variable $x_j$ free	$\leftrightarrow$	Constraint $j$: $=$

Read the table right to left to dualize a minimization problem. Mnemonic: "natural" constraints ($\le$ in a max problem, $\ge$ in a min problem) correspond to "natural" variables ($\ge 0$). Anything unnatural maps to unnatural, and equalities map to free variables. The dual of the dual is the primal.

In matrix form, the symmetric pair is:

$$ \max{c^\top x : Ax \le b,\ x \ge 0} \quad \longleftrightarrow \quad \min{b^\top y : A^\top y \ge c,\ y \ge 0}. $$

import numpy as np
from scipy.optimize import linprog

def build_dual(c, A, b, row_senses, var_signs):
    """Primal: max c^T x  with row senses ('<=', '>=', '=') and variable signs
    ('>=0', '<=0', 'free'). Returns the dual: min b^T y  with its own senses and signs."""
    y_sign = {"<=": ">=0", ">=": "<=0", "=": "free"}         # row sense  -> dual variable sign
    dual_sense = {">=0": ">=", "<=0": "<=", "free": "="}     # var sign   -> dual row sense
    return b, A.T, c, [dual_sense[s] for s in var_signs], [y_sign[s] for s in row_senses]

def solve(obj, A, rhs, senses, signs, maximize):
    """Small wrapper around linprog for mixed senses and signs."""
    bound = {">=0": (0, None), "<=0": (None, 0), "free": (None, None)}
    A_ub, b_ub, A_eq, b_eq = [], [], [], []
    for row, r, s in zip(A, rhs, senses):
        if s == "<=":
            A_ub.append(row); b_ub.append(r)
        elif s == ">=":
            A_ub.append(-row); b_ub.append(-r)       # flip to <=
        else:
            A_eq.append(row); b_eq.append(r)
    sgn = -1 if maximize else 1
    res = linprog(sgn * np.asarray(obj), A_ub=A_ub or None, b_ub=b_ub or None,
                  A_eq=A_eq or None, b_eq=b_eq or None, bounds=[bound[s] for s in signs])
    return res.x, sgn * res.fun

# A primal with every kind of row and variable:
#   max 3x1 + 2x2 + 8x3
#   s.t.  x1 + x2 + x3 <= 8,   x1 - x2 >= -2,   x1 + 2x3 = 4
#         x1 >= 0, x2 free, x3 <= 0
c = np.array([3., 2, 8])
A = np.array([[1., 1, 1], [1, -1, 0], [1, 0, 2]])
b = np.array([8., -2, 4])
senses, signs = ["<=", ">=", "="], [">=0", "free", "<=0"]

x, primal_value = solve(c, A, b, senses, signs, maximize=True)
d_obj, d_A, d_rhs, d_senses, d_signs = build_dual(c, A, b, senses, signs)
y, dual_value = solve(d_obj, d_A, d_rhs, d_senses, d_signs, maximize=False)

print("dual rows :", d_senses, " dual vars:", d_signs)
print("primal x =", x, "value", primal_value)
print("dual   y =", y, "value", dual_value)
# dual rows : ['>=', '=', '<=']  dual vars: ['>=0', '<=0', 'free']
# primal x = [4. 4. 0.] value 20.0
# dual   y = [2. 0. 1.] value 20.0      <- equal optimal values (strong duality)

Key takeaways

• The dual looks for the best bound on the primal obtained from nonnegative combinations of its constraints.
• Primal max with $Ax \le b$, $x \ge 0$ corresponds to dual min with $A^\top y \ge c$, $y \ge 0$. Rows become columns and $b$ and $c$ swap roles.
• Use the correspondence table for mixed senses: $=$ maps to a free variable, and unnatural senses map to nonpositive variables.
• One dual variable per primal constraint, one dual constraint per primal variable. The dual of the dual is the primal.

Practice

1. (Beginner) Write the dual of the diet problem from Section 2.1 ($\min c^\top x$ s.t. $Nx \ge r$, $x \ge 0$) and interpret its variables.
2. (Intermediate) Show by hand that dualizing the dual of the symmetric pair gives back the primal. Write the dual as a max problem first.
3. (Interview) "Why does an equality constraint in the primal give a free dual variable?" Hint: in the bounding argument, an equality may be multiplied by a negative weight without breaking the inequality direction.

Weak Duality Theorem

Why it matters. Weak duality is a one-line result with large practical consequences. Any feasible dual solution bounds the primal optimum, so solvers can measure how far a current solution is from optimal (the duality gap) and stop early with a guarantee. Branch-and-bound for integer programs (Chapter 4) relies on exactly this kind of bound.

Theory. For the pair $\max{c^\top x : Ax \le b, x \ge 0}$ and $\min{b^\top y : A^\top y \ge c, y \ge 0}$: if $x$ is primal feasible and $y$ is dual feasible, then

$$ c^\top x ;\le; b^\top y . $$

Proof. $c^\top x \le (A^\top y)^\top x = y^\top A x \le y^\top b$. The first inequality uses $A^\top y \ge c$ and $x \ge 0$. The second uses $Ax \le b$ and $y \ge 0$. $\blacksquare$

Corollaries:

1. Certificate of optimality: if $c^\top x = b^\top y$ for feasible $x, y$, then both are optimal.
2. Primal unbounded $\Rightarrow$ dual infeasible, since no finite bound can exist. Symmetrically, dual unbounded implies primal infeasible.
3. Duality gap: for any feasible pair, $b^\top y - c^\top x \ge 0$ bounds how far $x$ is from optimal.

   primal objective values             dual objective values
   (feasible x, maximize)              (feasible y, minimize)

   ------------------------->  z*  <------------------------
   0      12     27     30     36     37.6     45     ...
          |<---- gap for x=(4,3), y=(0,1.5,1.5): 27 vs 45 ---->|

Example: random feasible points. Sampling feasible points of both problems shows that every primal value lies below every dual value:

import numpy as np

c = np.array([3, 5])
A = np.array([[1, 0], [0, 2], [3, 2]])
b = np.array([4, 12, 18])
rng = np.random.default_rng(0)

best_primal, best_dual = -np.inf, np.inf
for _ in range(2000):
    x = rng.uniform(0, 7, size=2)
    if np.all(A @ x <= b):                  # primal feasible -> a lower bound on z*
        best_primal = max(best_primal, c @ x)
    y = rng.uniform(0, 4, size=3)
    if np.all(A.T @ y >= c):                # dual feasible -> an upper bound on z*
        best_dual = min(best_dual, b @ y)

print(f"best primal value found: {best_primal:.2f}")    # 35.47
print(f"best dual value found  : {best_dual:.2f}")      # 37.59
print(f"so z* lies in [{best_primal:.2f}, {best_dual:.2f}]; true value 36")

Without solving anything, random sampling already brackets the optimum within about 6%.

Key takeaways

• Any primal feasible value $\le$ any dual feasible value (for a max primal).
• Equal values for a feasible pair certify that both solutions are optimal.
• If one problem is unbounded, the other is infeasible.
• The duality gap $b^\top y - c^\top x$ measures suboptimality and is the standard stopping criterion in solvers.

Practice

1. (Beginner) Verify weak duality for $x = (4, 3)$ and $y = (0, 1.5, 1.5)$. Is $y$ dual feasible? What is the gap?
2. (Intermediate) Write the weak duality proof for a minimization primal $\min{c^\top x : Ax \ge b, x \ge 0}$. Which inequalities change direction?
3. (Interview) "A solver stops with 'relative gap 0.01%'. What exactly does it guarantee?" Hint: the reported primal solution is within 0.01% of a proven bound, regardless of the true optimum.

Strong Duality Theorem

Why it matters. Weak duality gives bounds. Strong duality says the best bound is exact: there is never a gap between the optimal values. This makes duality an exact alternative view of the same problem. You can solve whichever of the two is easier, and the dual solution comes for free with every primal optimum.

Theory. Theorem. If the primal has an optimal solution $x^$, then the dual has an optimal solution $y^$, and

$$ c^\top x^* = b^\top y^* . $$

Proof idea (via simplex). At an optimal basis $B$ in standard form, all reduced costs satisfy $c_j - c_B^\top B^{-1}A_j \le 0$ (maximization). Set $y^\top = c_B^\top B^{-1}$. Then $A^\top y \ge c$, so $y$ is dual feasible (the slack columns give $y \ge 0$), and $b^\top y = c_B^\top B^{-1}b = c^\top x^*$. By weak duality, $y$ is optimal. (A proof without simplex uses Farkas' lemma, below.)

Combined with weak duality, only four of the nine status combinations can occur:

Primal \ Dual	Optimal	Unbounded	Infeasible
Optimal	Possible (equal values)	Impossible	Impossible
Unbounded	Impossible	Impossible	Possible
Infeasible	Impossible	Possible	Possible

Both problems can be infeasible at once, e.g. $\max{x_1 : x_1 - x_2 \le -1,\ -x_1 + x_2 \le -1,\ x \ge 0}$.

Example: duals from a solver. Solvers return the dual solution together with the primal one. SciPy's HiGHS interface reports marginals: the derivative of its (minimization) objective with respect to each right-hand side.

import numpy as np
from scipy.optimize import linprog

c = np.array([3, 5])
A = np.array([[1, 0], [0, 2], [3, 2]])
b = np.array([4, 12, 18])

res = linprog(-c, A_ub=A, b_ub=b, method="highs")   # max c^T x  ==  min -c^T x
x = res.x
y = -res.ineqlin.marginals      # marginals of the minimization form; negate for the max form

print("primal x* =", x, " c^T x* =", c @ x)       # [2. 6.]   36.0
print("dual   y* =", y + 0.0, " b^T y* =", b @ y) # [0.  1.5 1. ]  36.0
print("dual feasible:", np.all(A.T @ y >= c - 1e-9), np.all(y >= 0))   # True True

Sign conventions differ between solvers (and between max/min forms), so always sanity-check by verifying $b^\top y^* = c^\top x^*$.

Key takeaways

• If either problem has a finite optimum, both do, and the optimal values are equal (zero duality gap).
• The optimal dual is $y^{*\top} = c_B^\top B^{-1}$ from the optimal simplex basis, so solvers obtain it at no extra cost.
• Only four status combinations are possible: both optimal, (unbounded, infeasible), (infeasible, unbounded), both infeasible.
• Check solver sign conventions by verifying $b^\top y^* = c^\top x^*$ and dual feasibility.

Practice

1. (Beginner) Verify by hand that $y^* = (0, 1.5, 1)$ satisfies both dual constraints and gives value 36.
2. (Intermediate) Solve the bakery dual with linprog directly and read the marginals of its constraints. Which primal quantities do you recover?
3. (Interview) "Strong duality holds for LP. Does it hold for integer programs?" Hint: compare the integer optimum with the LP relaxation's dual bound in the classification example of Chapter 1.3 (40 vs 41.25).

Complementary Slackness

Why it matters. Complementary slackness turns optimality into a simple check: "a resource with leftover capacity is worth nothing, and a product that is made earns exactly what its resources are worth". It lets you verify a proposed solution without solving the LP, recover the dual from the primal (or the reverse), and it is the basis of primal-dual and interior point algorithms.

Theory. For the symmetric pair, feasible $x$ and $y$ are both optimal if and only if

$$ \underbrace{y_i ,\big(b_i - a_i^\top x\big) = 0 \ \ \forall i}{\text{primal slack} \times \text{dual variable}} \qquad \text{and} \qquad \underbrace{x_j ,\big(A_j^\top y - c_j\big) = 0 \ \ \forall j}{\text{primal variable} \times \text{dual slack}} . $$

In words:

• If a constraint is not binding (slack $> 0$), its dual price is 0.
• If a dual price is positive, its constraint is binding.
• If a variable is positive, its dual constraint is tight (profit = imputed cost of resources).
• If a dual constraint is slack (resources are worth more than the profit), the variable is 0.

Why: in the weak duality chain $c^\top x \le y^\top Ax \le y^\top b$, equality (optimality) requires both inequalities to be equalities. Each is a sum of nonnegative products, so every product must be zero.

Example: verifying and rejecting candidate plans.

import numpy as np

c = np.array([3, 5])
A = np.array([[1, 0], [0, 2], [3, 2]], dtype=float)
b = np.array([4, 12, 18], dtype=float)

def certify(x):
    """Use complementary slackness to build a dual y for x, then check dual feasibility."""
    x = np.asarray(x, dtype=float)
    slack = b - A @ x
    binding = np.isclose(slack, 0)                   # non-binding row -> y_i = 0
    tight = x > 1e-9                                 # positive x_j -> (A^T y)_j = c_j
    # Unknowns: y_i of the binding rows. Equations: the tight dual constraints.
    free = np.where(binding)[0]
    M = A[:, tight].T[:, free]
    y = np.zeros(len(b))
    y[free] = np.linalg.lstsq(M, c[tight], rcond=None)[0]
    ok = np.all(y >= -1e-9) and np.all(A.T @ y >= c - 1e-9)
    print(f"x = {x}: slack = {slack}, derived y = {np.round(y, 3) + 0.0}, "
          f"{'OPTIMAL' if ok else 'not optimal'}")

certify([2, 6])   # slack [2 0 0]  y = [0, 1.5, 1]    -> dual feasible -> OPTIMAL
certify([4, 3])   # slack [0 6 0]  y = [-4.5, 0, 2.5] -> y1 < 0        -> not optimal

For $(4, 3)$, complementary slackness forces $y_2 = 0$ (cake oven has slack), and making both products requires $y_1 + 3y_3 = 3$ and $2y_3 = 5$, so $y_1 = -4.5 < 0$. No valid dual exists, which proves $(4, 3)$ is not optimal.

Key takeaways

• Optimality $\iff$ primal feasibility + dual feasibility + complementary slackness.
• Slack constraint $\Rightarrow$ zero price. Positive variable $\Rightarrow$ tight dual constraint.
• Given a candidate primal solution, complementary slackness determines the dual (under nondegeneracy) and checks optimality without solving an LP.
• Interior point methods solve exactly these conditions, with $x_j s_j = \mu$ relaxed toward 0 (Section 2.6).

Practice

1. (Beginner) For the furniture workshop of Section 2.1 (solution: 2 desks, 0 tables, 8 chairs), which dual prices must be zero according to complementary slackness?
2. (Intermediate) Use certify on the vertex $(0, 6)$. Explain the result in terms of the bread product.
3. (Interview) "Complementary slackness says $x_j > 0 \Rightarrow$ the dual constraint is tight. Is the converse true, i.e. tight dual constraint $\Rightarrow x_j > 0$?" Hint: consider degenerate solutions and multiple optima. Strict complementarity is a separate, stronger property.

Economic Interpretation of Dual Variables

Why it matters. Dual variables turn an LP from a calculator into a decision-support tool. They tell managers what an extra hour of machine time, an extra ton of material or a relaxed regulation is worth, which is the information needed for investment, pricing and negotiation decisions.

Theory. For a resource-allocation primal ($\max$ profit, $\le$ resources):

• $y_i$ is the marginal value (shadow price) of resource $i$: the increase in optimal profit per extra unit of $b_i$, valid for small changes (Section 2.5 computes the valid range).
• The dual problem is the resource buyer's problem. An outside firm wants to buy all your resources. It offers price $y_i$ per unit of resource $i$ and wants to pay as little as possible, $\min, b^\top y$. You only sell if, for every product, the resources it uses are worth at least the profit it would earn: $A_j^\top y \ge c_j$. Strong duality says the fair total price equals your best production profit.
• $A_j^\top y$ is the imputed cost of product $j$ (its resources valued at shadow prices). Its reduced cost is $c_j - A_j^\top y$.

Situation	Interpretation
$y_i = 0$	Resource $i$ is not scarce. More of it is worthless at the margin
$y_i > 0$	Resource $i$ is a bottleneck. Pay up to $y_i$ per extra unit
$c_j = A_j^\top y$	Product $j$ breaks even at shadow prices, and it can be in the plan
$c_j < A_j^\top y$	Product $j$ uses resources worth more than it earns. Do not make it

Analogy: shadow prices are the "internal exchange rates" of scarce resources. A product is worth making only if it converts resources into profit at least at those rates.

Example: furniture workshop. Shadow prices explain why tables are not made:

import numpy as np
from scipy.optimize import linprog

profit = np.array([60, 30, 20])                   # desk, table, chair
A = np.array([[8, 6, 1], [4, 2, 1.5], [2, 1.5, 0.5]])
b = np.array([48, 20, 8])                          # lumber, finishing, carpentry

res = linprog(-profit, A_ub=A, b_ub=b, bounds=[(0, None), (0, 5), (0, None)])
y = -res.ineqlin.marginals + 0.0
print("shadow prices (lumber, finishing, carpentry):", y)      # [ 0. 10. 10.]

for name, j in [("desk", 0), ("table", 1), ("chair", 2)]:
    imputed = A[:, j] @ y                          # value of the resources one unit uses
    print(f"{name:<6} profit {profit[j]:>3}  imputed cost {imputed:>5.1f}  "
          f"reduced cost {profit[j] - imputed:>5.1f}")
# desk   profit  60  imputed cost  60.0  reduced cost   0.0   <- breaks even, produced
# table  profit  30  imputed cost  35.0  reduced cost  -5.0   <- loses 5 at shadow prices
# chair  profit  20  imputed cost  20.0  reduced cost   0.0   <- breaks even, produced

print("resource buyer pays:", b @ y, "= optimal profit:", -res.fun)   # 280.0 = 280.0

A table uses 2 finishing hours and 1.5 carpentry hours worth $2 \cdot 10 + 1.5 \cdot 10 = 35$, but earns only 30. Lumber has a shadow price of 0 because 24 board feet are unused.

Key takeaways

• Dual variables are shadow prices: the marginal value of one more unit of each constraint's right-hand side.
• The dual is the resource buyer's problem, and strong duality makes the fair resource price equal to the optimal profit.
• Imputed cost $A_j^\top y$ values a product's resource consumption. Reduced cost $c_j - A_j^\top y$ tells whether it pays to make it.
• Non-binding resources have zero value. Investment should target resources with high shadow prices.

Practice

1. (Beginner) A supplier offers extra mixer hours for the bakery at 0.8 per hour. Should you buy? What about 1.2 per hour?
2. (Intermediate) In the furniture example, at what price would renting extra finishing hours stop being worthwhile? Verify by adding a "rent finishing hour" variable to the LP with that cost.
3. (Interview) "A plant manager wants to double the capacity of the machine with the highest utilization. How do you respond?" Hint: utilization is not value. Look at the shadow price, and remember it is only valid within a range (Section 2.5).

Dual Simplex Method

Why it matters. The dual simplex method is the workhorse for re-optimization. After adding a constraint (a cutting plane, a branching bound in Chapter 4) or changing right-hand sides, the old optimal basis is still dual feasible but no longer primal feasible, and the dual simplex method repairs it in a few pivots instead of starting over. Many commercial solvers use it as their default LP algorithm.

Theory. Compare the two methods:

	Primal simplex	Dual simplex
Maintains	Primal feasibility (rhs $\ge 0$)	Dual feasibility (z-row $\ge 0$ for max)
Works toward	Dual feasibility (optimality)	Primal feasibility
Chooses first	Entering variable (improving z-row entry)	Leaving variable (negative rhs)
Then	Ratio test on positive column entries	Ratio test on negative row entries
Objective	Improves monotonically	Worsens monotonically toward $z^*$ (bound tightens)

Dual simplex iteration (tableau in z-row form, maximization, all z-row entries $\ge 0$):

1. Leaving row: pick a row $r$ with negative rhs (most negative, for instance). If none, the tableau is optimal.
2. Entering column: among columns $j$ with $\bar a_{rj} < 0$, choose $q$ minimizing $\left|\bar z_j / \bar a_{rj}\right|$ (the dual ratio test). This keeps every z-row entry nonnegative. If no $\bar a_{rj} < 0$, the primal is infeasible.
3. Pivot on $(r, q)$ and repeat.

The dual simplex method on the primal performs the same pivots as the primal simplex method would on the dual. It just works on the primal tableau.

Example.

$$ \min\ 3x_1 + 2x_2 \quad \text{s.t.} \quad 3x_1 + x_2 \ge 3,\quad 4x_1 + 3x_2 \ge 6,\quad x_1 + x_2 \le 3,\quad x \ge 0 . $$

Multiplying the $\ge$ rows by $-1$ and adding slacks gives an all-slack basis with negative rhs (primal infeasible), while the z-row $(3, 2)$ is already nonnegative (dual feasible). This is the ideal starting point for the dual simplex method, and no artificial variables are needed.

from fractions import Fraction as F

# Reuses show() and pivot() from the Simplex Tableau topic (Section 2.3).

def dual_simplex(T, basis, names):
    """T must be dual feasible (z-row >= 0, maximization form)."""
    it = 0
    while True:
        r = min(range(1, len(T)), key=lambda i: T[i][-1])       # most negative rhs
        if T[r][-1] >= 0:
            return "optimal"
        candidates = [j for j in range(len(T[0]) - 1) if T[r][j] < 0]
        if not candidates:
            return "infeasible"            # entries >= 0 but rhs < 0: no x >= 0 fits this row
        q = min(candidates, key=lambda j: abs(T[0][j] / T[r][j]))   # dual ratio test
        it += 1
        print(f"iteration {it}: {names[basis[r - 1]]} leaves (rhs {T[r][-1]}), {names[q]} enters")
        pivot(T, basis, r, q)

names = ["x1", "x2", "e1", "e2", "s3"]
T = [[F(v) for v in row] for row in [
    [ 3,  2, 0, 0, 0,  0],    # z-row of max -3x1 - 2x2: z + 3x1 + 2x2 = 0  (dual feasible)
    [-3, -1, 1, 0, 0, -3],    # -3x1 -  x2 + e1 = -3   (from 3x1 + x2 >= 3)
    [-4, -3, 0, 1, 0, -6],    # -4x1 - 3x2 + e2 = -6
    [ 1,  1, 0, 0, 1,  3]]]   #   x1 +  x2 + s3 =  3
basis = [2, 3, 4]

print(dual_simplex(T, basis, names))
show(T, basis, names)
print("min cost =", -T[0][-1])

iteration 1: e2 leaves (rhs -6), x2 enters
iteration 2: e1 leaves (rhs -1), x1 enters
optimal
basis |     x1     x2     e1     e2     s3 |    rhs
    z |      0      0    1/5    3/5      0 |  -21/5
   x1 |      1      0   -3/5    1/5      0 |    3/5
   x2 |      0      1    4/5   -3/5      0 |    6/5
   s3 |      0      0   -1/5    2/5      1 |    6/5
min cost = 21/5

Two pivots reach the optimum $x = (3/5, 6/5)$ with cost $21/5$. The z-row entries under $e_1$ and $e_2$ ($1/5$ and $3/5$) are the dual values of the two $\ge$ constraints.

Key takeaways

• Dual simplex keeps optimality (dual feasibility) and works toward primal feasibility. The leaving variable is chosen first, then a dual ratio test picks the entering one.
• It needs no artificial variables when the initial z-row is already optimal, e.g. a min-cost problem with $c \ge 0$ and $\ge$ constraints.
• It is the standard method for re-optimizing after adding constraints or changing $b$, and is central in branch-and-bound and cutting planes.
• A negative-rhs row with no negative entries proves primal infeasibility.

Practice

1. (Beginner) Perform the first dual simplex pivot of the example by hand and check the intermediate tableau.
2. (Intermediate) Take the final bakery tableau, add the constraint $x_1 + x_2 \le 7$ as a new row (expressed in terms of the current nonbasic variables), and re-optimize with dual_simplex.
3. (Interview) "Why is dual simplex preferred inside branch-and-bound?" Hint: a branching bound like $x_j \le 2$ cuts off the current LP optimum but leaves the z-row unchanged.

Farkas' Lemma

Why it matters. Farkas' lemma is the "theorem of the alternative" underneath LP duality. It guarantees that when a system of linear constraints is infeasible, there is always a short, checkable proof of infeasibility. Solvers return such certificates (infeasibility rays), and modelers use them to locate contradictory constraints.

Theory. Farkas' lemma. For $A \in \mathbb{R}^{m \times n}$ and $b \in \mathbb{R}^m$, exactly one of the following systems has a solution:

$$ \text{(I)}\quad Ax = b,\ x \ge 0 \qquad\qquad \text{(II)}\quad A^\top y \ge 0,\ b^\top y < 0 . $$

Both cannot hold: if they did, $0 > b^\top y = (Ax)^\top y = x^\top (A^\top y) \ge 0$, a contradiction. The deep part is that one always holds. This is proven with a separating hyperplane argument.

Geometric meaning. ${Ax : x \ge 0}$ is the cone generated by the columns of $A$. Either $b$ lies in that cone (I), or there is a hyperplane through the origin, with normal $y$, that separates $b$ from the cone (II): all columns are on one side ($A^\top y \ge 0$) and $b$ is strictly on the other ($b^\top y < 0$).

            A_2
             ^
             |    cone {A x : x >= 0}
             |  ////////
             | ////////
             |////////
             +-----------> A_1
            /
           /   y  (normal of separating hyperplane)
          v
         b      <- b outside the cone: certificate y with A^T y >= 0, b^T y < 0

Inequality version (for $Ax \le b$, $x \ge 0$): either a feasible $x$ exists, or there is $y \ge 0$ with $A^\top y \ge 0$ and $b^\top y < 0$. The certificate $y$ is a nonnegative combination of the constraints that reads "$0 \le (\text{negative number})$", which is impossible.

Use in duality. Strong duality can be proven from Farkas: if no $y$ achieved the bound $c^\top x^$, the resulting system would be infeasible, and its Farkas certificate would contradict the optimality of $x^$.

import numpy as np
from scipy.optimize import linprog

# Infeasible system:  x1 + x2 <= 2   and   x1 + x2 >= 3  (written as -x1 - x2 <= -3),  x >= 0
A = np.array([[1, 1], [-1, -1]], dtype=float)
b = np.array([2, -3], dtype=float)

primal = linprog(np.zeros(2), A_ub=A, b_ub=b)
print("primal status:", primal.status)               # 2 = infeasible

# Search for a Farkas certificate: y >= 0, A^T y >= 0, b^T y <= -1 (normalized)
cert = linprog(np.zeros(2),
               A_ub=np.vstack([-A.T, b]),            # -A^T y <= 0  and  b^T y <= -1
               b_ub=np.r_[np.zeros(2), -1],
               bounds=[(0, None)] * 2)
y = cert.x
print("certificate y =", y)                          # [1. 1.]
print("A^T y =", A.T @ y, " b^T y =", b @ y)         # [0. 0.]   -1.0
# Read it: 1*(x1 + x2 <= 2) + 1*(-x1 - x2 <= -3)  gives  0 <= -1.  Impossible.

Key takeaways

• Exactly one holds: $Ax = b, x \ge 0$ is feasible, or there is $y$ with $A^\top y \ge 0$, $b^\top y < 0$.
• Geometrically, $b$ is either in the cone of the columns of $A$ or separated from it by a hyperplane.
• An infeasible LP always has a certificate, a combination of constraints that gives a contradiction such as $0 \le -1$.
• Farkas' lemma underlies strong duality, Phase I infeasibility proofs and solver "infeasibility rays".

Practice

1. (Beginner) Verify that $y = (1, 1)$ is a valid certificate for the example by writing out the weighted sum of the two constraints.
2. (Intermediate) Find a Farkas certificate for $x_1 + 2x_2 \le 4$, $x_1 \ge 3$, $x_2 \ge 1$, $x \ge 0$ (write all rows as $\le$), first by hand, then with the code.
3. (Interview) "Your model is infeasible and has 10,000 constraints. How can a Farkas certificate help you find the problem?" Hint: constraints with $y_i = 0$ play no role in the contradiction, so the support of $y$ points to a small conflicting subset (related to an IIS).

<a id="25-sensitivity-and-post-optimality-analysis"></a>

2.5 Sensitivity and Post-Optimality Analysis

An optimal solution is only as good as the data behind it. This subchapter shows how to answer "what if?" questions from the final simplex basis, without re-solving: how the optimum responds to changes in prices, capacities, new products and new rules, and over which ranges those answers remain valid.

All results here come from one object, the optimal basis $B$. For the bakery, $B$ consists of the columns of $(s_1, x_2, x_1)$:

$$ B = \begin{bmatrix} 1 & 0 & 1 \ 0 & 2 & 0 \ 0 & 2 & 3 \end{bmatrix}, \qquad B^{-1} = \begin{bmatrix} 1 & 1/3 & -1/3 \ 0 & 1/2 & 0 \ 0 & -1/3 & 1/3 \end{bmatrix}, \qquad x_B = B^{-1}b = \begin{bmatrix} 2 \ 6 \ 2 \end{bmatrix}, \qquad y^\top = c_B^\top B^{-1} = (0,\ 1.5,\ 1). $$

Change	Affects	Basis stays optimal while...
Objective coefficient $c_j$	Reduced costs (optimality)	All reduced costs keep their sign
Right-hand side $b_i$	$x_B = B^{-1}b$ (feasibility)	$x_B \ge 0$
New variable	One new reduced cost	Its reduced cost is not improving
New constraint	Feasibility of current $x^*$	$x^*$ satisfies it

Shadow Prices

Why it matters. The most frequent managerial question about an LP is "what is one more unit of this resource worth?". The shadow price answers it directly. It sets the maximum price worth paying for extra capacity, the cost of a regulation, or the value of relaxing a contract.

Theory. The shadow price of constraint $i$ is the rate of change of the optimal value with respect to its right-hand side:

$$ y_i^* = \frac{\partial z^}{\partial b_i}, \qquad \text{so } \ z^(b + \Delta e_i) = z^* + y_i^* \Delta \quad \text{for } \Delta \text{ within the RHS range.} $$

By strong duality, shadow prices equal the optimal dual variables, $y^{*\top} = c_B^\top B^{-1}$. In the final tableau, they are the z-row entries under the slack columns.

Problem	Constraint	Sign of the shadow price	Meaning
max	$\le$	$\ge 0$	More resource can only help
max	$\ge$	$\le 0$	Raising a requirement can only hurt
min	$\ge$	$\ge 0$	Raising a requirement raises cost
min	$\le$	$\le 0$	More capacity can only lower cost
any	$=$	any sign

Shadow prices are marginal values. They hold only while the optimal basis stays the same, i.e. inside the RHS ranging interval (a later topic). Beyond it, the value of an extra unit changes, usually decreasing (diminishing returns).

Analogy: a shadow price is like a store price for "one more unit". It is accurate for a few units, but if you want to buy in bulk, the price changes.

import numpy as np
from scipy.optimize import linprog

c = np.array([3, 5])
A = np.array([[1, 0], [0, 2], [3, 2]])
b = np.array([4, 12, 18], dtype=float)
names = ["bread oven", "cake oven", "mixer"]

def solve(b):
    res = linprog(-c, A_ub=A, b_ub=b)
    return -res.fun, -res.ineqlin.marginals + 0.0

z0, y = solve(b)
print("shadow prices:", y)                           # [0.  1.5 1. ]
for i, name in enumerate(names):
    for delta in (1, 10):
        b_new = b.copy()
        b_new[i] += delta
        gain = solve(b_new)[0] - z0
        print(f"{name:<10} +{delta:<2}: profit +{gain:<4g} (shadow price predicts +{y[i] * delta:g})")
# bread oven +1 : profit +0    (shadow price predicts +0)
# bread oven +10: profit +0    (shadow price predicts +0)
# cake oven  +1 : profit +1.5  (shadow price predicts +1.5)
# cake oven  +10: profit +9    (shadow price predicts +15)   <- outside the valid range
# mixer      +1 : profit +1    (shadow price predicts +1)
# mixer      +10: profit +6    (shadow price predicts +10)   <- outside the valid range

An extra cake-oven hour is worth 1.5 (hundred EUR), an extra mixer hour 1, an extra bread-oven hour nothing. Ten extra hours are worth less than ten times one hour, because the basis changes along the way.

Key takeaways

• Shadow price = marginal change in the optimal value per unit increase of a right-hand side = optimal dual variable.
• Non-binding constraints have shadow price 0. Binding constraints usually have nonzero prices.
• Signs depend on max/min and the constraint direction. Check the solver's convention.
• Shadow prices are valid only within the RHS range. Larger changes need re-solving or parametric analysis.

Practice

1. (Beginner) Read the shadow prices from the final bakery tableau in Section 2.3 and compare with the code output.
2. (Intermediate) For the diet problem in Section 2.1, compute the shadow prices of the four nutrient requirements and interpret the price of the sugar constraint in cents per ounce.
3. (Interview) "The shadow price of a constraint is zero. Can we remove the constraint?" Hint: zero means it is not binding at this optimum. After other data changes it may become binding, and degenerate optima can also have zero prices on binding constraints.

Reduced Costs

Why it matters. Reduced costs answer the mirror question for variables: "why is this product, route or investment not used, and what would it take to use it?". They give a precise target, such as a price increase or a cost reduction, rather than a vague "it is not profitable".

Theory. The reduced cost of variable $x_j$ is

$$ \bar c_j = c_j - y^{*\top} A_j = \text{(direct profit)} - \text{(value of the resources it consumes at shadow prices)} . $$

At an optimal solution of a max problem:

• Basic variables have $\bar c_j = 0$.
• Nonbasic variables at 0 have $\bar c_j \le 0$. Forcing one unit of $x_j$ into the solution changes the optimal value by $\bar c_j$ (for small amounts).
• $|\bar c_j|$ is how much $c_j$ must improve before $x_j$ becomes attractive. It is the allowable increase of a nonbasic cost coefficient.
• A nonbasic variable with $\bar c_j = 0$ signals alternative optima.

For a min problem, all signs flip: nonbasic reduced costs are $\ge 0$, and $\bar c_j$ is how much the cost must decrease.

Variable status	Reduced cost (max)	Interpretation
Basic (in the solution)	$0$	Earns exactly its resource value
Nonbasic at lower bound	$\le 0$	Loses $\lvert \bar c_j \rvert$ per unit if forced in
Nonbasic at upper bound	$\ge 0$	Would gain if its upper bound were raised

import numpy as np
from scipy.optimize import linprog

profit = np.array([60, 30, 20])                         # desk, table, chair
A = np.array([[8, 6, 1], [4, 2, 1.5], [2, 1.5, 0.5]])
b = np.array([48, 20, 8])
bounds = [(0, None), (0, 5), (0, None)]

res = linprog(-profit, A_ub=A, b_ub=b, bounds=bounds)
# lower.marginals = d(min objective)/d(lower bound) = reduced costs of the min form
reduced = -res.lower.marginals + 0.0                     # back to the max form
print("plan          :", res.x)                         # [2. 0. 8.]
print("reduced costs :", reduced)                       # [ 0. -5.  0.] -> tables lose 5 each

for table_profit in (34, 36):                           # just below / above 30 + 5
    p = profit.copy()
    p[1] = table_profit
    r = linprog(-p, A_ub=A, b_ub=b, bounds=bounds)
    print(f"table profit {table_profit}: plan {np.round(r.x, 2)}, profit {-r.fun:.1f}")
# table profit 34: plan [2. 0. 8.], profit 280.0       <- still not worth it
# table profit 36: plan [ 0.   1.6 11.2], profit 281.6  <- tables enter the plan

The reduced cost of $-5$ says exactly what the sales team needs to know: tables need at least 5 more profit per unit, through a price increase or cost reduction, before they belong in the plan.

Key takeaways

• Reduced cost = objective coefficient minus the shadow-price value of the resources the variable uses.
• Basic variables have reduced cost 0. Nonbasic ones show how much their coefficient must improve to enter.
• Forcing a nonbasic variable into the solution worsens the objective by about $|\bar c_j|$ per unit.
• A zero reduced cost on a nonbasic variable means alternative optimal solutions exist.

Practice

1. (Beginner) In the final bakery tableau, what are the reduced costs of $x_1$, $x_2$, $s_2$ and $s_3$? Interpret the one for $s_3$.
2. (Intermediate) Force the furniture workshop to make one table (bounds for tables = (1, 5)). Does the profit drop by exactly 5? Explain.
3. (Interview) "In a 5000-variable model, how would you use reduced costs to decide which variables to drop to speed things up?" Hint: variables with very unattractive reduced costs are unlikely to enter, but beware of how the duals change when data changes. Look up "variable fixing by reduced costs".

Ranging of Objective Coefficients

Why it matters. Profits and costs are estimates. Ranging tells you how wrong each estimate can be before the optimal plan changes. If the range is wide, the plan is robust and you can stop worrying about that number. If it is narrow, the estimate deserves more effort.

Theory. Change one coefficient $c_j \to c_j + \delta$. The feasible region does not change, so the current BFS stays feasible. The question is whether it stays optimal, i.e. whether all reduced costs keep their sign.

• Nonbasic $x_j$: only $\bar c_j$ changes, to $\bar c_j + \delta$. For max, the basis stays optimal while $\delta \le -\bar c_j = |\bar c_j|$, and $\delta$ can decrease without limit.
• Basic $x_j$ (row $k$ of the final tableau): the change spreads to every nonbasic reduced cost, $\bar c_q(\delta) = \bar c_q - \delta, \bar a_{kq}$. The basis stays optimal while all $\bar c_q(\delta) \le 0$:

$$ \max_{q:\ \bar a_{kq} > 0} \frac{\bar c_q}{\bar a_{kq}} ;\le; \delta ;\le; \min_{q:\ \bar a_{kq} < 0} \frac{\bar c_q}{\bar a_{kq}} . $$

Inside the range, $x^$ is unchanged and the optimal value changes linearly: $z^ + \delta x_j^*$.

Geometrically, changing $c_j$ rotates the iso-profit lines. The current vertex stays optimal until the lines become parallel to one of the edges leaving it.

  Rotating the objective around the bakery optimum (2, 6):

  x2              slope of iso-profit line = -c1/5
  6 *=====*          c1 = 0   : parallel to cake-oven edge   (lower limit)
    |      \         c1 = 3   : current
    |       \        c1 = 7.5 : parallel to mixer edge        (upper limit)
  3 +        *
    |        |       inside [0, 7.5] the optimum stays at (2, 6)
    *--------*---> x1

import numpy as np
from scipy.optimize import linprog

c = np.array([3, 5, 0, 0, 0], dtype=float)                    # max form, standard form
A = np.array([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=float)
names = ["x1", "x2", "s1", "s2", "s3"]
basis = [2, 1, 0]                                             # optimal basis (s1, x2, x1)

Binv = np.linalg.inv(A[:, basis])
tableau = Binv @ A                                            # rows of the final tableau
cbar = c - (c[basis] @ Binv) @ A                              # reduced costs (<= 0: optimal)
nonbasic = [j for j in range(5) if j not in basis]

for k, j in enumerate(basis):
    if names[j].startswith("s"):
        continue                                              # only decision variables
    lo, hi = -np.inf, np.inf
    for q in nonbasic:
        a = tableau[k, q]
        if a > 1e-12:
            lo = max(lo, cbar[q] / a)
        elif a < -1e-12:
            hi = min(hi, cbar[q] / a)
    print(f"c_{names[j]} = {c[j]:g}: stays optimal for c in [{c[j] + lo:g}, {c[j] + hi:g}]")
# c_x2 = 5: stays optimal for c in [2, inf]
# c_x1 = 3: stays optimal for c in [0, 7.5]

# Brute-force check: re-solve over a grid of bread profits
A_ub, b_ub = [[1, 0], [0, 2], [3, 2]], [4, 12, 18]
for c1 in (-0.5, 0.5, 7.4, 7.6):
    x = linprog([-c1, -5], A_ub=A_ub, b_ub=b_ub).x
    print(f"c1 = {c1:>4}: x* = {x}")
# c1 = -0.5: x* = [0. 6.]   <- below the range, the optimum moves
# c1 =  0.5: x* = [2. 6.]
# c1 =  7.4: x* = [2. 6.]
# c1 =  7.6: x* = [4. 3.]   <- above the range, the optimum moves

This reproduces exactly the threshold of 7.5 observed by brute force in Chapter 1.3, now computed from one tableau with no re-solving.

Key takeaways

• Changing $c_j$ never affects feasibility, only optimality (the signs of reduced costs).
• Nonbasic variable: allowable increase $= |\bar c_j|$, allowable decrease unlimited (for max).
• Basic variable: the range follows from ratios of reduced costs to the variable's row in the final tableau.
• Inside the range the plan is unchanged and the profit changes linearly. At the limits, alternative optima appear.

Practice

1. (Beginner) Using the ranges, decide without solving whether the plan changes if the cake profit drops from 5 to 2.5.
2. (Intermediate) Compute the ranges for desks, tables and chairs in the furniture workshop with the same code (add slack columns, find the optimal basis from the linprog solution).
3. (Interview) "Two profit coefficients are uncertain at the same time. Can you use the individual ranges?" Hint: look up the 100% rule: the sum of each change divided by its allowable change (in that direction) must not exceed 1 for a guarantee.

Ranging of Right-Hand Side Values

Why it matters. Capacities, budgets and demands also change: a machine breaks down, a budget is cut, a contract grows. RHS ranging tells you how far each right-hand side can move while the current shadow prices and the current set of binding constraints remain valid, which is the region where "each extra hour is worth $y_i$" is true.

Theory. Change $b_i \to b_i + \Delta$. Now the reduced costs are unaffected (they do not depend on $b$), so the basis stays optimal as long as it stays feasible:

$$ x_B(\Delta) = B^{-1}(b + \Delta e_i) = x_B + \Delta, (B^{-1})_{\cdot i} \ \ge\ 0 . $$

With $d = (B^{-1})_{\cdot i}$ (column $i$ of $B^{-1}$):

$$ \max_{k:\ d_k > 0} \left(-\frac{x_{B_k}}{d_k}\right) ;\le; \Delta ;\le; \min_{k:\ d_k < 0} \left(-\frac{x_{B_k}}{d_k}\right) . $$

Inside the range:

• The same constraints are binding, and the optimal values of the variables change linearly in $\Delta$.
• The optimal objective changes by exactly $y_i^* \Delta$.

At the limit, a basic variable hits 0 (degenerate vertex). Beyond it, a dual simplex pivot changes the basis and the shadow price.

import numpy as np

A = np.array([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=float)
b = np.array([4, 12, 18], dtype=float)
c = np.array([3, 5, 0, 0, 0], dtype=float)
basis = [2, 1, 0]                                  # optimal basis (s1, x2, x1)
names = ["bread oven", "cake oven", "mixer"]

Binv = np.linalg.inv(A[:, basis])
x_B = Binv @ b                                      # [2, 6, 2]
y = c[basis] @ Binv                                 # shadow prices [0, 1.5, 1]

for i, name in enumerate(names):
    d = Binv[:, i]                                  # how x_B moves per unit of b_i
    lo = max((-x_B[k] / d[k] for k in range(3) if d[k] > 1e-12), default=-np.inf)
    hi = min((-x_B[k] / d[k] for k in range(3) if d[k] < -1e-12), default=np.inf)
    print(f"{name:<10} b = {b[i]:>4g}  range [{b[i] + lo:g}, {b[i] + hi:g}]  "
          f"shadow price {y[i]:g}")
# bread oven b =    4  range [2, inf]  shadow price 0
# cake oven  b =   12  range [6, 18]  shadow price 1.5
# mixer      b =   18  range [12, 24]  shadow price 1

Why the mixer range is $[12, 24]$: below 12 hours, the mixer cannot even support 6 cakes, and the plan stops making bread. Above 24 hours, bread production hits the bread-oven limit $x_1 = 4$, and extra mixer time becomes worthless. The earlier finding that 10 extra mixer hours are worth only 6, not 10, is explained by this range: only 6 of those 10 hours fall inside it.

Key takeaways

• Changing $b_i$ never affects optimality (reduced costs), only feasibility ($x_B = B^{-1}b \ge 0$).
• The range comes from column $i$ of $B^{-1}$: how much $b_i$ can move before some basic variable hits 0.
• Inside the range, the shadow price is constant and $z^*$ changes linearly. The plan values do change.
• Non-binding constraints can increase without limit, and can decrease down to their current usage.

Practice

1. (Beginner) Explain in words why the bread-oven range is $[2, \infty)$.
2. (Intermediate) For $b_3 = 21$ (inside the range), compute the new plan as $x_B + 3 \cdot (B^{-1})_{\cdot 3}$ and verify it with linprog.
3. (Interview) "Your sensitivity report says a capacity can increase by 'infinity' at shadow price 0. A manager reads this as 'capacity doesn't matter'. Is that right?" Hint: it doesn't matter upward. Check the allowable decrease.

Adding New Variables or Constraints

Why it matters. Models evolve. Marketing proposes a new product, engineering adds a safety rule, a new supplier appears. Re-solving from scratch wastes the information in the current optimal basis. Post-optimality analysis evaluates the change in one step and re-optimizes from a warm start in a few pivots.

Theory.

New variable $x_{n+1}$ with profit $c_{n+1}$ and column $A_{n+1}$:

1. Price it out using the current duals: $\bar c_{n+1} = c_{n+1} - y^{*\top} A_{n+1}$.
2. If $\bar c_{n+1} \le 0$ (max), the current solution stays optimal, and the new product is not worth making.
3. Otherwise, add the column $B^{-1}A_{n+1}$ to the final tableau and continue with primal simplex. The basis is still feasible and the new column enters.

New constraint $a^\top x \le \beta$:

1. If the current $x^$ satisfies it, $x^$ remains optimal (adding a constraint can never improve the objective).
2. Otherwise, add the row to the tableau (expressed in the current basis) with its slack as a new basic variable. The slack is negative, so the tableau is dual feasible but primal infeasible, and dual simplex restores feasibility.

Change	Check	Re-optimize with
New variable	Reduced cost $c_{n+1} - y^{*\top} A_{n+1}$	Primal simplex
New constraint	Does $x^*$ satisfy it?	Dual simplex
Changed $c_j$ outside range	Reduced costs	Primal simplex
Changed $b_i$ outside range	$B^{-1}b \ge 0$?	Dual simplex

Example: a new pastry and a new labor rule. Pastry earns 4 per batch and uses 1 hour of cake oven and 2 hours of mixer. Separately, a labor agreement limits total batches: $x_1 + x_2 \le 7$.

import numpy as np
from scipy.optimize import linprog

y = np.array([0, 1.5, 1])                     # current shadow prices (bread, cake oven, mixer)
x_star = np.array([2, 6])

# 1) New variable: pastry, profit 4, column (0, 1, 2)
col, profit = np.array([0, 1, 2]), 4
cbar = profit - y @ col                       # 4 - (1.5*1 + 1*2) = 0.5 > 0 -> worth adding
print("pastry reduced cost:", cbar)
res = linprog([-3, -5, -4], A_ub=[[1, 0, 0], [0, 2, 1], [3, 2, 2]], b_ub=[4, 12, 18])
print("new plan (bread, cake, pastry):", res.x, "profit", -res.fun)   # [0. 3. 6.] profit 39.0

# 2) New constraint: x1 + x2 <= 7
a, beta = np.array([1, 1]), 7
print("current plan uses", a @ x_star, "> 7 -> violated, re-optimize")
res = linprog([-3, -5], A_ub=[[1, 0], [0, 2], [3, 2], [1, 1]], b_ub=[4, 12, 18, 7])
print("new plan:", res.x, "profit", -res.fun)                       # [1. 6.] profit 33.0
print("new shadow prices:", -res.ineqlin.marginals + 0.0)           # [0. 1. 0. 3.]

The pastry passes the pricing test (reduced cost $+0.5$), and re-optimizing replaces bread entirely: profit rises from 36 to 39. The labor rule cuts off $(2, 6)$, and the new optimum $(1, 6)$ loses 3. The mixer is no longer a bottleneck, while the labor rule now has a shadow price of 3.

Key takeaways

• A new variable is evaluated by its reduced cost. It is worth adding only if it improves at current shadow prices.
• A new constraint matters only if the current optimum violates it, and then the dual simplex method re-optimizes.
• Warm starts from the old basis typically need few pivots, which is far cheaper than solving from scratch.
• Pricing out new columns is the core idea of column generation (Section 2.6).

Practice

1. (Beginner) Would a pastry with profit 3.4 be worth adding? What is the break-even profit?
2. (Intermediate) Add the pastry column $B^{-1}A_{\text{pastry}}$ to the final bakery tableau from Section 2.3 and perform the primal simplex pivots by hand with pivot. Check that you reach profit 39.
3. (Interview) "Branch-and-bound adds a bound like $x_1 \le 1$ to an LP thousands of times. Why is that fast?" Hint: a new constraint keeps the old basis dual feasible, so dual simplex warm-starts from it.

Parametric Programming

Why it matters. Ranging answers "how far can one number move?". Parametric programming answers "what happens over the entire range?". It traces the optimal value and plan as a parameter sweeps continuously: the full value-of-capacity curve for an investment decision, or the full trade-off curve between two objectives (Chapter 7).

Theory. Two standard cases, with a scalar parameter $\theta$:

• Parametric RHS: $b(\theta) = b + \theta d$. The optimal value $z^(\theta)$ (for max) is piecewise linear and concave in $\theta$. Each linear piece corresponds to one optimal basis, its slope is $y^{\top} d$, and slopes decrease at each breakpoint (diminishing returns).
• Parametric objective: $c(\theta) = c + \theta d$. The optimal value is piecewise linear and convex in $\theta$ (for max). Each piece is one optimal vertex, and its slope is $d^\top x^*$.

The algorithm moves from one breakpoint to the next: at the end of the current range, perform one dual simplex pivot (RHS case) or one primal simplex pivot (objective case), compute the new range, and repeat.

  z*(b3): value of the mixer capacity b3 in the bakery

  z*
  42 +                         *-------------------   slope 0   (bread oven full)
     |                    _.-'
  36 +               _.-*  <- b3 = 18 (today)        slope 1   (y3 = 1)
  30 +          _.-*
     |        .'   b3 = 12
     |      .'                                        slope 2.5 (only cake fits)
     |    .'
   0 *--'----------+-----------+-----------------> b3
     0            12          24

import numpy as np
from scipy.optimize import linprog

c = np.array([3, 5])
A = np.array([[1, 0], [0, 2], [3, 2]])
b = np.array([4, 12, 18], dtype=float)

thetas = np.arange(0, 31, 1.0)                       # mixer capacity b3 from 0 to 30
values = []
for t in thetas:
    res = linprog(-c, A_ub=A, b_ub=[b[0], b[1], t])
    values.append(-res.fun)

slopes = np.diff(values)                             # marginal value of each extra hour
print("b3 range      slope (= shadow price of mixer)")
start = 0
for k in range(1, len(slopes) + 1):
    if k == len(slopes) or not np.isclose(slopes[k], slopes[start]):
        print(f"[{thetas[start]:>4g}, {thetas[k]:>4g}]   {slopes[start]:g}")
        start = k
# b3 range      slope (= shadow price of mixer)
# [   0,   12]   2.5
# [  12,   24]   1
# [  24,   30]   0

The slopes fall from 2.5 to 1 to 0, matching the concavity theorem. The middle piece, $[12, 24]$ with slope 1, is exactly the RHS range and shadow price computed earlier. Parametric programming connects all the pieces.

Key takeaways

• Parametric programming traces the optimal solution as data moves continuously along a direction.
• For max: $z^*$ is piecewise linear concave in the RHS, and piecewise linear convex in the objective coefficients.
• Each linear piece corresponds to one optimal basis. Breakpoints are where the basis changes, and slopes are shadow prices (RHS) or variable values (objective).
• Solvers move between breakpoints with single dual or primal simplex pivots, which is much cheaper than a grid of re-solves.

Practice

1. (Beginner) Using the curve, what is the total value of increasing the mixer from 18 to 28 hours?
2. (Intermediate) Repeat the sweep for the objective $c(\theta) = (\theta, 5)$ with $\theta \in [0, 12]$. Plot $z^(\theta)$, find the breakpoints, and verify that the slopes are the values of $x_1^$.
3. (Interview) "Why is the optimal value concave in the right-hand side for a maximization LP?" Hint: if $x_1$ is optimal for $b_1$ and $x_2$ for $b_2$, what can you say about $\tfrac12(x_1 + x_2)$ for $\tfrac12(b_1 + b_2)$?

<a id="26-advanced-lp-algorithms"></a>

2.6 Advanced LP Algorithms

This subchapter looks beyond the textbook simplex method: how fast LP can be solved in theory and in practice, the polynomial-time algorithms (ellipsoid and interior point), and the decomposition techniques that make LPs with millions of variables or constraints tractable by exploiting their structure.

Computational Complexity of Simplex

Why it matters. The simplex method has solved real LPs efficiently for over 75 years, yet in the worst case it takes exponential time. Understanding this gap between theory and practice explains why solvers offer both simplex and interior point methods, and it is a classic example of why worst-case analysis alone can mislead.

Theory.

• Worst case. Klee and Minty (1972) built a slightly squashed $n$-dimensional cube on which Dantzig's rule visits all $2^n$ vertices. Similar constructions defeat essentially every deterministic pivoting rule studied so far (Bland's, greatest improvement on other variants, steepest edge).
• Practice. On real models the number of pivots typically grows roughly linearly with the number of constraints, often between $m$ and $3m$.
• Explaining the gap. Average-case analysis (Borgwardt, 1982) shows polynomial expected running time on random LPs. Smoothed analysis (Spielman and Teng, 2001) shows that adding tiny random noise to any LP makes the expected number of pivots polynomial, so exponential instances are fragile, like a pencil balanced on its tip.
• Polynomial-time LP. LP is in class P. This was first proven with the ellipsoid method (Khachiyan, 1979), then with interior point methods (Karmarkar, 1984). Their running time is polynomial in $n$, $m$ and $L$, the number of bits needed to encode the data.
• Open problems. Is there a pivoting rule with a polynomial number of pivots? Is there a strongly polynomial LP algorithm, whose number of arithmetic operations does not depend on $L$ (Smale's 9th problem)? The polynomial Hirsch conjecture on the diameter of polytopes is also open (the original Hirsch bound was disproved by Santos in 2010).

Algorithm	Worst case	In practice	Typical use
Simplex (any known rule)	Exponential	~ $m$ to $3m$ pivots, fast	Small/medium LPs, re-optimization, MIP
Ellipsoid	Polynomial, $O(n^2 L)$ iterations	Very slow	Theory (complexity proofs)
Interior point	Polynomial, $O(\sqrt{n}, L)$ iterations	20 to 80 iterations, almost independent of size	Large LPs, solved from scratch

Example: counting pivots on Klee-Minty cubes. The $n$-dimensional Klee-Minty LP is

$$ \max \sum_{j=1}^n 2^{n-j} x_j \quad \text{s.t.} \quad 2\sum_{j=1}^{i-1} 2^{i-j} x_j + x_i \le 5^i \ \ (i = 1..n), \quad x \ge 0 . $$

from fractions import Fraction as F

# Reuses pivot() and simplex() from Section 2.3 (Simplex Tableau and Pivoting Rules topics).

def klee_minty(n):
    """z-row tableau of the n-dimensional Klee-Minty cube (all-slack starting basis)."""
    T = [[F(0)] * (2 * n + 1) for _ in range(n + 1)]
    for j in range(n):
        T[0][j] = F(-2 ** (n - 1 - j))               # objective coefficients 2^(n-j)
    for i in range(n):
        for j in range(i):
            T[i + 1][j] = F(2 ** (i - j + 1))        # 2 * 2^(i-j)
        T[i + 1][i] = F(1)
        T[i + 1][n + i] = F(1)                       # slack
        T[i + 1][-1] = F(5 ** (i + 1))
    return T, list(range(n, 2 * n))

print(" n   Dantzig   Bland   greatest-improvement")
for n in range(1, 9):
    counts = []
    for rule in ("dantzig", "bland", "greatest"):
        T, basis = klee_minty(n)
        status, pivots = simplex(T, basis, rule, max_iter=10**4)
        counts.append(pivots)
    print(f"{n:>2} {counts[0]:>9} {counts[1]:>7} {counts[2]:>10}")
#  n   Dantzig   Bland   greatest-improvement
#  1         1       1          1
#  2         3       3          1
#  3         7       5          1
#  4        15       9          1
#  5        31      15          1
#  6        63      25          1
#  7       127      41          1
#  8       255      67          1

Dantzig's rule needs exactly $2^n - 1$ pivots, visiting every vertex. Bland's rule also grows exponentially here (each count is roughly 1.6 times the previous one). Greatest improvement jumps straight to the optimum on this cube, but other "deformed cubes" defeat it too.

Key takeaways

• Simplex is exponential in the worst case (Klee-Minty) but typically needs a number of pivots roughly linear in $m$.
• Smoothed analysis explains the good practical behavior: bad instances disappear under small random perturbations.
• LP is solvable in polynomial time (ellipsoid, interior point), but a polynomial pivoting rule and a strongly polynomial algorithm are still unknown.
• Worst-case complexity and practical performance are different questions. Benchmark on your own models.

Practice

1. (Beginner) Solve the 2D Klee-Minty LP graphically and trace the 3 Dantzig pivots on the drawing.
2. (Intermediate) Time linprog(method="highs-ds") (dual simplex) and linprog(method="highs-ipm") (interior point) on Klee-Minty cubes with $n$ up to 50. Do modern solvers suffer from the exponential behavior? Why not?
3. (Interview) "If the simplex method is exponential in the worst case, why is it still the default in many solvers?" Hint: consider warm starts, basic (vertex) solutions, and the needs of branch-and-bound.

Ellipsoid Method

Why it matters. The ellipsoid method was the first algorithm to prove that LP can be solved in polynomial time (Khachiyan, 1979), a result that made front-page news at the time. It is too slow for practical use, but it has a powerful theoretical consequence: you can optimize over any convex set for which you can efficiently answer "is this point feasible, and if not, which constraint is violated?". This is used to prove that many combinatorial problems are polynomially solvable.

Theory. The method shrinks an ellipsoid that is guaranteed to contain the solution, like a "hot and cold" search that halves the search region at each step:

1. Start with a large ellipsoid $E_0 = {x : (x - x_0)^\top P_0^{-1}(x - x_0) \le 1}$ that contains the feasible region (or its relevant part).
2. At the center $x_k$:
   • If some constraint $a_i^\top x \le b_i$ is violated, the feasible region lies in the half-ellipsoid ${x \in E_k : a_i^\top x \le a_i^\top x_k}$ (feasibility cut).
   • If $x_k$ is feasible, record it, and keep only points at least as good: ${x \in E_k : c^\top x \ge c^\top x_k}$ (objective cut).
3. Replace the half-ellipsoid by the smallest ellipsoid containing it. With the cut normal $a$ and $g = P_k a / \sqrt{a^\top P_k a}$:

$$ x_{k+1} = x_k - \frac{1}{n+1}, g, \qquad P_{k+1} = \frac{n^2}{n^2 - 1}\left(P_k - \frac{2}{n+1}, g g^\top\right). $$

4. The volume shrinks by a factor $< e^{-1/(2(n+1))}$ per iteration, so after $O(n^2 L)$ iterations the ellipsoid is too small to contain any vertex not already found.

   E_k (contains the solution)        cut through the center        E_{k+1}: smallest ellipsoid
                                                                   around the kept half
        .-------.                         .---|---.                       .----.
      /           \                     /     |###  \                    /  ###  \
     |      x_k    |       ===>        |      x_k####|      ===>        |   ####  |
      \           /                     \     |###  /                    \  ###  /
        '-------'                         '---|---'                       '----'
                                    keep only the side that can
                                    contain the solution

The method needs only a separation oracle: a procedure that either confirms $x$ is feasible or returns one violated inequality. The constraints never have to be listed explicitly, which is why the method works even for LPs with exponentially many constraints (Grötschel, Lovász and Schrijver, 1981).

import numpy as np

# Bakery as A x <= b (with nonnegativity rows), maximize c^T x
A = np.array([[1, 0], [0, 2], [3, 2], [-1, 0], [0, -1]], dtype=float)
b = np.array([4, 12, 18, 0, 0], dtype=float)
c = np.array([3, 5], dtype=float)

n = 2
x = np.zeros(n)                     # center of the initial ellipsoid
P = np.eye(n) * 100.0**2            # ball of radius 100: surely contains the region
best = None
for k in range(1, 201):
    violation = A @ x - b
    i = int(np.argmax(violation))
    if violation[i] > 0:            # separation oracle: a violated constraint
        a = A[i]                    # feasibility cut  a^T x <= a^T x_k
    else:
        if best is None or c @ x > c @ best:
            best = x.copy()         # feasible center: remember it
        a = -c                      # objective cut  c^T x >= c^T x_k
    g = P @ a / np.sqrt(a @ P @ a)
    x = x - g / (n + 1)
    P = n**2 / (n**2 - 1) * (P - 2 / (n + 1) * np.outer(g, g))
    if k in (10, 25, 50, 100, 200):
        size = np.sqrt(np.linalg.det(P))          # proportional to the area of the ellipsoid
        print(f"iter {k:>3}: best x = {np.round(best, 4)}, value {c @ best:.4f}, area ~ {size:.1e}")
# iter  10: best x = [3.0566 1.7818], value 18.0790, area ~ 7.3e+02
# iter  25: best x = [2.4232 4.696 ], value 30.7498, area ~ 1.4e+01
# iter  50: best x = [2.0085 5.9485], value 35.7678, area ~ 2.1e-02
# iter 100: best x = [2. 6.], value 35.9999, area ~ 4.3e-08
# iter 200: best x = [2. 6.], value 36.0000, area ~ 1.9e-19

Convergence is steady but slow: each iteration cuts the volume by a constant factor, regardless of how "easy" the problem is. In $n$ dimensions, this factor gets closer to 1, so the method needs many more iterations than simplex or interior point methods.

Key takeaways

• The ellipsoid method keeps an ellipsoid containing the solution, cuts it through the center, and re-encloses the kept half.
• It runs in polynomial time and proved that LP is in P, but it is far too slow and numerically delicate for practice.
• It only needs a separation oracle, so it works for LPs with exponentially many (implicit) constraints.
• Its main legacy is theoretical: separation and optimization are polynomially equivalent.

Practice

1. (Beginner) Why does a feasible center lead to an objective cut instead of stopping immediately?
2. (Intermediate) Record the volume ratio $\sqrt{\det P_{k+1} / \det P_k}$ in the code and compare it with the bound $e^{-1/(2(n+1))}$ for $n = 2$.
3. (Interview) "The ellipsoid method is polynomial and the simplex method is not. Why does nobody use the ellipsoid method?" Hint: compare the constants and the dependence on $L$, and think about numerical precision with nearly flat ellipsoids.

Interior Point Methods

Why it matters. Interior point methods (IPMs) are the polynomial-time algorithms that also work in practice. On large, sparse LPs they are often much faster than simplex, and their number of iterations barely grows with problem size. Every major solver (Gurobi, CPLEX, HiGHS, MOSEK) includes one, and the same machinery extends to convex quadratic and conic optimization (Chapter 5).

Theory. Simplex walks along the boundary from vertex to vertex. IPMs cut through the inside of the feasible region and approach the optimal vertex from within.

The key tool is the logarithmic barrier. For $\max{c^\top x : Ax \le b}$, define

$$ \max_x \ \ c^\top x + \mu \sum_{i=1}^m \ln\big(b_i - a_i^\top x\big), \qquad \mu > 0 . $$

• The log term goes to $-\infty$ at the boundary, so it acts as a force field that keeps $x$ strictly inside the region.
• For each $\mu$, the barrier problem is smooth and strictly concave, with a unique maximizer $x(\mu)$. The curve ${x(\mu) : \mu > 0}$ is the central path.
• As $\mu \to 0$, $x(\mu) \to x^$, and the suboptimality is bounded: $c^\top x^ - c^\top x(\mu) \le m\mu$.

Algorithm (barrier method): start inside, take a few Newton steps to approximately maximize the barrier problem for the current $\mu$, reduce $\mu$ (e.g. by a factor of 10), and repeat. Short-step variants provably need $O(\sqrt{m}, \log(1/\varepsilon))$ Newton steps, and in practice 20 to 80 steps suffice for almost any LP.

  Central path inside the bakery polygon (approaching the optimum from within)

  x2
  6 *-----------* (2,6)  <- x(mu) -> x* as mu -> 0
    |         .'|\
    |       .'  | \           mu = 0.1 : (1.99, 5.97)
    |     o     |  \          mu = 1   : (1.89, 5.68)
    |   .'      |   * (4,3)   mu = 10  : (1.59, 3.68)
    |  o        |   |
    | (analytic center, mu -> infinity)
    *-----------+---*----> x1

import numpy as np

A = np.array([[1, 0], [0, 2], [3, 2], [-1, 0], [0, -1]], dtype=float)   # includes x >= 0
b = np.array([4, 12, 18, 0, 0], dtype=float)
c = np.array([3, 5], dtype=float)

x = np.array([1.0, 1.0])                         # strictly feasible start
for mu in (10, 1, 0.1, 0.01, 0.001):
    for _ in range(50):                          # Newton's method on the barrier problem
        s = b - A @ x                            # slacks, must stay > 0
        grad = c - mu * A.T @ (1 / s)
        hess = -mu * A.T @ np.diag(1 / s**2) @ A
        dx = -np.linalg.solve(hess, grad)        # Newton direction
        t = 1.0
        while np.any(b - A @ (x + t * dx) <= 0): # backtrack to stay strictly inside
            t *= 0.5
        x = x + t * dx
        if np.linalg.norm(dx) < 1e-10:
            break
    print(f"mu = {mu:<6} x(mu) = {np.round(x, 4)}  value {c @ x:.4f}  gap bound m*mu = {5 * mu:g}")
# mu = 10     x(mu) = [1.5939 3.6782]  value 23.1725  gap bound m*mu = 50
# mu = 1      x(mu) = [1.8854 5.6813]  value 34.0627  gap bound m*mu = 5
# mu = 0.1    x(mu) = [1.9888 5.9668]  value 35.8006  gap bound m*mu = 0.5
# mu = 0.01   x(mu) = [1.9989 5.9967]  value 35.9800  gap bound m*mu = 0.05
# mu = 0.001  x(mu) = [1.9999 5.9997]  value 35.9980  gap bound m*mu = 0.005

IPMs never land exactly on a vertex. Solvers finish with a crossover step that moves from the IPM solution to a nearby optimal basis, which is needed for sensitivity analysis and warm starts.

	Simplex	Interior point
Path	Along edges, vertex to vertex	Through the interior, along the central path
Iterations	Many cheap pivots	Few expensive Newton steps (one sparse linear system each)
Large sparse LPs	Can be slow	Usually faster
Warm start after small changes	Excellent	Poor
Output	Exact vertex + basis	Near-optimal interior point (crossover needed for a basis)

Key takeaways

• IPMs follow the central path of a log-barrier problem through the interior toward the optimum.
• Each iteration is a Newton step (a linear system). Typically 20 to 80 iterations are needed, regardless of size.
• The barrier parameter $\mu$ bounds the suboptimality ($\le m\mu$), and driving $\mu \to 0$ gives the optimum.
• IPMs are best for large LPs solved from scratch. Simplex is best for re-optimization and when a basis is needed.

Practice

1. (Beginner) Explain why $\ln(b_i - a_i^\top x)$ keeps iterates strictly feasible.
2. (Intermediate) In the code, count the total number of Newton steps for each value of $\mu$. Then try reducing $\mu$ by a factor of 100 per stage instead of 10 and compare the counts.
3. (Interview) "When would you choose an interior point method over simplex for an LP?" Hint: consider size, sparsity, whether you need a basis, and whether you will re-solve after small changes.

Primal-Dual Path-Following Algorithms

Why it matters. The interior point method in modern solvers is not the pure barrier method of the previous topic. It is a primal-dual method that updates $x$, $y$ and the dual slacks together, with Mehrotra's predictor-corrector heuristic. It is the most robust and efficient IPM, and understanding it explains solver logs ("primal infeasibility, dual infeasibility, complementarity") line by line.

Theory. For standard form $\min{c^\top x : Ax = b, x \ge 0}$ with dual $\max{b^\top y : A^\top y + s = c, s \ge 0}$, the optimality conditions are

$$ Ax = b, \qquad A^\top y + s = c, \qquad x_j s_j = 0 \ \ \forall j, \qquad x, s \ge 0 . $$

The last equation is complementary slackness, and it is the only nonlinear part. Path-following methods relax it to $x_j s_j = \mu$ for all $j$. For each $\mu > 0$ this defines the (primal-dual) central path, and as $\mu \to 0$ it converges to an optimal primal-dual pair. Each iteration:

1. Measure the duality measure $\mu = x^\top s / n$ (the average complementarity).
2. Compute a Newton step $(\Delta x, \Delta y, \Delta s)$ toward the point with target $\sigma\mu$, where the centering parameter $\sigma \in (0, 1)$:

$$ \begin{bmatrix} A & 0 & 0 \ 0 & A^\top & I \ S & 0 & X \end{bmatrix} \begin{bmatrix} \Delta x \ \Delta y \ \Delta s \end{bmatrix}

\begin{bmatrix} b - Ax \ c - A^\top y - s \ \sigma\mu e - XSe \end{bmatrix} . $$

3. Take the largest step $\alpha \le 1$ that keeps $x > 0$ and $s > 0$ (scaled by 0.99 to stay strictly inside).

"Infeasible" variants, used in practice, start from any $x, s > 0$ and reduce the primal and dual residuals along with $\mu$. Mehrotra's predictor-corrector first computes an affine-scaling step ($\sigma = 0$) to estimate how much $\mu$ can drop, then chooses $\sigma$ adaptively and adds a second-order correction. Both steps reuse one factorization.

import numpy as np

# Bakery in standard form: min -3x1 - 5x2  s.t.  A x = b, x >= 0  (x includes slacks)
A = np.array([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=float)
b = np.array([4, 12, 18], dtype=float)
c = np.array([-3, -5, 0, 0, 0], dtype=float)
m, n = A.shape

x, s, y = np.ones(n), np.ones(n), np.zeros(m)   # infeasible start, only x, s > 0 needed
sigma = 0.1                                     # centering parameter

def max_step(v, dv):
    """Largest alpha in (0, 1] with v + alpha * dv >= 0."""
    neg = dv < 0
    return min(1.0, np.min(-v[neg] / dv[neg])) if neg.any() else 1.0

for k in range(1, 50):
    r_p, r_d, mu = b - A @ x, c - A.T @ y - s, x @ s / n
    print(f"{k:>2}  mu = {mu:8.1e}  primal inf = {np.linalg.norm(r_p):7.1e}  "
          f"dual inf = {np.linalg.norm(r_d):7.1e}  obj = {c @ x:9.4f}")
    if max(mu, np.linalg.norm(r_p), np.linalg.norm(r_d)) < 1e-8:
        break
    K = np.block([[A, np.zeros((m, m)), np.zeros((m, n))],          # Newton (KKT) matrix
                  [np.zeros((n, n)), A.T, np.eye(n)],
                  [np.diag(s), np.zeros((n, m)), np.diag(x)]])
    d = np.linalg.solve(K, np.r_[r_p, r_d, sigma * mu - x * s])
    dx, dy, ds = d[:n], d[n:n + m], d[n + m:]
    alpha = min(1.0, 0.99 * min(max_step(x, dx), max_step(s, ds)))
    x, y, s = x + alpha * dx, y + alpha * dy, s + alpha * ds

print("x =", np.round(x, 6), " y =", np.round(y, 6) + 0.0)

 1  mu =  1.0e+00  primal inf = 1.5e+01  dual inf = 7.4e+00  obj =   -8.0000
 2  mu =  6.3e-01  primal inf = 1.2e+01  dual inf = 6.0e+00  obj =  -12.8735
 3  mu =  1.6e-01  primal inf = 5.9e+00  dual inf = 2.9e+00  obj =  -25.5614
 4  mu =  3.2e-02  primal inf = 1.4e+00  dual inf = 6.6e-01  obj =  -33.5987
 5  mu =  3.1e-03  primal inf = 1.4e-02  dual inf = 6.6e-03  obj =  -35.9695
 6  mu =  3.4e-04  primal inf = 1.4e-04  dual inf = 6.6e-05  obj =  -35.9991
 ...
11  mu =  5.2e-09  primal inf = 1.4e-14  dual inf = 7.2e-15  obj =  -36.0000
x = [2. 6. 2. 0. 0.]  y = [ 0.  -1.5 -1. ]

The log has exactly the columns of a commercial solver's barrier log. All three residuals shrink by about a factor of 10 per iteration, and the dual solution $y$ (here in minimization signs) comes out together with $x$.

Key takeaways

• Primal-dual methods apply Newton's method to the optimality conditions with complementarity relaxed to $x_j s_j = \sigma\mu$.
• Each iteration reduces primal infeasibility, dual infeasibility and the duality measure $\mu = x^\top s / n$ together.
• Step lengths keep $x$ and $s$ strictly positive. Mehrotra's predictor-corrector chooses the centering adaptively.
• The main cost per iteration is solving one sparse linear system (normal equations $AXS^{-1}A^\top$). Exploiting sparsity is what makes IPMs scale.

Practice

1. (Beginner) Identify which block row of the Newton system corresponds to primal feasibility, dual feasibility and complementarity.
2. (Intermediate) Eliminate $\Delta s$ and $\Delta x$ to derive the normal equations $A (X S^{-1}) A^\top \Delta y = \dots$, then implement them in place of the full KKT solve.
3. (Interview) "A barrier solve stalls at relative gap $10^{-6}$ and reports numerical trouble. What is happening as $\mu \to 0$?" Hint: look at the ratios $x_j / s_j$. Some go to 0 and some to infinity, so $AXS^{-1}A^\top$ becomes very ill-conditioned.

Column Generation

Why it matters. Some LPs have far too many variables to write down: every possible cutting pattern for paper rolls, every possible crew schedule for an airline, every feasible route for a delivery vehicle. Column generation solves them by working with a small subset of columns and generating new ones only when they would improve the solution. It powers airline crew scheduling, vehicle routing and cutting stock software.

Theory. The idea is pricing (Section 2.5) applied systematically:

1. Restricted master problem (RMP): solve the LP using only a manageable subset of columns.
2. Pricing subproblem: with the RMP's dual values $y$, find the column with the most attractive reduced cost, $\min_j (c_j - y^\top A_j)$ for a min problem. Columns are not enumerated. Instead, the pricing problem is a structured optimization problem (knapsack, shortest path, ...) whose solutions are columns.
3. If the best reduced cost is negative (min problem), add that column to the RMP and repeat. Otherwise, no column can improve the RMP, and the RMP solution is optimal for the full LP.

Cutting stock problem (Gilmore and Gomory, 1961). Rolls of width $W = 100$ must be cut into pieces of widths $w_i$ with demands $d_i$. A pattern $p$ says how many pieces of each width are cut from one roll ($a_{ip}$). With $x_p$ = number of rolls cut with pattern $p$:

$$ \min \sum_p x_p \quad \text{s.t.} \quad \sum_p a_{ip} x_p \ge d_i \ \ \forall i, \qquad x \ge 0 . $$

A pattern's reduced cost is $1 - \sum_i y_i a_{ip}$, so the most attractive pattern solves a knapsack problem: $\max \sum_i y_i a_i$ s.t. $\sum_i w_i a_i \le W$, $a_i$ nonnegative integers.

  Column generation loop

  +--------------------------+   duals y     +------------------------------+
  | restricted master LP     | ------------> | pricing: knapsack            |
  | (few patterns)           |               | max sum y_i a_i, sum w_i a_i <= W
  +--------------------------+ <------------ +------------------------------+
            ^                  new pattern if 1 - value < 0
            |
            +---- stop when no pattern has negative reduced cost

import numpy as np
from scipy.optimize import linprog

W = 100
widths = np.array([45, 36, 31, 14])
demand = np.array([97, 610, 395, 211])

def knapsack(values, weights, cap):
    """Unbounded integer knapsack by dynamic programming. Returns (value, pattern)."""
    best = np.zeros(cap + 1)
    take = [-1] * (cap + 1)
    for c in range(1, cap + 1):
        best[c], take[c] = best[c - 1], -1                # option: waste one unit of width
        for i, w in enumerate(weights):
            if w <= c and best[c - w] + values[i] > best[c] + 1e-12:
                best[c], take[c] = best[c - w] + values[i], i
    pattern, c = np.zeros(len(weights), dtype=int), cap
    while c > 0:                                          # rebuild the chosen pattern
        if take[c] == -1:
            c -= 1
        else:
            pattern[take[c]] += 1
            c -= weights[take[c]]
    return best[cap], pattern

# Initial columns: homogeneous patterns (as many pieces of one width as fit)
patterns = [np.eye(len(widths), dtype=int)[i] * (W // widths[i]) for i in range(len(widths))]

for it in range(1, 50):
    P = np.array(patterns).T                              # columns = patterns
    rmp = linprog(np.ones(P.shape[1]), A_ub=-P, b_ub=-demand)   # P x >= demand
    y = -rmp.ineqlin.marginals                            # duals of the demand rows
    value, new = knapsack(y, widths, W)
    print(f"iter {it}: rolls = {rmp.fun:.2f}, best new pattern {new}, "
          f"reduced cost {1 - value:+.4f}")
    if 1 - value >= -1e-9:                                # no improving column: optimal
        break
    patterns.append(new)

for p, x in zip(patterns, rmp.x):
    if x > 1e-9:
        print(f"cut {x:7.2f} rolls with pattern {p}")
# iter 1: rolls = 515.31, best new pattern [0 2 0 2], reduced cost -0.2857
# iter 2: rolls = 485.17, best new pattern [0 1 2 0], reduced cost -0.1667
# iter 3: rolls = 452.25, best new pattern [0 2 0 0], reduced cost +0.0000
# cut   48.50 rolls with pattern [2 0 0 0]
# cut  206.25 rolls with pattern [0 2 0 2]
# cut  197.50 rolls with pattern [0 1 2 0]

Out of hundreds of possible patterns, only 6 were ever needed. The LP bound of 452.25 proves that any integer cutting plan needs at least 453 rolls. Rounding the LP solution gives a good integer plan, and branch-and-price (Chapter 4) closes the remaining gap exactly.

Key takeaways

• Column generation solves LPs with huge numbers of columns by adding only columns with improving reduced costs.
• The restricted master gives duals, the pricing subproblem finds the best new column, and the process stops when no column prices out.
• Pricing is a structured optimization problem (knapsack, shortest path) rather than an enumeration.
• Convergence often has a long "tail" with small improvements, so stabilization techniques and early stopping with Lagrangian bounds are used in practice.

Practice

1. (Beginner) Verify that the pattern $[0, 2, 0, 2]$ fits in a roll of width 100 and compute its reduced cost with the iteration-1 duals $y = (1/2, 1/2, 1/3, 1/7)$.
2. (Intermediate) Round the final LP solution up to integers and count the rolls. How far is it from the LP bound? Can you do better by adjusting a few patterns by hand?
3. (Interview) "How would you model airline crew scheduling with column generation? What are the columns and what is the pricing problem?" Hint: a column is a legal sequence of flights for one crew (a pairing), and pricing is a resource-constrained shortest path in a flight network.

Dantzig-Wolfe Decomposition

Why it matters. Many large LPs consist of independent pieces (plants, regions, time periods, commodities) tied together by a few shared constraints, such as a common budget or raw material. Dantzig-Wolfe decomposition solves each piece separately and coordinates them through prices. It mirrors how a company coordinates divisions with internal transfer prices, and it is the theoretical foundation of branch-and-price.

Theory. Consider a block-angular LP with $K$ blocks:

                 x_1     x_2    ...    x_K
             +-------+-------+-----+-------+
  linking    |  L_1  |  L_2  | ... |  L_K  |  <= b_0     (few shared constraints)
             +-------+-------+-----+-------+
  block 1    |  D_1  |       |     |       |  <= e_1
  block 2    |       |  D_2  |     |       |  <= e_2
   ...       |       |       | ... |       |
  block K    |       |       |     |  D_K  |  <= e_K
             +-------+-------+-----+-------+

$$ \max \sum_k c_k^\top x_k \quad \text{s.t.} \quad \sum_k L_k x_k \le b_0, \qquad x_k \in X_k = {x_k \ge 0 : D_k x_k \le e_k}. $$

By the Minkowski-Weyl theorem, every point of a bounded polytope $X_k$ is a convex combination of its vertices $x_k^p$: $x_k = \sum_p \lambda_{kp} x_k^p$ with $\sum_p \lambda_{kp} = 1$ and $\lambda \ge 0$. Substituting gives the master problem in the variables $\lambda$:

$$ \max \sum_{k,p} (c_k^\top x_k^p), \lambda_{kp} \quad \text{s.t.} \quad \sum_{k,p} (L_k x_k^p), \lambda_{kp} \le b_0 \ \ (\text{duals } \pi), \qquad \sum_p \lambda_{kp} = 1 \ \ (\text{duals } \mu_k), \qquad \lambda \ge 0 . $$

It has few rows (linking constraints + one convexity row per block) but a huge number of columns (all vertices), so it is solved by column generation. The pricing problem decomposes by block:

$$ \text{for each } k: \quad \max_{x_k \in X_k}\ (c_k - L_k^\top \pi)^\top x_k - \mu_k . $$

If the value is positive, the maximizer is a new column for block $k$. The economic reading: the master sets prices $\pi$ for the shared resources, each block (division) optimizes its own plan with those resources charged at $\pi$, and proposes the plan to the master.

Example: two plants sharing a raw material. Plant 1 makes products with profits $(5, 4)$ under local limits $2x_1 + x_2 \le 8$, $x_1 + 3x_2 \le 9$. Plant 2 makes products with profits $(4, 7)$ under $x_3 + x_4 \le 5$, $x_3 + 2x_4 \le 8$. Both draw on a raw material with 12 units available: $x_1 + 2x_2 + 2x_3 + 3x_4 \le 12$.

import numpy as np
from scipy.optimize import linprog

c = [np.array([5., 4.]), np.array([4., 7.])]                        # profits per block
D = [np.array([[2., 1.], [1., 3.]]), np.array([[1., 1.], [1., 2.]])] # local constraints
e = [np.array([8., 9.]), np.array([5., 8.])]
L = [np.array([[1., 2.]]), np.array([[2., 3.]])]                     # raw material use
b0 = np.array([12.])
K = 2

columns = [[np.zeros(2)] for _ in range(K)]      # start: each plant proposes "produce nothing"
for it in range(1, 20):
    # Restricted master over the proposals collected so far
    profit, usage, owner = [], [], []
    for k in range(K):
        for xp in columns[k]:
            profit.append(c[k] @ xp); usage.append(L[k] @ xp); owner.append(k)
    convexity = np.array([[1.0 if o == k else 0.0 for o in owner] for k in range(K)])
    master = linprog(-np.array(profit), A_ub=np.array(usage).T, b_ub=b0,
                     A_eq=convexity, b_eq=np.ones(K))
    pi = -master.ineqlin.marginals                # price of the raw material
    mu = -master.eqlin.marginals                  # value of each plant's current mix
    print(f"iter {it}: master profit {-master.fun + 0.0:.3f}, raw material price {pi[0]:.3f}")

    improved = False
    for k in range(K):                            # each plant optimizes at price pi
        sub = linprog(-(c[k] - pi @ L[k]), A_ub=D[k], b_ub=e[k])
        reduced = -sub.fun - mu[k]
        if reduced > 1e-9:
            columns[k].append(sub.x)
            improved = True
            print(f"   plant {k + 1} proposes {sub.x}, reduced cost {reduced:.3f}")
    if not improved:
        break

# Recover x_k as the convex combination of the proposals
weights, j = master.x, 0
for k in range(K):
    xk = sum(weights[j + p] * xp for p, xp in enumerate(columns[k]))
    j += len(columns[k])
    print(f"plant {k + 1} plan: {np.round(xk, 3)}")

full = linprog(-np.r_[c[0], c[1]], A_ub=np.block([[D[0], np.zeros((2, 2))],
               [np.zeros((2, 2)), D[1]], [L[0], L[1]]]), b_ub=np.r_[e[0], e[1], b0])
print("check, full LP profit:", round(-full.fun, 3))

iter 1: master profit 0.000, raw material price 0.000
   plant 1 proposes [3. 2.], reduced cost 23.000
   plant 2 proposes [2. 3.], reduced cost 29.000
iter 2: master profit 34.154, raw material price 2.231
   plant 1 proposes [4. 0.], reduced cost 3.692
   plant 2 proposes [0. 4.], reduced cost 1.231
iter 3: master profit 38.667, raw material price 2.333
plant 1 plan: [4. 0.]
plant 2 plan: [0.    2.667]
check, full LP profit: 38.667

The master never sees the plants' local constraints, and the plants never see each other. Prices on the shared raw material are enough to coordinate them to the global optimum.

Key takeaways

• Dantzig-Wolfe applies to block-angular LPs: independent blocks linked by a few coupling constraints.
• Each block is represented by convex combinations of its vertices. The master LP coordinates them, and column generation adds vertices as needed.
• Pricing splits into independent block subproblems that use the master's dual prices on the linking constraints, so they can be solved in parallel.
• The master's LP bound is often much stronger for integer problems, which is why Dantzig-Wolfe underlies branch-and-price.

Practice

1. (Beginner) In iteration 2, the price of raw material is 2.231. Compute plant 1's "adjusted profits" $c_1 - L_1^\top \pi$ and explain its new proposal.
2. (Intermediate) Add a third plant with profits $(6, 3)$, local constraint $x_5 + x_6 \le 4$ and raw material use $(1, 1)$. Update the code (it is written for general $K$) and compare with the full LP.
3. (Interview) "How does Dantzig-Wolfe relate to how a large company with independent divisions could set internal prices?" Hint: the master is headquarters setting transfer prices for shared resources, and each division maximizes its own profit net of those prices.

Benders Decomposition

Why it matters. Benders decomposition handles the opposite structure to Dantzig-Wolfe: a few complicating variables (which facilities to open, how much capacity to build, a first-stage decision under uncertainty) that, once fixed, leave an easy LP. It is the standard method for two-stage stochastic programs (Chapter 9) and for many network design and facility location models.

Theory.

  Dantzig-Wolfe: linking CONSTRAINTS            Benders: linking VARIABLES

  [ L1  L2  L3 ]  <- shared rows                  y   x1  x2  x3
  [ D1         ]                                [ F1  D1          ]
  [     D2     ]                                [ F2      D2      ]
  [         D3 ]                                [ F3          D3  ]
                                                  ^ shared columns

Consider $\min_y\ \big(f^\top y + Q(y)\big)$ over $y \in Y$, where the recourse function is

$$ Q(y) = \min_x \ { q^\top x : Dx \ge h - Fy,\ x \ge 0 }. $$

• $Q(y)$ is an LP's optimal value as a function of its right-hand side, so it is piecewise linear and convex in $y$ (Section 2.5, parametric programming).
• By duality, $Q(y) = \max_u { u^\top (h - Fy) : D^\top u \le q,\ u \ge 0}$. The dual feasible set does not depend on $y$, so solving the subproblem at $\bar y$ with dual solution $\bar u$ gives a valid linear lower bound for every $y$:

$$ Q(y) \ \ge\ \bar u^\top (h - Fy) \qquad (\text{optimality cut}). $$

Algorithm:

1. Master: $\min f^\top y + \theta$ over $y \in Y$, subject to all cuts found so far ($\theta \ge \bar u_t^\top (h - Fy)$). Its value is a lower bound on the true optimum.
2. Subproblem: fix $y = \bar y$ from the master and solve the LP for $Q(\bar y)$. Then $f^\top \bar y + Q(\bar y)$ is an upper bound.
3. Add the cut from the subproblem's dual solution and repeat until the bounds meet.

If the subproblem is infeasible for some $\bar y$, a dual ray gives a feasibility cut that removes $\bar y$. With "complete recourse" (the subproblem is always feasible) this case cannot occur.

Example: how much capacity to build. Building capacity $y \in [0, 20]$ costs 4 per unit. Afterwards, two customers with demands 6 and 8 are served from that capacity at 1 per unit, or outsourced at 10 (customer 1) or 4.5 (customer 2) per unit.

import numpy as np
from scipy.optimize import linprog

f, demand, outsource = 4.0, np.array([6., 8.]), np.array([10., 4.5])

def subproblem(y):
    """Q(y) = min x1 + x2 + 10 u1 + 4.5 u2  s.t.  x1 + x2 <= y,  x_i + u_i >= d_i,  all >= 0."""
    res = linprog(np.r_[1, 1, outsource],
                  A_ub=[[1, 1, 0, 0], [-1, 0, -1, 0], [0, -1, 0, -1]],
                  b_ub=np.r_[y, -demand])
    slope = res.ineqlin.marginals[0]       # dQ/dy: dual value of the capacity row
    return res.fun, slope

cuts, y = [], 0.0                          # start with no capacity
for it in range(1, 20):
    Q, slope = subproblem(y)
    upper = f * y + Q                      # cost of a real, implementable plan
    cuts.append((Q, slope, y))             # cut: theta >= Q + slope * (y' - y)
    # Master: min f y + theta  s.t. all cuts, 0 <= y <= 20, theta >= 0
    A_ub = [[s, -1] for (_, s, _) in cuts]
    b_ub = [s * yb - q for (q, s, yb) in cuts]
    master = linprog([f, 1], A_ub=A_ub, b_ub=b_ub, bounds=[(0, 20), (0, None)])
    lower = master.fun
    print(f"iter {it}: y = {y:6.3f}  Q(y) = {Q:6.2f}  cut slope {slope:5.1f}  "
          f"bounds [{lower:.2f}, {upper:.2f}]")
    if upper - lower < 1e-6:
        break
    y = master.x[0]
# iter 1: y =  0.000  Q(y) =  96.00  cut slope  -9.0  bounds [42.67, 96.00]
# iter 2: y = 10.667  Q(y) =  25.67  cut slope  -3.5  bounds [66.00, 68.33]
# iter 3: y =  6.000  Q(y) =  42.00  cut slope  -3.5  bounds [66.00, 66.00]

The optimum builds $y = 6$: enough capacity for customer 1 only. Each unit for customer 1 saves $10 - 1 = 9$ at a building cost of 4. Customer 2 saves only $4.5 - 1 = 3.5 < 4$ per unit, so outsourcing is cheaper. The cut slopes $-9$ and $-3.5$ are exactly these marginal savings, which the master learns one cut at a time.

Key takeaways

• Benders targets complicating variables: fix them, solve an easy LP, and send back a cut built from its dual.
• The recourse function $Q(y)$ is convex and piecewise linear, and optimality cuts are its supporting hyperplanes.
• The master gives lower bounds and the subproblems give upper bounds, so the gap measures convergence.
• When the master has integer variables (facility opening, network design), Benders is the standard decomposition. With many scenarios, the subproblems split per scenario (L-shaped method, Chapter 9).

Practice

1. (Beginner) Compute $Q(0)$, $Q(6)$ and $Q(14)$ by hand and sketch $4y + Q(y)$ on $[0, 20]$.
2. (Intermediate) Make the outsourcing cost of customer 2 equal to 6 and re-run. Predict the new optimal capacity before running the code.
3. (Interview) "Compare Benders and Dantzig-Wolfe decomposition. When would you choose each?" Hint: look at the matrix structure. Shared rows point to Dantzig-Wolfe, shared columns to Benders, and the two are dual to each other.

<a id="3-network-optimization"></a>