Prev
Unit 11 · Numerical Optimisation
Constrained Optimisation and Linear Programming
65m100 XP
Next
⛰️Numerical Optimisation

Constrained Optimisation and Linear Programming

Lagrange multipliers and KKT conditions, penalty and SQP methods, linear programming, minimum-mass truss layout

The Pentagon, 1947
George Dantzig is a mathematical adviser to the US Air Force Comptroller. His task, under a project code-named SCOOP, is to mechanise planning: how to allocate aircraft, crews, fuel and training across thousands of activities, subject to hundreds of limits, to meet the Air Force's objectives at least cost. Until now such plans have been assembled by hand, by experience and negotiation. Dantzig formulates the problem with every relationship linear — so much fuel per sortie, so many crews per aircraft — and invents an algorithm to solve it: the simplex method, which walks from corner to corner of the region of feasible plans, improving the objective at each step, until no neighbouring corner is better. Unknown to him, Leonid Kantorovich in Leningrad had formulated similar problems in 1939 for a plywood trust, but his work was little known outside the Soviet Union. Linear programming spread from military logistics to oil refineries, airlines, steel mills and electricity grids. Kantorovich shared the 1975 Nobel Prize in economics for it with Tjalling Koopmans. Today's solvers handle millions of variables, and an engineer meets them in structural design, scheduling and power dispatch.
A constrained problem minimises f(x)f(\mathbf x) subject to equality constraints hi(x)=0h_i(\mathbf x) = 0 and inequality constraints gj(x)≥0g_j(\mathbf x) \ge 0. The solution often lies on the boundary of the feasible set, where some inequality constraints are active (hold with equality). Setting the gradient to zero no longer works: the objective may still be decreasing in a direction that leaves the feasible set.
💡Lagrange multipliers and the KKT conditions
At a constrained minimum, ∇f\nabla f must be a combination of the gradients of the active constraints — otherwise some feasible direction would still decrease ff: ∇f(x∗)=∑iλi∇hi(x∗)+∑jμj∇gj(x∗),μj≥0,μjgj(x∗)=0.\nabla f(\mathbf x^*) = \sum_i\lambda_i\nabla h_i(\mathbf x^*) + \sum_j\mu_j\nabla g_j(\mathbf x^*), \qquad \mu_j \ge 0,\quad \mu_jg_j(\mathbf x^*) = 0. These are the Karush–Kuhn–Tucker (KKT) conditions (Karush 1939, Kuhn and Tucker 1951). The multipliers have a practical meaning: λi\lambda_i or μj\mu_j is the rate at which the optimal objective would improve if that constraint were relaxed by one unit — a shadow price. A zero multiplier means the constraint is not limiting the design.
Derivation
The multiplier condition, from geometry and on an example
One equality constraint. On the surface h(x)=0h(\mathbf x) = 0, a feasible move d\mathbf d must keep hh unchanged to first order: ∇hTd=0\nabla h^{\mathsf T}\mathbf d = 0, so feasible directions are perpendicular to ∇h\nabla h. At a constrained minimum, no feasible direction may decrease ff, so ∇fTd=0\nabla f^{\mathsf T}\mathbf d = 0 for all of them too. A vector perpendicular to everything that is perpendicular to ∇h\nabla h must be parallel to it:
∇f=λ ∇h.\nabla f = \lambda\,\nabla h.
The closest point on the unit circle to (2,1)(2, 1). Minimise f=(x−2)2+(y−1)2f = (x - 2)^2 + (y - 1)^2 with h=x2+y2−1=0h = x^2 + y^2 - 1 = 0. The condition reads 2(x−2)=2λx2(x - 2) = 2\lambda x and 2(y−1)=2λy2(y - 1) = 2\lambda y, so
x=21−λ,y=11−λ.x = \frac{2}{1 - \lambda}, \qquad y = \frac{1}{1 - \lambda}.
The point lies on the line through the origin and (2,1)(2, 1), as geometry demands. The constraint x2+y2=5/(1−λ)2=1x^2 + y^2 = 5/(1 - \lambda)^2 = 1 gives 1−λ=51 - \lambda = \sqrt5, so (x,y)=(2,1)/5(x, y) = (2, 1)/\sqrt5 and λ=1−5=−1.236\lambda = 1 - \sqrt5 = -1.236.
What the multiplier measures. Change the constraint to x2+y2=bx^2 + y^2 = b. The best distance squared becomes f∗(b)=(5−b)2f^*(b) = (\sqrt5 - \sqrt b)^2, and at b=1b = 1
df∗db=−5−bb=1−5=λ.\frac{df^*}{db} = -\frac{\sqrt5 - \sqrt b}{\sqrt b} = 1 - \sqrt5 = \lambda.
The multiplier is the rate at which the optimum changes as the constraint is relaxed — the shadow price. (The quadratic penalty of the lab satisfies ∇f=−μh ∇h\nabla f = -\mu h\,\nabla h at its minimum instead; comparing with ∇f=λ∇h\nabla f = \lambda\nabla h shows its constraint violation is h≈−λ/μ=1.236/μh \approx -\lambda/\mu = 1.236/\mu, the +1.22×10−2+1.22 \times 10^{-2} the lab prints at μ=100\mu = 100.)
Inequalities. For g(x)≥0g(\mathbf x) \ge 0 the feasible side is where gg increases. If the constraint is active and ∇f=μ∇g\nabla f = \mu\nabla g with μ<0\mu < 0, then stepping into the feasible region (along ∇g\nabla g) would lower ff — not a minimum. So μ≥0\mu \ge 0. And if the constraint is not active, it cannot hold the design back, so μ=0\mu = 0: in every case μg=0\mu g = 0.
📐The lightest pressure tank
A steel tank must hold 2 m³ at 3 MPa with allowable stress 200 MPa; cylinder with hemispherical heads; shell thickness from hoop stress, heads half as thick. In formulas: the shell needs t=pr/σt = pr/\sigma and the heads pr/(2σ)pr/(2\sigma), so the steel volume is 2πrL⋅pr/σ+4πr2⋅pr/(2σ)2\pi rL\cdot pr/\sigma + 4\pi r^2\cdot pr/(2\sigma) and the contents πr2L+43πr3=2\pi r^2L + \tfrac43\pi r^3 = 2 m³. Without other limits the lightest tank is a sphere of radius 0.782 m, 353 kg. But it must fit on a road vehicle: diameter at most 1.2 m. scipy.optimize.minimize(method='SLSQP') gives radius 0.6000 m (the limit), cylinder length 0.968 m, mass 417.7 kg. Checking KKT numerically: λ=235.5\lambda = 235.5 kg/m³ for the volume constraint and μ=133.2\mu = 133.2 kg/m for the diameter limit — positive, so the limit is active, and every extra metre of allowable diameter would save about 133 kg at the margin.
How constrained solvers work. A penalty method adds μ2∑hi2\tfrac\mu2\sum h_i^2 (and similar terms for violated inequalities) to ff and minimises without constraints, increasing μ\mu. It is simple but imperfect: the constraint is satisfied only as μ→∞\mu \to \infty, and the problem's Hessian condition number grows like μ\mu — at μ=104\mu = 10^4 the error is still 6×10−56 \times 10^{-5} with κ≈104\kappa \approx 10^4. Augmented Lagrangian methods add multiplier estimates and avoid infinite μ\mu. Sequential quadratic programming (SQP, e.g. SLSQP) solves at each step a quadratic model of the Lagrangian with linearised constraints — Newton's method applied to the KKT conditions — and reaches 10−910^{-9} directly. Interior-point methods keep iterates strictly feasible with a logarithmic barrier and dominate large problems.
💡Linear programming
When the objective and every constraint are linear — minimise cTx\mathbf c^{\mathsf T}\mathbf x subject to Ax≤bA\mathbf x \le \mathbf b, Aeqx=beqA_{eq}\mathbf x = \mathbf b_{eq}, x≥0\mathbf x \ge 0 — the feasible set is a convex polyhedron and an optimum, if one exists, lies at a vertex. The simplex method moves between adjacent vertices; interior-point methods (Karmarkar, 1984) cut through the middle. Both are inside scipy.optimize.linprog, which calls the HiGHS solver. Linear programmes are convex, so a solution found is the global optimum, and duality theory provides shadow prices for every constraint.
Linear programming is viewed as a revolutionary development giving man the ability to state general objectives and to find, by means of the simplex method, optimal policy decisions for a broad class of practical decision problems of great complexity.
– George Dantzig, Reminiscences about the origins of linear programming (1983)
📐The lightest truss that can carry a load
Place six nodes on a grid, connect every pair by a candidate bar (15 bars), pin two nodes and hang 100 kN from another. Split each bar force into tension t≥0t \ge 0 and compression c≥0c \ge 0; equilibrium at the free nodes is linear in (t,c)(t, c); the steel volume ∑L(t+c)/σ\sum L(t + c)/\sigma is linear too. linprog finds the minimum: a volume of 7,200 cm³ (56.5 kg of steel) using only 2 of the 15 bars — one diagonal tie in tension (223.6 kN) and one horizontal strut in compression (200 kN). The optimiser discovered the structure, not just its sizes. This 'ground structure' formulation (Dorn, Gomory and Greenberg, 1964) descends from Anthony Michell's 1904 theory of optimal frames and is the seed of modern topology optimisation.
💡Where you will use this
Optimal Guidance: powered-descent guidance for reusable rockets is solved on board as a convex optimisation problem every fraction of a second; sequential convex programming handles the nonconvex parts. Structures: sizing, shape and topology optimisation under stress and displacement constraints. Power Systems: economic dispatch and optimal power flow are LPs and SQPs. Geotechnics: limit analysis of slopes and foundations can be written as linear programmes.
📚Further reading
  • Jorge Nocedal and Stephen Wright, Numerical Optimization (2nd ed., Springer, 2006), chapters 12, 13 and 18 — the KKT conditions, the simplex method and sequential quadratic programming.
  • Stephen Boyd and Lieven Vandenberghe, Convex Optimization (Cambridge University Press, 2004), chapter 5, 'Duality' — Lagrange multipliers and shadow prices; free online.
  • George Dantzig, Linear Programming and Extensions (Princeton University Press, 1963) — by the inventor of the simplex method.
  • W. S. Dorn, R. E. Gomory and H. J. Greenberg, 'Automatic design of optimal structures', Journal de Mécanique 3, 25–52 (1964) — the ground-structure method of the truss example.
↓But there is a problem
Every number in this course so far — loads, material properties, initial conditions, measurements — has been treated as exactly known. Even least squares (Unit 5) treated noise only through the residuals it left behind. Reality is uncertain. Loads vary, materials vary, sensors drift, and a navigation system must decide where it is from measurements it knows are wrong. The last unit is about computing with uncertainty: generating randomness on a deterministic machine, propagating it through models, and combining uncertain predictions with uncertain measurements to know more than either alone. It starts with a warning from the 1960s: one of the most widely used random-number generators in the world produced numbers that were not random at all.