Equality Constrained Optimization

In previous sections, we focused on unconstrained minimization: finding $\mathbf{x} \in \mathbb{R}^n$ that minimizes an objective function $f(\mathbf{x})$ .

Many practical engineering problems require finding an optimum that satisfies physical restrictions or conservation laws expressed as equality constraints:

\min_{\mathbf{x}} f(\mathbf{x}) \quad \text{subject to} \quad \mathbf{g}(\mathbf{x}) = \mathbf{0}

(1)

Here, $\mathbf{g}(\mathbf{x}) = \mathbf{0}$ represents $m$ simultaneous scalar constraints ( $m < n$ ). For example, constraining a 2D trajectory to the unit circle:

\begin{aligned} x^2 + y^2 &= 1 \\ g(x, y) = x^2 + y^2 - 1 &= 0 \end{aligned}

(2)

import numpy as np
import scipy as sp
from scipy.optimize import minimize, NonlinearConstraint
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

from executable_engineering import plot_constrained_optimum_3d

Penalty Method¶

A straightforward heuristic for enforcing constraints is to augment the objective function with a penalty term proportional to the constraint violation:

\min_{\mathbf{x}} f(\mathbf{x}) = f(\mathbf{x}) + \lambda \|\mathbf{g}(\mathbf{x})\|^2

(3)

where $\lambda > 0$ is a penalty parameter. When the constraint is violated ( $\mathbf{g}(\mathbf{x}) \ne \mathbf{0}$ ), the penalty increases the effective cost.

While intuitive, penalty methods face severe practical challenges:

Too small $\lambda$ : The optimizer drifts away from the constraint surface, yielding infeasible solutions.
Too large $\lambda$ : The penalty introduces steep, narrow canyon walls into the objective landscape, resulting in an ill-conditioned Hessian ( $\mathbf{H} \sim \lambda$ ) that stalls gradient-based and Newton solvers.

Sequential unconstrained minimization techniques (SUMT) adjust $\lambda$ iteratively, but balancing feasibility against numerical conditioning remains challenging.

Method of Undetermined Lagrange Multipliers¶

At the constrained optimum, the constraint is satisfied exactly ( $\mathbf{g}(\mathbf{x}^*) = \mathbf{0}$ ). Instead of choosing an arbitrary scalar penalty $\lambda$ , we treat the multiplier $\boldsymbol{\lambda}$ as an unknown variable to be determined simultaneously with $\mathbf{x}$ .

We formulate the Lagrangian function:

\mathcal{L}(\mathbf{x}, \boldsymbol{\lambda}) = f(\mathbf{x}) + \boldsymbol{\lambda}^T \mathbf{g}(\mathbf{x})

(4)

The necessary condition for a constrained stationary point is that the gradient of the Lagrangian with respect to both $\mathbf{x}$ and $\boldsymbol{\lambda}$ vanishes:

\nabla \mathcal{L}(\mathbf{x}, \boldsymbol{\lambda}) = \mathbf{0}

(5)

Differentiating with respect to the multipliers $\boldsymbol{\lambda}$ recovers the constraint equations:

\frac{\partial \mathcal{L}}{\partial \boldsymbol{\lambda}} = \mathbf{g}(\mathbf{x}) = \mathbf{0}

(6)

Remarkably, the Lagrangian formulation converts a constrained optimization problem into an unconstrained problem $\nabla_{\mathbf{x,\lambda}} \mathcal{L} = \mathbf{0}$ which implies our normal optimization methods apply (or one of them will anyway)!

Geometric Interpretation of the Optimality Condition¶

Differentiating the Lagrangian with respect to $\mathbf{x}$ yields:

\frac{\partial \mathcal{L}}{\partial \mathbf{x}} = \nabla f(\mathbf{x}) + \sum_{i=1}^m \lambda_i \nabla g_i(\mathbf{x}) = \mathbf{0} \implies \nabla f(\mathbf{x}) = -\sum_{i=1}^m \lambda_i \nabla g_i(\mathbf{x})

