1"""
2Constrained and Unconstrained Optimization
3
4This script demonstrates:
51. Unconstrained optimization using scipy.optimize
62. Lagrange multipliers for equality constraints
73. KKT conditions for inequality constraints
84. Convex vs non-convex optimization comparison
9
10Author: Math for AI Examples
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.optimize import minimize, NonlinearConstraint
16from typing import Callable, Tuple
17
18
19def unconstrained_optimization_demo():
20 """Demonstrate unconstrained optimization with scipy.optimize."""
21 print("\n" + "="*60)
22 print("1. Unconstrained Optimization")
23 print("="*60)
24
25 # Objective: Minimize f(x,y) = x^2 + y^2 + 2x + 4y + 5
26 # Analytical solution: x = -1, y = -2, f_min = 0
27 def objective(x):
28 return x[0]**2 + x[1]**2 + 2*x[0] + 4*x[1] + 5
29
30 def gradient(x):
31 return np.array([2*x[0] + 2, 2*x[1] + 4])
32
33 x0 = np.array([5.0, 5.0])
34
35 # Using different optimization methods
36 methods = ['BFGS', 'CG', 'Newton-CG', 'L-BFGS-B']
37 results = {}
38
39 for method in methods:
40 if method == 'Newton-CG':
41 result = minimize(objective, x0, method=method, jac=gradient)
42 else:
43 result = minimize(objective, x0, method=method)
44
45 results[method] = result
46 print(f"\n{method}:")
47 print(f" Optimal point: ({result.x[0]:.6f}, {result.x[1]:.6f})")
48 print(f" Minimum value: {result.fun:.6f}")
49 print(f" Iterations: {result.nit}")
50
51 print("\nAnalytical solution: x = -1.0, y = -2.0, f_min = 0.0")
52
53
54def lagrange_multipliers_demo():
55 """
56 Demonstrate Lagrange multipliers for equality constraints.
57
58 Problem: Minimize f(x,y) = x^2 + y^2
59 Subject to: g(x,y) = x + y - 1 = 0
60
61 Using Lagrange multipliers:
62 L(x,y,λ) = f(x,y) - λ*g(x,y) = x^2 + y^2 - λ(x + y - 1)
63
64 KKT conditions:
65 ∂L/∂x = 2x - λ = 0 => x = λ/2
66 ∂L/∂y = 2y - λ = 0 => y = λ/2
67 ∂L/∂λ = -(x + y - 1) = 0 => x + y = 1
68
69 Solution: x = y = 0.5, f_min = 0.5
70 """
71 print("\n" + "="*60)
72 print("2. Lagrange Multipliers (Equality Constraints)")
73 print("="*60)
74 print("\nProblem: Minimize f(x,y) = x^2 + y^2")
75 print(" Subject to: x + y = 1")
76
77 # Manual solution using Lagrange multipliers
78 print("\n--- Analytical Solution (Lagrange Multipliers) ---")
79 print("Setting up Lagrangian: L(x,y,λ) = x^2 + y^2 - λ(x + y - 1)")
80 print("∂L/∂x = 2x - λ = 0")
81 print("∂L/∂y = 2y - λ = 0")
82 print("∂L/∂λ = -(x + y - 1) = 0")
83 print("\nFrom first two equations: x = y = λ/2")
84 print("From constraint: λ/2 + λ/2 = 1 => λ = 1")
85 print("Solution: x = 0.5, y = 0.5, f_min = 0.5")
86
87 # Numerical solution using scipy
88 print("\n--- Numerical Solution (scipy.optimize) ---")
89
90 def objective(x):
91 return x[0]**2 + x[1]**2
92
93 def constraint(x):
94 return x[0] + x[1] - 1
95
96 constraints = {'type': 'eq', 'fun': constraint}
97 x0 = np.array([2.0, 2.0])
98
99 result = minimize(objective, x0, method='SLSQP', constraints=constraints)
100
101 print(f"Optimal point: ({result.x[0]:.6f}, {result.x[1]:.6f})")
102 print(f"Minimum value: {result.fun:.6f}")
103 print(f"Constraint satisfied: {constraint(result.x):.6e}")
104
105 # Visualization
106 fig, ax = plt.subplots(figsize=(8, 8))
107
108 # Plot objective function contours
109 x = np.linspace(-0.5, 2, 300)
110 y = np.linspace(-0.5, 2, 300)
111 X, Y = np.meshgrid(x, y)
112 Z = X**2 + Y**2
113
114 contour = ax.contour(X, Y, Z, levels=15, cmap='viridis', alpha=0.6)
115 ax.clabel(contour, inline=True, fontsize=8)
116
117 # Plot constraint line
118 x_line = np.linspace(-0.5, 2, 100)
119 y_line = 1 - x_line
120 ax.plot(x_line, y_line, 'r-', linewidth=3, label='Constraint: x + y = 1')
121
122 # Plot optimal point
123 ax.plot(0.5, 0.5, 'r*', markersize=20, label='Optimal Point (0.5, 0.5)')
124
125 ax.set_xlabel('x', fontsize=12)
126 ax.set_ylabel('y', fontsize=12)
127 ax.set_title('Lagrange Multipliers: Min x² + y² s.t. x + y = 1',
128 fontsize=14, fontweight='bold')
129 ax.legend(fontsize=11)
130 ax.grid(True, alpha=0.3)
131 ax.set_aspect('equal')
132
133 plt.tight_layout()
134 plt.savefig('lagrange_multipliers.png', dpi=150, bbox_inches='tight')
135 print("\nSaved visualization to 'lagrange_multipliers.png'")
136
137
138def kkt_conditions_demo():
139 """
140 Demonstrate KKT conditions for inequality constraints.
141
142 Problem: Minimize f(x,y) = (x-2)^2 + (y-1)^2
143 Subject to: x + y <= 1 (g1)
144 x >= 0 (g2)
145 y >= 0 (g3)
146
147 KKT conditions:
148 1. Stationarity: ∇f(x*) + Σμ_i*∇g_i(x*) = 0
149 2. Primal feasibility: g_i(x*) <= 0
150 3. Dual feasibility: μ_i >= 0
151 4. Complementary slackness: μ_i * g_i(x*) = 0
152 """
153 print("\n" + "="*60)
154 print("3. KKT Conditions (Inequality Constraints)")
155 print("="*60)
156 print("\nProblem: Minimize f(x,y) = (x-2)^2 + (y-1)^2")
157 print(" Subject to: x + y <= 1, x >= 0, y >= 0")
158
159 def objective(x):
160 return (x[0] - 2)**2 + (x[1] - 1)**2
161
162 # Constraints: g(x) <= 0
163 constraints = [
164 {'type': 'ineq', 'fun': lambda x: 1 - x[0] - x[1]}, # x + y <= 1
165 {'type': 'ineq', 'fun': lambda x: x[0]}, # x >= 0
166 {'type': 'ineq', 'fun': lambda x: x[1]} # y >= 0
167 ]
168
169 x0 = np.array([0.5, 0.5])
170 result = minimize(objective, x0, method='SLSQP', constraints=constraints)
171
172 print(f"\nOptimal point: ({result.x[0]:.6f}, {result.x[1]:.6f})")
173 print(f"Minimum value: {result.fun:.6f}")
174
175 # Check which constraints are active
176 print("\nConstraint activity:")
177 print(f" x + y - 1 = {result.x[0] + result.x[1] - 1:.6f} (active if ≈ 0)")
178 print(f" x = {result.x[0]:.6f} (active if ≈ 0)")
179 print(f" y = {result.x[1]:.6f} (active if ≈ 0)")
180
181 # Visualization
182 fig, ax = plt.subplots(figsize=(8, 8))
183
184 x = np.linspace(-0.5, 2.5, 300)
185 y = np.linspace(-0.5, 2.5, 300)
186 X, Y = np.meshgrid(x, y)
187 Z = (X - 2)**2 + (Y - 1)**2
188
189 contour = ax.contour(X, Y, Z, levels=20, cmap='viridis', alpha=0.6)
190 ax.clabel(contour, inline=True, fontsize=8)
191
192 # Plot feasible region
193 x_constraint = np.linspace(0, 1, 100)
194 y_constraint = 1 - x_constraint
195
196 ax.fill([0, 1, 0, 0], [0, 0, 1, 0], color='lightblue',
197 alpha=0.3, label='Feasible Region')
198 ax.plot(x_constraint, y_constraint, 'r-', linewidth=2, label='x + y = 1')
199 ax.axhline(y=0, color='k', linestyle='--', linewidth=1)
200 ax.axvline(x=0, color='k', linestyle='--', linewidth=1)
201
202 # Plot unconstrained minimum
203 ax.plot(2, 1, 'ko', markersize=10, label='Unconstrained Min (2, 1)')
204
205 # Plot constrained minimum
206 ax.plot(result.x[0], result.x[1], 'r*', markersize=20,
207 label=f'Constrained Min ({result.x[0]:.2f}, {result.x[1]:.2f})')
208
209 ax.set_xlabel('x', fontsize=12)
210 ax.set_ylabel('y', fontsize=12)
211 ax.set_title('KKT Conditions: Inequality Constraints',
212 fontsize=14, fontweight='bold')
213 ax.legend(fontsize=10)
214 ax.grid(True, alpha=0.3)
215 ax.set_xlim(-0.5, 2.5)
216 ax.set_ylim(-0.5, 2.5)
217
218 plt.tight_layout()
219 plt.savefig('kkt_conditions.png', dpi=150, bbox_inches='tight')
220 print("\nSaved visualization to 'kkt_conditions.png'")
221
222
223def convex_vs_nonconvex():
224 """Compare optimization on convex vs non-convex functions."""
225 print("\n" + "="*60)
226 print("4. Convex vs Non-Convex Optimization")
227 print("="*60)
228
229 # Convex function: f(x,y) = x^2 + y^2
230 def convex_func(x):
231 return x[0]**2 + x[1]**2
232
233 # Non-convex function: Himmelblau's function
234 # f(x,y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2
235 # Has 4 local minima, all with same value
236 def nonconvex_func(x):
237 return (x[0]**2 + x[1] - 11)**2 + (x[0] + x[1]**2 - 7)**2
238
239 # Test multiple starting points
240 starting_points = [
241 np.array([0.0, 0.0]),
242 np.array([2.0, 2.0]),
243 np.array([-2.0, 2.0]),
244 np.array([2.0, -2.0])
245 ]
246
247 print("\n--- Convex Function: f(x,y) = x^2 + y^2 ---")
248 print("Expected: All starting points converge to (0, 0)")
249
250 for i, x0 in enumerate(starting_points):
251 result = minimize(convex_func, x0, method='BFGS')
252 print(f"Start {i+1} {x0}: -> ({result.x[0]:.4f}, {result.x[1]:.4f}), "
253 f"f = {result.fun:.4f}")
254
255 print("\n--- Non-Convex Function: Himmelblau's Function ---")
256 print("Expected: Different starting points may converge to different minima")
257 print("Known minima: (3, 2), (-2.805, 3.131), (-3.779, -3.283), (3.584, -1.848)")
258
259 for i, x0 in enumerate(starting_points):
260 result = minimize(nonconvex_func, x0, method='BFGS')
261 print(f"Start {i+1} {x0}: -> ({result.x[0]:.4f}, {result.x[1]:.4f}), "
262 f"f = {result.fun:.4f}")
263
264 # Visualization
265 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
266
267 x = np.linspace(-5, 5, 300)
268 y = np.linspace(-5, 5, 300)
269 X, Y = np.meshgrid(x, y)
270
271 # Convex function
272 Z1 = X**2 + Y**2
273 contour1 = ax1.contour(X, Y, Z1, levels=20, cmap='viridis')
274 ax1.clabel(contour1, inline=True, fontsize=8)
275 ax1.plot(0, 0, 'r*', markersize=20, label='Global Minimum')
276 ax1.set_xlabel('x', fontsize=12)
277 ax1.set_ylabel('y', fontsize=12)
278 ax1.set_title('Convex Function: x² + y²', fontsize=14, fontweight='bold')
279 ax1.legend(fontsize=11)
280 ax1.grid(True, alpha=0.3)
281
282 # Non-convex function
283 Z2 = (X**2 + Y - 11)**2 + (X + Y**2 - 7)**2
284 contour2 = ax2.contour(X, Y, Z2, levels=30, cmap='viridis')
285 ax2.clabel(contour2, inline=True, fontsize=8)
286
287 # Plot all four minima
288 minima = [(3, 2), (-2.805, 3.131), (-3.779, -3.283), (3.584, -1.848)]
289 for xm, ym in minima:
290 ax2.plot(xm, ym, 'r*', markersize=15)
291
292 ax2.set_xlabel('x', fontsize=12)
293 ax2.set_ylabel('y', fontsize=12)
294 ax2.set_title("Non-Convex: Himmelblau's Function", fontsize=14, fontweight='bold')
295 ax2.grid(True, alpha=0.3)
296
297 plt.tight_layout()
298 plt.savefig('convex_vs_nonconvex.png', dpi=150, bbox_inches='tight')
299 print("\nSaved comparison plot to 'convex_vs_nonconvex.png'")
300
301
302if __name__ == "__main__":
303 print("="*60)
304 print("Constrained and Unconstrained Optimization Examples")
305 print("="*60)
306
307 # Run demonstrations
308 unconstrained_optimization_demo()
309 lagrange_multipliers_demo()
310 kkt_conditions_demo()
311 convex_vs_nonconvex()
312
313 print("\n" + "="*60)
314 print("Key Takeaways:")
315 print("="*60)
316 print("1. Unconstrained: Use gradient-based methods (BFGS, CG, L-BFGS)")
317 print("2. Equality constraints: Lagrange multipliers convert to unconstrained")
318 print("3. Inequality constraints: KKT conditions generalize Lagrange multipliers")
319 print("4. Convex problems: Guarantee global optimum, efficient algorithms")
320 print("5. Non-convex problems: May have multiple local minima, need global search")