(7)

Along the constraint curve $\mathbf{g}(\mathbf{x}) = \mathbf{0}$ , the gradient of the objective function $\nabla f$ must be parallel (or anti-parallel) to the normal of the constraint curve $\nabla g$ . If $\nabla f$ had any non-zero component tangential to the constraint curve, we could decrease $f(\mathbf{x})$ further by moving along that curve without violating the constraint!

Hiking Analogy: Imagine walking along a narrow mountain trail defined by $g(x, y) = 0$ while the mountain height is $f(x, y)$ . To find the lowest point on the trail, you hike downhill along the path. At the lowest path elevation, the ground slopes perpendicularly away from the trail ( $\nabla f \parallel \nabla g$ ). To prevent falling down the slope, you must push against the slope with constraint force $\lambda$ .

Quadratic Programming with Linear Equality Constraints¶

If we blindly apply Lagrange multipliers, our methods will likely fail. To understand why, consider a quadratic objective and linear constraints (i.e.: a 2nd order Taylor expansion). The Lagrangian is quadratic:

f(\mathbf{x}) = \frac{1}{2} \mathbf{x}^T \mathbf{H} \mathbf{x} + \mathbf{c}^T \mathbf{x}, \quad \mathbf{g}(\mathbf{x}) = \mathbf{A} \mathbf{x} - \mathbf{b} = \mathbf{0}

(8)

\mathcal{L} = \frac{1}{2} \mathbf{x}^T \mathbf{H} \mathbf{x} + \mathbf{c}^T \mathbf{x} +\boldsymbol{\lambda} (\mathbf{A} \mathbf{x} - \mathbf{b})

(9)

The stationary conditions $\nabla_{\mathbf{x}} \mathcal{L} = \mathbf{0}$ and $\nabla_{\boldsymbol{\lambda}} \mathcal{L} = \mathbf{0}$ become:

\begin{aligned} \mathbf{H} \mathbf{x} + \mathbf{c} + \mathbf{A}^T \boldsymbol{\lambda} &= \mathbf{0} \\ \mathbf{A} \mathbf{x} - \mathbf{b} &= \mathbf{0} \end{aligned}

(10)

Writing this as a block system yields the symmetric Karush-Kuhn-Tucker (KKT) linear system: $$\begin{bmatrix} \mathbf{H} & \mathbf{A}^T \ \mathbf{A} & \mathbf{0} \end{bmatrix}

\begin{bmatrix} \mathbf{x} \\ \boldsymbol{\lambda} \end{bmatrix}

(11)

=

\begin{bmatrix} -\mathbf{c} \\ \mathbf{b} \end{bmatrix}

(12)

Remarkably, solving a linearly constrained quadratic program reduces to solving a single square linear system!

Example: Quadratic KKT System¶

Minimize the objective $f(x, y, z) = xy + yz$ subject to two linear constraints:

x + 2y = 6 \quad \text{and} \quad x - 3z = 0

(13)

Here $\mathbf{x} = [x, y, z]^T$ , $\mathbf{c} = [0, 0, 0]^T$ , and:

\mathbf{H} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \quad \mathbf{A} = \begin{bmatrix} 1 & 2 & 0 \\ 1 & 0 & -3 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} 6 \\ 0 \end{bmatrix}

(14)

Assembling the $5 \times 5$ KKT matrix $\mathbf{K} = \begin{bmatrix} \mathbf{H} & \mathbf{A}^T \\ \mathbf{A} & \mathbf{0} \end{bmatrix}$ , we solve $\mathbf{K} \begin{bmatrix} \mathbf{x} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} -\mathbf{c} \\ \mathbf{b} \end{bmatrix}$ .

K = np.array([
    [0, 1, 0, 1, 1],
    [1, 0, 1, 2, 0],
    [0, 1, 0, 0, -3],
    [1, 2, 0, 0, 0],
    [1, 0, -3, 0, 0]
], dtype=float)

rhs = np.array([0, 0, 0, 6, 0], dtype=float)
sol = sp.linalg.solve(K, rhs)

print(f"Optimal coordinates: x = {sol[0]:.4f}, y = {sol[1]:.4f}, z = {sol[2]:.4f}")
print(f"Lagrange multipliers: lambda_1 = {sol[3]:.4f}, lambda_2 = {sol[4]:.4f}")
print(f"Objective value f(x, y, z): {sol[0]*sol[1] + sol[1]*sol[2]:.4f}")

Saddle-Point Structure of the KKT Matrix¶

Inspecting the eigenvalues of the KKT matrix reveals a fundamental property of constrained optimization.

Because the eigenvalues have opposite signs (both positive and negative), the KKT system is indefinite. The constrained optimum is not a local minimum of $\mathcal{L}(\mathbf{x}, \boldsymbol{\lambda})$ , but a saddle point (a minimum with respect to $\mathbf{x}$ , and a maximum with respect to $\boldsymbol{\lambda}$ ).

Standard unconstrained descent methods fail when applied naively to the Lagrangian. However, the block structure of the KKT matrix makes it directly amenable to block elimination and specialized constrained solvers.

evals = sp.linalg.eigvals(K)
print(f"KKT Eigenvalues: {np.sort(np.round(np.real(evals), 2))}")

Meaning of the Lagrange Multipliers¶

Lagrange multipliers $\lambda_i$ represent the sensitivity (marginal rate of change) of the optimal objective value $f^*$ with respect to relaxing the $i$ -th constraint:

\lambda_i^* = -\frac{\partial f^*}{\partial b_i}

(15)

Depending on the domain, $\lambda_i^*$ carries profound physical and economic meaning:

Structural & Mechanical Engineering: $\lambda_i$ corresponds to the constraint reaction force (or internal stress) required to maintain the boundary condition $g_i(\mathbf{x}) = 0$ .
Economics & Operations Research: $\lambda_i$ is the shadow price or marginal cost of relaxing a resource limit by one unit.
Thermodynamics & Statistical Physics: Multipliers represent conjugate thermodynamic potentials (e.g., temperature $\beta = 1/k_B T$ enforcing average energy conservation, chemical potential $\mu$ enforcing particle conservation).

General Constrained Optimization¶

For nonlinear objective functions subject to nonlinear equality or inequality constraints:

\min_{\mathbf{x}} f(\mathbf{x}) \quad \text{subject to} \quad \mathbf{g}(\mathbf{x}) = \mathbf{0}, \quad \mathbf{h}(\mathbf{x}) \le \mathbf{0}

(16)

Modern numerical solvers employ two primary families of algorithms:

Sequential Quadratic Programming (SQP): Solves a sequence of quadratic programming subproblems, iteratively approximating the Hessian of the Lagrangian and linearizing constraints (e.g., SLSQP in SciPy).
Interior-Point (Barrier) / Trust-Region Methods: Traverses the interior of the feasible region using a barrier function for inequalities while enforcing equality constraints via trust-region steps (e.g., trust-constr in SciPy).

Example: Minimizing $x^4 + y^4$ on the Unit Circle¶

Consider minimizing the quartic objective $f(x, y) = x^4 + y^4$ subject to the equality constraint $g(x, y) = x^2 + y^2 - 1 = 0$ (the unit circle).

The unconstrained minimum is at the origin $(0, 0)$ , where $f(0, 0) = 0$ . On the unit circle $x^2 + y^2 = 1$ , the function achieves four symmetric minima at $(\pm 1/\sqrt{2}, \pm 1/\sqrt{2}) \approx (\pm 0.707, \pm 0.707)$ with $f^* = 0.5$ .

We solve this using SciPy’s Sequential Least Squares Programming (SLSQP) optimizer and visualize the surface, constraint ring, and optimum.

def objective(x):
    return x[0]**4 + x[1]**4

def constraint(x):
    return x[0]**2 + x[1]**2 - 1.0

x0 = [0.5, 0.5]
cons = {'type': 'eq', 'fun': constraint}
sol = minimize(objective, x0, method='SLSQP', constraints=cons)

print(f"Status: {sol.message}")
print(f"Optimal solution: x = {sol.x[0]:.4f}, y = {sol.x[1]:.4f}")
print(f"Objective value: f(x*) = {sol.fun:.4f}")

theta = np.linspace(0, 2 * np.pi, 200)
circle = np.column_stack([np.cos(theta), np.sin(theta)])
plot_constrained_optimum_3d(
    objective, sol.x, constraint_curve=circle, x0=x0,
    title="SLSQP: Minimizing x^4 + y^4 on Unit Circle"
)

Example: Constraining an Indefinite Saddle System¶

Next, consider the hyperbolic saddle objective $f(x, y) = x^2 - 6xy + y^2$ subject to the unit circle constraint $x^2 + y^2 = 1$ .

In unconstrained space, the origin $(0, 0)$ is a saddle point where the Hessian has eigenvalues of opposite signs ( $\lambda_1 = 8, \lambda_2 = -4$ ). When restricted to the unit circle, the objective attains constrained minima at $(\pm 1/\sqrt{2}, \pm 1/\sqrt{2})$ with $f^* = -2$ .

We solve this using SciPy’s trust-region constrained optimizer (trust-constr).

def saddle_obj(x):
    return x[0]**2 - 6*x[0]*x[1] + x[1]**2

cons_saddle = {'type': 'eq', 'fun': constraint}
sol_saddle = minimize(saddle_obj, [0.1, 0.1], method='trust-constr', constraints=cons_saddle)

print(f"Status: {sol_saddle.message}")
print(f"Optimal solution: x = {sol_saddle.x[0]:.4f}, y = {sol_saddle.x[1]:.4f}")
print(f"Objective value: f(x*) = {sol_saddle.fun:.4f}")

plot_constrained_optimum_3d(
    saddle_obj, sol_saddle.x, constraint_curve=circle, x0=[0.1, 0.1],
    title="trust-constr: Saddle x^2 - 6xy + y^2 on Unit Circle"
)

Example: Inequality Constraints on the Unit Disk¶

Finally, consider minimizing $f(x, y) = \sin(x)\cos(y) + x^2 + y^2$ subject to an inequality constraint:

x^2 + y^2 \le 1 \iff x^2 + y^2 - 1 \le 0

(17)

The feasible region is the closed unit disk. When optimizing with inequalities, the constraint may be:

Active: The unconstrained minimum lies outside the disk; the constrained optimum sits on the boundary $x^2 + y^2 = 1$ .
Inactive: The unconstrained minimum lies strictly inside the disk; the multiplier $\lambda = 0$ .

We formulate this using SciPy’s NonlinearConstraint(fun, -np.inf, 0) with trust-constr.

def nonlinear_obj(x):
    return np.sin(x[0]) * np.cos(x[1]) + x[0]**2 + x[1]**2

disk_constraint = NonlinearConstraint(lambda x: x[0]**2 + x[1]**2 - 1.0, -np.inf, 0.0)
sol_disk = minimize(nonlinear_obj, [0.5, 0.5], method='trust-constr', constraints=disk_constraint)

print(f"Status: {sol_disk.message}")
print(f"Optimal solution: x = {sol_disk.x[0]:.4f}, y = {sol_disk.x[1]:.4f}")
print(f"Distance from origin: r = {np.linalg.norm(sol_disk.x):.4f}")
print(f"Objective value: f(x*) = {sol_disk.fun:.4f}")

plot_constrained_optimum_3d(
    nonlinear_obj, sol_disk.x, constraint_curve=circle, x0=[0.5, 0.5],
    title="trust-constr: sin(x)cos(y) + x^2 + y^2 on Unit Disk"
)