diff --git a/class01/Manifest.toml b/class01/Manifest.toml index dee6afe..dbb8e20 100644 --- a/class01/Manifest.toml +++ b/class01/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.5" manifest_format = "2.0" -project_hash = "214e4ad4bd7d700b3029139dbd6c6e64d5c86ac2" +project_hash = "33fab7715768642c9fee44610a16a04e852547fe" [[deps.ADTypes]] git-tree-sha1 = "7927b9af540ee964cc5d1b73293f1eb0b761a3a1" diff --git a/class01/Project.toml b/class01/Project.toml index 4eb6976..caf9860 100644 --- a/class01/Project.toml +++ b/class01/Project.toml @@ -24,4 +24,5 @@ ShortCodes = "f62ebe17-55c5-4640-972f-b59c0dd11ccf" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +WebIO = "0f1e0344-ec1d-5b48-a673-e5cf874b6c29" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/class01/background_materials/optimization_basics.jl b/class01/background_materials/optimization_basics.jl index 9757985..513da66 100644 --- a/class01/background_materials/optimization_basics.jl +++ b/class01/background_materials/optimization_basics.jl @@ -19,7 +19,7 @@ begin using PlutoTeachingTools using ShortCodes, MarkdownLiteral using Random - using Plots + using CairoMakie Random.seed!(8803) end @@ -28,18 +28,31 @@ begin using JuMP using HiGHS # Solver for LPs and MILPs using Ipopt # Solver for NLPs - using Test # For quick validation checks - # const MOI = JuMP.MathOptInterface end -# ╔═╡ 0df8b65a-0527-4545-bf11-00e9912bced0 +# ╔═╡ b879c4bf-8292-45a0-87d0-b1f0fa456585 md""" - | | | | |-----------:|:--|:------------------| | Lecturer | : | Rosemberg, Andrew | | Date | : | 28 of July, 2025 | +""" + +# ╔═╡ eeceb82e-abfb-4502-bcfb-6c9f76a0879d +md""" +## Prelude + +Research and make sure you understand the following concepts/questions: + - What is a convex maximization optimziation problem? and a minimization one? Why are convex problems desirable (compared to non-convex ones)? + - What is a linear program? + - What are Lagrangian multipliers? + - What are Integer problems? and Mixed-Integer problems? + - What is duality? Can you construct the dual of at least a linear program? +""" + +# ╔═╡ 0df8b65a-0527-4545-bf11-00e9912bced0 +md""" # Background – Modeling Optimization Problems in JuMP 🏗️ This short Pluto notebook walks you through three small optimisation models of increasing @@ -54,7 +67,7 @@ For every task you will: * Write down the mathematical formulation. * Translate it into a JuMP model. * Solve it. -* Run the provided `@testset` to make sure your implementation is correct. +* Pluto will run tests to check your answer. When the tests are green ✅ you can be confident that your model is producing the expected answer. """ @@ -64,9 +77,9 @@ md""" ## 1. Linear program – Production planning -A workshop makes **widgets** \(w\) and **gadgets** \(g\). +A workshop makes **widgets** \(w\) and **gadgets** \(g\). The table below specifies how many machine--hours and human labour hours are required to produce each product and how much profit the company makes after selling the product. -| | Machine‑hours | Labour‑hours | Profit (\$) | +| | Machine‑-hours | Labour‑-hours | Profit (\$) | |------------|---------------|--------------|--------------| | Widget (\(w\)) | 2 | 3 | 3 | | Gadget (\(g\)) | 4 | 2 | 5 | @@ -74,9 +87,8 @@ A workshop makes **widgets** \(w\) and **gadgets** \(g\). Resources available this week: **100 machine‑hours** and **90 labour‑hours**. ### 1.1  Your tasks -1. Write the mathematical model *(maximization)*. +1. Write the mathematical model that maximizes profit! Let's assume we can sell continuous fractions of our produts. 2. Fill in the JuMP code in the next cell. -3. Run the tests. """ # ╔═╡ 49042d6c-cf78-46d3-bfee-a8fd7ddf3aa0 @@ -85,31 +97,238 @@ begin # Replace the contents of this cell with your own model. model_lp = Model(HiGHS.Optimizer) -# Suggested variable names -# @variable(model_lp, w >= 0) -# @variable(model_lp, g >= 0) +# Required variable names (used for testing) +@variable(model_lp, w >= 0) # number of Widgets +@variable(model_lp, g >= 0) # number of Gadgets # --- YOUR CODE HERE --- -# optimize!(model_lp) +# optimize!(model_lp) # uncomment to optimize + +# Let's look at our model +println(model_lp) end -# ╔═╡ 6fb672d0-5a18-4ccc-b7b3-184839c2401b +# ╔═╡ 1d3edbdd-7747-4651-b650-c6b9bf87b460 +md"Tests will automatically fetch the optimal values from your solved model." + +# ╔═╡ 248b398a-0cf5-4c2b-8752-7b9cc4e765d6 +question_box(md"Did we get partial products?") + +# ╔═╡ 245eb671-84e1-447b-8045-e9eb04966d80 +md""" +### 1.2 Other Modeling Tricks +#### Epigraph reformulation for absolute-value expressions +A convex **epigraph** turns a non-linear term into a linear one by treating it +as the upper envelope of an auxiliary variable. + +*Goal:* model $u \;\ge\; |x|$ (or set $u=|x|$ in an objective). + +**Epigraph form** + +```math +\begin{aligned} +u &\ge \phantom{-}x,\\[-2pt] +u &\ge -x. +\end{aligned} +``` + +* Both constraints are linear, so $u$ is the **least upper bound** on $|x|$. +* If the objective minimises $u$ (or if $u$ appears on the left-hand side of + other $\le$ constraints) the optimum forces $u=|x|$ automatically. +* No binaries are required – this is purely continuous and exploits the convex + shape of $|x|$. + +--- + +""" + +# ╔═╡ 6a823649-04fa-4322-a028-2fb29dffb08b +md""" +**Example** – Weighted-Median Warehouse Location on a Highway + +A logistics company must decide where to build a single cross-dock hub along a straight 50-km stretch of highway. +Eight supermarkets lie on the same axis at known kilometre posts $d_i$ and generate a daily demand (truck-loads) $w_i$. + +| Supermarket $i$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | +|------------------------------|---|---|---|---|---|---|---|---| +| Location $d_i\,[\text{km}]$| 3 | 8 | 15 | 20 | 25 | 29 | 35 | 40 | +| Daily demand $w_i$ | 2 | 1 | 3 | 4 | 2 | 3 | 5 | 1 | + +Let $x \in [0,50]$ be the kilometre post chosen for the hub. +Travel occurs on the highway, so distance is **one-dimensional**. +Devise an optimization problem that minimizes the work done by the truck (distance times weight). + +""" + +# ╔═╡ c369ab46-b416-4c12-83fe-65040a0c47c8 begin -# === Quick check === -@testset "LP check" begin - @test termination_status(model_lp) == MOI.OPTIMAL - @test isapprox(objective_value(model_lp), 135.0; atol = 1e-3) -end +# === Your LP model goes below === +# Replace the contents of this cell with your own model. +model_lp2 = Model(HiGHS.Optimizer) + +# Required variable names (used for testing) +@variable(model_lp2, x_d) + +# --- YOUR CODE HERE --- + +# optimize!(model_lp2) # uncomment to optimize + +# Let's look at our model +println(model_lp2) end +# ╔═╡ b13f9775-68c2-4646-9b67-c69ee23a4ea0 +md""" +#### Cutting-plane (outer-approximation) method for nonlinear functions +When a constraint contains a smooth convex function $g(x)\le 0$, one can +iteratively approximate it by linear cuts that **support** the graph of $g$. + +**Idea** + +1. **Relax** the nonlinear constraint and solve the LP/MILP master problem. +2. **Check** the candidate solution $\bar x$: + * If $g(\bar x)\le0$ it is feasible (and the current best). + * Otherwise, generate a tangent plane that cuts off $\bar x$ but remains + valid for all feasible points (by convexity). + + The gradient $\nabla g(\bar x)$ gives the supporting hyperplane + +```math + g(\bar x) + \nabla g(\bar x)^{\!\top}(x-\bar x)\;\le\;0 . +``` + +3. **Add** this cut to the master problem and repeat. + +Because each cut is a *global* under-estimator of $g$, the feasible region is +shrunk safely from the outside—hence *outer approximation*. For a convex NLP +the procedure converges finitely to the exact optimum; for non-convex problems +it provides a lower bound similar to branch-and-cut. + +*Common uses* + +* **Piecewise-linearisation** of $\sqrt{\cdot},\;\log(\cdot),\;\exp(\cdot)$ or + power losses in energy networks. +* **Benders-like** decomposition where cuts approximate a difficult sub-problem + (e.g. chance-constraints, stochastic recourse). + +--- + +Combined with the earlier tricks, epigraphs give **exact linear models** for +many convex terms, while cutting planes let you tackle the truly nonlinear bits +progressively, keeping the master problem linear and solver-friendly. +""" + +# ╔═╡ ea3ea95a-58cb-4d0d-a167-aa68b8bc2645 +md""" +**Example** -- Cutting-plane Exercise – Approximating a Circular Constraint with Linear Cuts + +A small drone must carry a **payload** of mass $m\,[\mathrm{kg}]$ and an +**on-board battery** of capacity $c\,[\mathrm{Ah}]$. +Because of fuselage geometry the two design variables must lie +*inside* a circular envelope in the $(m,c)$-plane: + +```math +m^{2} + c^{2} \; \le \; 100 . +``` + +The engineering goal is to minimise the total cost + +```math +\min\; 200\,m + 80\,c , +``` + +subject to the circle above and simple bounds + +```math +0 \;\le\; m \;\le\; 12, \qquad +0 \;\le\; c \;\le\; 12 . +``` + +Since the solver available to you handles only linear (M)ILPs, you will +replace the quadratic constraint by a sequence of cutting planes +(outer approximation). + +**Cutting-plane algorithm outline** + +1. **Master LP** + * Start with the bounding box only. + * Solve for a tentative point $(\bar m,\bar c)$. + +2. **Feasibility check** + * If $\bar m^{2}+\bar c^{2} \le 100$ the point is feasible $\Rightarrow$ **done**. + * Otherwise generate a *supporting hyperplane* for the circle at $(\bar m,\bar c)$: + +```math +2\,\bar m\,(m - \bar m) \;+\; 2\,\bar c\,(c - \bar c) \;\le\; 100 \;-\; \bar m^{2} \;-\; \bar c^{2}. +``` + +(Derived from the gradient of $g(m,c)=m^{2}+c^{2}-100$.) + +3. **Add Cut** + * Add cut & repeat. + +At any point in the solution process, the current LP solution looks like (for a given set of cuts indexed by $j$): +```math +\begin{aligned} +\min_{m,c} \quad & 200\,m + 80\,c \\ +\text{s.t.}\quad & 0 \le m \le 12, \\ + & 0 \le c \le 12, \\ + & 2\,\bar m_j\,(m - \bar m_j) + 2\,\bar c_j\,(c - \bar c_j) \le 100 - \bar m_j^{2} - \bar c_j^{2}, \quad j=1,\dots,J. +\end{aligned} +``` + +Your tasks + +1. **Implement** the iterative procedure in JuMP + HiGHS: + - master model with variables `m`, `c`; + - a loop that solves, checks feasibility, and adds cuts when needed. + +2. **Record** and (optionally) plot the objective value after each iteration. + +""" + +# ╔═╡ 7d855c60-41dd-40af-b00a-60b3e779ad13 +question_box(md"How can you check your answer?") + +# ╔═╡ ae263aac-1668-4f18-8104-ac25953a4503 +question_box(md""" +### How do you solve two-stage problems? + +How can you use cutting planes to solve a linear Two--Stage problem + +```math +\begin{aligned} +\min_{x}\; & c^{\top}x + \;+\; + \mathbb{E}_{\xi}\!\Bigl[ + \min_{y(\xi)} \; q^{\top}y(\xi) + \Bigr] \\[6pt] +\text{s.t.}\; & A\,x = b, \quad x \ge 0 + \qquad\qquad\text{(first-stage decisions)} +\end{aligned} +``` + +For each realisation $\xi = \bigl(h(\xi),\,T(\xi)\bigr)$ of the uncertain data, +the **second-stage (recourse) problem** is + +```math +\begin{aligned} +\min_{y(\xi)}\; & q^{\top}y(\xi) \\[4pt] +\text{s.t.}\; & T(\xi)\,y(\xi) = h(\xi) - W\,x, \\ + & y(\xi) \ge 0 + \qquad\qquad\text{(second-stage decisions)} +\end{aligned} +``` +""") + # ╔═╡ 808c505d-e10d-42e3-9fb1-9c6f384b2c3c md""" --- ## 2. MILP – 0‑1 Knapsack -You have a backpack that can carry at most **10 kg**. +You have a backpack that can carry at most **10 kg**. There are three items: | Item | Value | Weight | @@ -119,7 +338,7 @@ There are three items: | 3 | 5 | 2 | ### 2.1  Your tasks -1. Write the mathematical model with **binary** decision variables \(x_i \in \{0,1\}\). +1. Write the mathematical model with **binary** decision variables $x_i \in \{0,1\}$ that maximizes the acumulated value of the items in your bag. No partial items allowed. 2. Complete the JuMP model and solve it. 3. Pass the tests. """ @@ -130,64 +349,284 @@ begin # Replace the contents of this cell with your own model. model_milp = Model(HiGHS.Optimizer) -# Example: -# @variable(model_milp, x[1:3], Bin) +# Variables should be a vector named x (used for testing) +# @variable(model_milp, x[1:3] ... # --- YOUR CODE HERE --- # optimize!(model_milp) + +# Let's look at our model +println(model_milp) end # ╔═╡ 01367096-3971-4e79-ace2-83600672fbde begin -# === Quick check === -@testset "MILP check" begin - @test termination_status(model_milp) == MOI.OPTIMAL - @test isapprox(objective_value(model_milp), 22.0; atol = 1e-3) + ground_truth_2 = (x = [1.0, 1.0, 1.0], obj = 22.0) + + ans2 = missing + try + ans2 = ( + x = haskey(model_milp, :x) ? JuMP.value.(model_milp[:x]) : missing, + obj = objective_value(model_milp), + ) + catch + ans2 = missing + end + + good = !ismissing(ans2) && + all(isapprox.(ans2.x, ground_truth_2.x; atol=1e-3)) && + isapprox(ans2.obj, ground_truth_2.obj; atol=1e-3) + + if ismissing(ans2) + still_missing() + elseif good + correct() + else + keep_working() + end end + +# ╔═╡ 38b3a8f3-35ae-46da-91ce-0e4ba27ae098 +question_box(md"Is the answer the same if we allow partial products?") + +# ╔═╡ bca712e4-3f1c-467e-9209-e535aed5ab0a +md""" +### 2.2  Sudoku +(Credit to the similar JuMP Tutorial) + +Now let us solve the classic sudoku problem! + +1. Write the mathematical model with **binary** decision variables $x_i \in \{0,1\}$ that solves the following sudoku. +2. Complete the JuMP model and solve it. +3. Pass the tests. +""" + +# ╔═╡ 3997d993-0a31-435e-86cd-50242746c305 +@htl """ + +""" + +# ╔═╡ 3f56ec63-1fa6-403c-8d2a-1990382b97ae +begin +# === Your MILP model goes below === +# Replace the contents of this cell with your own model. +sudoku = Model(HiGHS.Optimizer) + +# Variables should be a vector named x_s (used for testing) +@variable(sudoku, x_s[i = 1:9, j = 1:9, k = 1:9], Bin); + +# --- YOUR CODE HERE --- + +# optimize!(sudoku) + +# Let's look at the stats of our model +sudoku +end + +# ╔═╡ 0e8ed625-df85-4bd2-8b16-b475a72df566 +begin + ground_truth_s = (x_ss = [[ 5 3 4 6 7 8 9 1 2]; + [6 7 2 1 9 5 3 4 8]; + [1 9 8 3 4 2 5 6 7]; + [8 5 9 7 6 1 4 2 3]; + [4 2 6 8 5 3 7 9 1]; + [7 1 3 9 2 4 8 5 6]; + [9 6 1 5 3 7 2 8 4]; + [2 8 7 4 1 9 6 3 5]; + [3 4 5 2 8 6 1 7 9]]) + + anss = missing + try + anss = ( + x_ss = haskey(sudoku, :x_s) ? JuMP.value.(sudoku[:x_s]) : missing, + ) + catch + anss = missing + end + + goods = !ismissing(anss) && + all(isapprox.(anss.x_ss, ground_truth_s.x_ss; atol=1e-3)) + + if ismissing(anss) + still_missing() + elseif goods + correct() + else + keep_working() + end end +# ╔═╡ fa5785a1-7274-4524-9e54-895d46e83861 +md""" +### 2.3 Other Modeling Tricks + +Below are some classic “tricks of the trade’’ for turning non-linear or logical +requirements into **mixed-integer linear programming (MILP)** form. +Throughout, $x\in\mathbb R^n$ are continuous decision variables and +$z\in\{0,1\}$ are binary (0–1) variables. + +--- + +#### Big-$M$ linearisation of conditional constraints +Suppose a constraint should apply **only if** a binary variable is 1: + +```math +z = 1 \;\Longrightarrow\; a^\top x \le b. +``` + +Introduce a sufficiently large constant $M>0$ and write + +```math +a^\top x \;\le\; b + M\,(1-z). +``` + +* If $z=1$ the right-hand side is $b$, so the original constraint is enforced. +* If $z=0$ the bound is relaxed by $M$ and becomes non-binding. + +> **Caveat:** pick $M$ as tight as possible. Over-large $M$ values hurt the LP +> relaxation and may cause numerical instability. + +--- + +#### Indicator constraints (a safer alternative) +Modern solvers allow *indicator* constraints that internally handle the +implication without an explicit $M$: + +```math +z = 1 \;\Longrightarrow\; a^\top x \le b. +``` + +In JuMP: + +```julia +@constraint(model, z --> a' * x <= b) +``` + +--- + +""" + # ╔═╡ 5e3444d0-8333-4f51-9146-d3d9625fe2e9 md""" --- -## 3. Non‑linear program – Rosenbrock valley +## 3. Non‑linear program – Modified Rosenbrock valley Non‑linear models dominate **optimal control** because discretising the differential equations that describe a physical system almost always yields a **non‑linear program (NLP)**. -A classic (and benign) test problem is the **Rosenbrock** function +Let's modify the classic (and benign) **Rosenbrock** function -\[ -\min_{x,\,y\in\mathbb R} \; f(x,y)= (1-x)^{2} + 100\,(y - x^{2})^{2}. -\] +```math +\begin{aligned} +\min_{x,\,y} \quad & f(x,y) \;=\; -(1-x)^{2} + 100\,(y - x^{2})^{2} + 10000x \\ +\text{s.t.}\quad & -10 \le x \le 10,\\ + & -60 \le y \le 60. +\end{aligned} +``` -It has a single global optimum at \((x^{\star},y^{\star}) = (1,1)\) with \(f^{\star}=0\). +It has a single global minimum within the feasible reagion defined by the box constraints on $x$ and $y$. ### 3.1  Your tasks -1. Build and solve the model with **Ipopt**. +1. Build a model to find this minimum and solve it with **Ipopt**. 2. Inspect the solution and objective. -3. Check your work below. +3. Pass the tests """ +# ╔═╡ 0e190de3-da60-41e9-9da5-5a0c7fefd1d7 +f(x, y) = -(1-x)^2 + 100 * (y-x^2)^2 + 10000*x + +# ╔═╡ cac18d70-b354-48c7-9f37-31ee0c585675 +begin +# 1. Domain grids +xs = range(-10, 10; length = 100) # 100 x-points +ys = range(-60, 60; length = 100) # 100 y-points + +# 3. Surface heights — matrix z[y, x] +zs = [f(x, y) for y in ys, x in xs] + +# 4. Create a 1×3 layout +fig = Figure(size = (1500, 500)) + +# 4a. 3-D surface +ax3d = Axis3(fig[1, 1]; xlabel = "y", ylabel = "x", zlabel = "f(x,y)", + title = "Surface") +surface!(ax3d, ys, xs, zs; colormap = :viridis) +ax3d.azimuth[] = deg2rad(-10) # ≈ camera = (-10°, 30°) +ax3d.elevation[] = deg2rad( 30) + +# 4b. Slice at y = +y0=40 +ax_x = Axis(fig[1, 2]; xlabel = "x", ylabel = "f(x,$(y0))", title = "y = $(y0) slice") +lines!(ax_x, xs, f.(xs, y0)) + +# 4c. Slice at x = +x_0 = 0 +ax_y = Axis(fig[1, 3]; xlabel = "y", ylabel = "f($(x_0),y)", title = "x = $(x_0) slice") +lines!(ax_y, ys, f.(x_0, ys)) + +fig +end + # ╔═╡ 00728de8-3c36-48c7-8520-4c9f408a7c5f begin # === Your NLP model goes below === # Replace the contents of this cell with your own model. model_nlp = Model(Ipopt.Optimizer) +# Required named variables +@variable(model_nlp, x) +@variable(model_nlp, y) + # --- YOUR CODE HERE --- # optimize!(model_nlp) + +println(model_nlp) end -# ╔═╡ 254b9a87-17f9-4fea-8b28-0e3873b58fe2 +# ╔═╡ 45541136-695a-4260-82c1-66d38ec44dcc +md""" +### 3.2 Intersecting-ellipse constraints +Find the minimum of the same modified Rosenbrock objective, **but now the feasible +region is the intersection of three ellipses** defined only through their focal points +and the constant sum-of-distances $2a$: + +| Ellipse | Focal points $(p_1,\,p_2)$ | Required sum of distances $2a$ | +|---------|---------------------------|--------------------------------| +| $E_1$ | $(-4,\,0),\;(4,\,0)$ | $12$ | +| $E_2$ | $(0,\,-5),\;(0,\,5)$ | $14$ | +| $E_3$ | $(-3,\,-3),\;(3,\,3)$ | $12$ | + +Recall that a point $(x,y)$ lies **inside** an ellipse if the **sum of its Euclidean +distances to the two foci is _no greater_ than $2a$**. +Formulate these three nonlinear constraints and use **Ipopt** to locate the optimal +$(x^*,y^*)$ and corresponding objective value. + +1. Implement the model with the three ellipse constraints. +2. Solve it and report the optimal point and objective. +3. Verify that the solver stopped at a feasible point (all three distance sums $\le 2a$). +""" + +# ╔═╡ b107bcd7-60ca-4f09-aa42-f8335e13233e begin -# === Quick check === -@testset "NLP check" begin - @test termination_status(model_nlp) == MOI.LOCALLY_SOLVED || termination_status(model_nlp) == MOI.OPTIMAL - @test isapprox(objective_value(model_nlp), 0.0; atol = 1e-6) -end + # ── NLP model ──────────────────────────────────────────────────────────────── + model_nlp2 = Model(Ipopt.Optimizer) + + # Decision variables with box bounds + @variable(model_nlp2, x2) + @variable(model_nlp2, y2) + + # --- YOUR CODE HERE --- + + # Solve + # optimize!(model_nlp2) + + # Quick report + println(model_nlp2) end # ╔═╡ 147fe732-fe65-4226-af43-956b33a75bff @@ -215,18 +654,160 @@ Getting comfortable with modelling and debugging small nonlinear examples like R when you step up to thousands of variables in real control problems! """ +# ╔═╡ 87ffc247-3769-4002-a584-c687fd813125 +begin + hidden_answers = Dict( + :lp => (w = 20.0, g = 15.0, obj = 135.0), + :milp => (x = [1,1,1], obj = 22.0), + :nlp => (x = 1.0, y = 1.0, obj = 0.0), + ) + safeval(model, sym) = haskey(model, sym) ? JuMP.value(model[sym]) : missing + md"" +end + +# ╔═╡ 6fb672d0-5a18-4ccc-b7b3-184839c2401b +begin + # ground truth + ground_truth = (w = 20.0, g = 15.0, obj = 135.0) + + # student answer + ans = missing + try + ans = ( + w = safeval(model_lp, :w), + g = safeval(model_lp, :g), + obj = objective_value(model_lp), + ) + catch # objective_value will throw if model_lp not ready + ans = missing + end + + # Decide which badge to show + if ismissing(ans) # nothing yet + still_missing() + elseif isapprox(ans.w, ground_truth.w; atol=1e-3) && + isapprox(ans.g, ground_truth.g; atol=1e-3) && + isapprox(ans.obj, ground_truth.obj; atol=1e-3) + correct() + else + keep_working() + end +end + +# ╔═╡ 20aef3e9-47b5-4f60-9726-7db77f7c3e47 +begin + # student answer + ansd = missing + try + ansd = safeval(model_lp2, :x_d) + catch # objective_value will throw if model_lp not ready + ansd = missing + end + + # Decide which badge to show + if ismissing(ansd) # nothing yet + still_missing() + elseif x == 25.0 + correct() + else + keep_working() + end +end + +# ╔═╡ 254b9a87-17f9-4fea-8b28-0e3873b58fe2 +begin + ground_truth_3 = (x = -7.946795, y = 60.00, obj = -78554.7682) + + ans3=missing + try + ans3 = ( + x = safeval(model_nlp, :x), + y = safeval(model_nlp, :y), + obj = objective_value(model_nlp), + ) + catch + ans3 = missing + end + + good_2 = !ismissing(ans3) && + isapprox(ans3.x, ground_truth_3.x; atol=1e-3) && + isapprox(ans3.y, ground_truth_3.y; atol=1e-3) && + isapprox(ans3.obj, ground_truth_3.obj; atol=1e-3) + + if ismissing(ans3) + still_missing() + elseif good_2 + correct() + else + keep_working() + end +end + +# ╔═╡ 1fcaedb8-34d0-4faf-9052-fc074d2edda3 +begin + ground_truth_32 = (x = -3.121657, y = 2.875823, obj = -26515.3545) + + ans4=missing + try + ans4 = ( + x = safeval(model_nlp2, :x2), + y = safeval(model_nlp2, :y2), + obj = objective_value(model_nlp2), + ) + catch + ans4 = missing + end + + good_3 = !ismissing(ans4) && + isapprox(ans4.x, ground_truth_32.x; atol=1e-3) && + isapprox(ans4.y, ground_truth_32.y; atol=1e-3) && + isapprox(ans4.obj, ground_truth_32.obj; atol=1e-3) + + if ismissing(ans4) + still_missing() + elseif good_3 + correct() + else + keep_working() + end +end + # ╔═╡ Cell order: # ╟─881eed45-e7f0-4785-bde8-530e378d7050 # ╟─9f5675a3-07df-4fb1-b683-4c5fd2a85002 +# ╟─b879c4bf-8292-45a0-87d0-b1f0fa456585 +# ╟─eeceb82e-abfb-4502-bcfb-6c9f76a0879d # ╟─0df8b65a-0527-4545-bf11-00e9912bced0 # ╠═9ce52307-bc22-4f66-a4af-a4e4ac382212 # ╟─6f67ca7c-1391-4cb9-b692-cd818e037587 # ╠═49042d6c-cf78-46d3-bfee-a8fd7ddf3aa0 -# ╠═6fb672d0-5a18-4ccc-b7b3-184839c2401b -# ╠═808c505d-e10d-42e3-9fb1-9c6f384b2c3c +# ╟─1d3edbdd-7747-4651-b650-c6b9bf87b460 +# ╟─6fb672d0-5a18-4ccc-b7b3-184839c2401b +# ╟─248b398a-0cf5-4c2b-8752-7b9cc4e765d6 +# ╟─245eb671-84e1-447b-8045-e9eb04966d80 +# ╟─6a823649-04fa-4322-a028-2fb29dffb08b +# ╠═c369ab46-b416-4c12-83fe-65040a0c47c8 +# ╟─20aef3e9-47b5-4f60-9726-7db77f7c3e47 +# ╟─b13f9775-68c2-4646-9b67-c69ee23a4ea0 +# ╟─ea3ea95a-58cb-4d0d-a167-aa68b8bc2645 +# ╟─7d855c60-41dd-40af-b00a-60b3e779ad13 +# ╟─ae263aac-1668-4f18-8104-ac25953a4503 +# ╟─808c505d-e10d-42e3-9fb1-9c6f384b2c3c # ╠═39617561-bbbf-4ef6-91e2-358dfe76581c -# ╠═01367096-3971-4e79-ace2-83600672fbde -# ╠═5e3444d0-8333-4f51-9146-d3d9625fe2e9 +# ╟─01367096-3971-4e79-ace2-83600672fbde +# ╟─38b3a8f3-35ae-46da-91ce-0e4ba27ae098 +# ╟─bca712e4-3f1c-467e-9209-e535aed5ab0a +# ╟─3997d993-0a31-435e-86cd-50242746c305 +# ╠═3f56ec63-1fa6-403c-8d2a-1990382b97ae +# ╟─0e8ed625-df85-4bd2-8b16-b475a72df566 +# ╟─fa5785a1-7274-4524-9e54-895d46e83861 +# ╟─5e3444d0-8333-4f51-9146-d3d9625fe2e9 +# ╠═0e190de3-da60-41e9-9da5-5a0c7fefd1d7 +# ╟─cac18d70-b354-48c7-9f37-31ee0c585675 # ╠═00728de8-3c36-48c7-8520-4c9f408a7c5f -# ╠═254b9a87-17f9-4fea-8b28-0e3873b58fe2 +# ╟─254b9a87-17f9-4fea-8b28-0e3873b58fe2 +# ╟─45541136-695a-4260-82c1-66d38ec44dcc +# ╠═b107bcd7-60ca-4f09-aa42-f8335e13233e +# ╟─1fcaedb8-34d0-4faf-9052-fc074d2edda3 # ╟─147fe732-fe65-4226-af43-956b33a75bff +# ╟─87ffc247-3769-4002-a584-c687fd813125 diff --git a/class01/class01_intro.jl b/class01/class01_intro.jl index 7e3f76f..781a102 100644 --- a/class01/class01_intro.jl +++ b/class01/class01_intro.jl @@ -49,6 +49,26 @@ md" | Date | : | 28 of July, 2025 | " +# ╔═╡ ced1b968-3ba6-4e58-9bcd-bbc6bee2b93c +md"#### Reference Material" + +# ╔═╡ 97994ed8-5606-46ef-bd30-c5343c1d99cf +begin + MarkdownLiteral.@markdown( +""" + +[^cmu]: Zachary Manchester et al. Optimal Control and Reinforcement Learning at Carnegie Mellon University - [CMU 16-745]("https://optimalcontrol.ri.cmu.edu/") + +[^OptProx]: Van Hentenryck, P., 2024. Fusing Artificial Intelligence and Optimization with Trustworthy Optimization Proxies. Collections, 57(02). + +[^ArmManip]: Guechi, E.H., Bouzoualegh, S., Zennir, Y. and Blažič, S., 2018. MPC control and LQ optimal control of a two-link robot arm: A comparative study. Machines, 6(3), p.37. + +[^ZachMIT]: Zachary Manchester talk at MIT - [MIT Robotics - Zac Manchester - Composable Optimization for Robotic Motion Planning and Control]("https://www.youtube.com/watch?v=eSleutHuc0w&ab_channel=MITRobotics"). + +""" +) +end + # ╔═╡ 1f774f46-d57d-4668-8204-dc83d50d8c94 md"# Intro - Optimal Control and Learning @@ -146,14 +166,200 @@ in stochastic programming. # ╔═╡ 5d7a4408-21ff-41ec-b004-4b0a9f04bb4f question_box(md"Can you name a few ways to try and/or solve this problem?") +# ╔═╡ 7e487ebc-8327-4f3e-a8ca-1e07fb39991a +md""" +### Solution Methods + +There are a few ways to solve these problems: + +```math +(\mathbf{x}_{t-1}, w_t)\xrightarrow[\pi_t^{*}(\mathbf{x}_{t-1}, w_t)]{ +\begin{align} + &\min_{\mathbf{x}_t, \mathbf{u}_t} \quad \! \! c(\mathbf{x}_t, \mathbf{u}_t) + \mathbf{E}_{t+1}[V_{t+1}(\mathbf{x}_t, w_{t+1})] \\ + & \text{ s.t. } \quad\mathbf{x}_t = f(\mathbf{x}_{t-1}, w_t, \mathbf{u}_t) \nonumber \\ + & \quad \quad \quad \! \! h(\mathbf{x}_t, \mathbf{u}_t) \geq 0. \nonumber +\end{align} +} (\mathbf{x}_t^{*}, \mathbf{u}_t^{*}) +``` + +**Exact Methods:** + - Deterministic Equivalent: Explicitly model all decisions of all possible scenarios. (Good Luck!) + - Stochastic Dual Dynamic Programming, Progressive Hedging, ... (Hard but doable for some class of problems.) + +**Approximate Methods**: + - Approximate Dynamic Programming, (model-free and model-based)Reinforcement Learning, Two-Stage Decision Rules, ... + - **Optimization Proxies**: + +```math +\theta^{\star} +\;=\; +\operatorname*{arg\,min}_{\theta \in \Theta} +\; +\mathbb{E}\Bigl[\bigl\|\,\pi_t^{\ast}-\pi_t(\,\cdot\,;\theta)\bigr\|_{\mathcal F}\Bigr], +``` + +""" + +# ╔═╡ bd623016-24ce-4c10-acb3-b2b80d4facc8 +md"[^OptProx]" + +# ╔═╡ 2d211386-675a-4223-b4ca-124edd375958 +@htl """ + + + +""" + +# ╔═╡ 45275d44-e268-43cb-8156-feecd916a6da +@htl """ +
+ + +

LearningToOptimize Organization

+ +

+ LearningToOptimize (L2O) is a collection of open-source tools + focused on the emerging paradigm of amortized optimization—using machine-learning + methods to accelerate traditional constrained-optimization solvers. + L2O is a work-in-progress; existing functionality is considered experimental and may + change. +

+ + +

Open-Source Repositories

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ LearningToOptimize.jl + + Flagship Julia package that wraps data generation, training loops and evaluation + utilities for fitting surrogate models to parametric optimization problems. +
+ DecisionRules.jl + + Build decision rules for multistage stochastic programs, as proposed in + Efficiently + Training Deep-Learning Parametric Policies using Lagrangian Duality. +
+ L2OALM.jl + + Implementation of the primal-dual learning method ALM, + introduced in + + Self-Supervised Primal-Dual Learning for Constrained Optimization. +
+ L2ODLL.jl + + Implementation of the dual learning method DLL, + proposed in + + Dual Lagrangian Learning for Conic Optimization. +
+ L2ODC3.jl + + Implementation of the primal learning method DC3, as described in + + DC3: A Learning Method for Optimization with Hard Constraints. +
+ BatchNLPKernels.jl + + GPU kernels that evaluate objectives, Jacobians and Hessians for + batches of + ExaModels, + useful when defining loss functions for large-batch ML predictions. +
+ BatchConeKernels.jl + + GPU kernels for batched cone operations (projections, distances, etc.), + enabling advanced architectures such as repair layers. +
+ LearningToControlClass + + Course repository for Special Topics on Optimal Control & Learning + (Fall 2025, Georgia Tech). +
+ + +

Open Datasets and Weights

+ +

+ The + + LearningToOptimize 🤗 Hugging Face organization + hosts datasets and pre-trained weights that can be used with L2O packages. +

+ +
+""" + # ╔═╡ c08f511e-b91d-4d17-a286-96469c31568a md"## Example: Robotic Arm Manipulation" # ╔═╡ b3129bcb-c24a-4faa-a5cf-f69ce518ea87 -load(joinpath(class_dir, "nlp_robot_arm.png")) +begin + load(joinpath(class_dir, "nlp_robot_arm.png")) +end # ╔═╡ c1f43c8d-0616-4572-bb48-dbb71e40adda md""" +[^ArmManip] + The tip of the second link is computed using the direct geometric model: ```math @@ -169,6 +375,7 @@ y = L_{1}\,\cos\theta_{1} \;+\; L_{2}\,\cos\!\bigl(\theta_{1}+\theta_{2}\bigr). # ╔═╡ 57d896ca-221a-4cfc-b37a-be9898fac923 begin md""" + **State** ```math \mathbf{x}_t=\begin{bmatrix}\theta_{1,t}&\theta_{2,t}&\dot\theta_{1,t}&\dot\theta_{2,t}\end{bmatrix}^{\!\top} @@ -215,9 +422,6 @@ h(\mathbf{x}_t,\mathbf{u}_t)\ge 0\;:\; """ end -# ╔═╡ e2d3d160-d3b6-41f2-a8bc-2878ba71e78c - - # ╔═╡ 52005382-177b-4a11-a914-49a5ffc412a3 section_outline(md"A Crash Course:",md" (Continuous-Time) Dynamics ") @@ -463,11 +667,22 @@ Foldable(md"All mechanical systems can be written this way. Why?", md""" Manipulator Dynamics Equations are a way of rewriting the Euler--Lagrange equations. + +> The equations were discovered in the 1750s by Swiss mathematician Leonhard Euler and Italian mathematician Joseph-Louis Lagrange. + +""") + +# ╔═╡ 5a691d10-44f7-4d44-a2c9-a7d4d720b7cc +begin +md""" #### 🚀 Detour: The Principle of Least Action 🚀 -In the calculus of variations and classical mechanics, the Euler–Lagrange equations are a system of second-order ordinary differential equations whose solutions are stationary points of the given action functional. +In the calculus of variations and classical mechanics, the Euler–Lagrange equations are a system of second-order ordinary differential equations whose solutions are stationary points of the given action functional: -> The equations were discovered in the 1750s by Swiss mathematician Leonhard Euler and Italian mathematician Joseph-Louis Lagrange. +```math +\mathcal{S}[q(\cdot)] \;=\; +\int_{t_0}^{t_f} L\!\bigl(q(t),\; \dot q(t)\bigr)\,dt, +``` In classical mechanics: @@ -475,6 +690,14 @@ In classical mechanics: L = \underbrace{\frac{1}{2} v^{\top}M(q)v}_{\text{Kinematic Energy}} - \underbrace{U(q)}_{\text{Potential Energy}} ``` +""" +end + +# ╔═╡ f3d155c6-5384-481a-8373-582e753ea8d6 +question_box(md"What can you say about $M(q)$? When do we have a problem inverting it?") + +# ╔═╡ ee5c5e2e-e9f1-4f94-95c9-21d506281ae1 +md""" A curve ($q^\star(t)$) is physically realised iff it is a stationary point of ($\mathcal{S}$) : @@ -487,16 +710,92 @@ point of ($\mathcal{S}$) : M(q)\,\ddot q + C(q,\dot q)\,\dot q + \nabla U(q)=0 . ``` -""") - -# ╔═╡ f3d155c6-5384-481a-8373-582e753ea8d6 -question_box(md"What can you say about $M(q)$? When do we have a problem inverting it?") +""" # ╔═╡ b9aeab8a-f8ea-4310-8568-5d6bda0bb4d3 question_box(md"Can you derive the stationary condition?") -# ╔═╡ e1dc6ecf-4e62-415a-a620-0731953c5ab4 +# ╔═╡ 30a013a8-c02e-4816-af0d-9280473c916b +md""" +In most cases: +```math +q^{*} \in \arg \min_{q} +\int_{t_0}^{t_f} L\!\bigl(q(t),\; \dot q(t)\bigr)\,dt, +``` + +Now, suppose the configuration must satisfy a *gap function* +$\phi(q)\ge 0$ (e.g. **contact with the ground**, obstacle avoidance, joint +limits). +The variational problem becomes +```math +q^{*} \;\in\; +\arg\!\min_{q(\cdot)} +\int_{t_0}^{t_f} L\!\bigl(q(t),\dot q(t)\bigr)\,dt +\quad\text{s.t.}\quad +\phi\!\bigl(q(t)\bigr)\;\ge 0 \;\;\;\forall\,t. +``` + +Let $(t_k = t_0 + k,\Delta t)$ with $(k=0,\dots,N)$ and +$(q_k \approx q(t_k))$. +Using the midpoint rule we approximate the action by + +```math +S_N(q_{0:N}) +\;=\; +\sum_{k=0}^{N-1} +L\!\Bigl( + \tfrac12\bigl(q_k+q_{k+1}\bigr),\; + \tfrac{q_{k+1}-q_k}{\Delta t} +\Bigr)\,\Delta t, +``` + +and obtain the finite‐dimensional problem + +```math +\begin{aligned} +\min_{q_1,\dots,q_{N}} +& \; S_N(q_{0:N}) \\[4pt] +\text{s.t.}\;&\; + \phi(q_{k+1}) \;\ge 0, + \qquad k = 0,\dots,N-1. +\end{aligned} +``` + +""" + +# ╔═╡ 2cc57795-717a-46f0-9bb5-67b601a766de +begin + gif_url = "https://raw.githubusercontent.com/dojo-sim/Dojo.jl/main/docs/src/assets/animations/atlas_drop.gif" + still_url = "https://encrypted-tbn2.gstatic.com/images?q=tbn:ANd9GcQkrtL7TCGzNxFlXIqYHW_cFP9pfLscwd7vLSH09nfRFEQCqX_J" + md"" +end + +# ╔═╡ 59f6167d-796c-4844-89c0-c796fb59aa2e +Columns(md"[^ZachMIT]", md"▶/⏸$(@bind playing CheckBox(default=false))") + +# ╔═╡ 58c2e1f2-819d-40fc-8e92-03a1a3019a3d +Columns(md""" +$(load(joinpath(class_dir, "rocket_physics.png"))) + +#### Dojo.jl + +A differentiable physics engine for robotics that simulates systems using optimization. + +- [ArXiv preprint](https://arxiv.org/abs/2203.00806) +- [GitHub](https://github.com/dojo-sim/Dojo.jl) + +""" +, +@htl """ + +""" +) + +# ╔═╡ 70690e72-c31e-4c91-b211-35c74d1d9973 +warning_box(md"But in general we need a *ReFeynman* of the these equations!") # ╔═╡ 5f35a169-887f-477f-b010-167627f7ce4c md"## Linear Systems @@ -1142,32 +1441,26 @@ end # ╔═╡ de4807ca-4e17-4020-9810-5f7c0fcae9a3 question_box(md"### Why most simulators use Backward--Euler?") -# ╔═╡ 97994ed8-5606-46ef-bd30-c5343c1d99cf -begin - MarkdownLiteral.@markdown( -""" - -[^cmu]: [CMU Course on Optimal Control]("https://optimalcontrol.ri.cmu.edu/") - -""" -) -end - # ╔═╡ Cell order: # ╟─13b12c00-6d6e-11f0-3780-a16e73360478 # ╟─ec473e69-d5ec-4d6a-b868-b89dadb85705 # ╟─8d7a34ef-5a2d-41a8-ac55-39ab00d7e432 +# ╟─ced1b968-3ba6-4e58-9bcd-bbc6bee2b93c +# ╟─97994ed8-5606-46ef-bd30-c5343c1d99cf # ╟─1f774f46-d57d-4668-8204-dc83d50d8c94 # ╟─a0f71960-c97c-40d1-8f78-4b1860d2e0a2 # ╟─1d7092cd-0044-4d38-962a-ce3214c48c24 # ╟─60ba261a-f2eb-4b45-ad6d-b6042926ccab # ╟─15709f7b-943e-4190-8f40-0cfdb8772183 # ╟─5d7a4408-21ff-41ec-b004-4b0a9f04bb4f +# ╟─7e487ebc-8327-4f3e-a8ca-1e07fb39991a +# ╟─bd623016-24ce-4c10-acb3-b2b80d4facc8 +# ╟─2d211386-675a-4223-b4ca-124edd375958 +# ╟─45275d44-e268-43cb-8156-feecd916a6da # ╟─c08f511e-b91d-4d17-a286-96469c31568a # ╟─b3129bcb-c24a-4faa-a5cf-f69ce518ea87 # ╟─c1f43c8d-0616-4572-bb48-dbb71e40adda # ╟─57d896ca-221a-4cfc-b37a-be9898fac923 -# ╠═e2d3d160-d3b6-41f2-a8bc-2878ba71e78c # ╟─52005382-177b-4a11-a914-49a5ffc412a3 # ╟─8ea866a6-de0f-4812-8f59-2aebec709243 # ╟─2be161cd-2d4c-4778-adca-d45f8ab05f98 @@ -1185,9 +1478,15 @@ end # ╟─f10927fe-d392-4374-bad1-ab5ac85b8116 # ╟─b8b206ef-cdc5-4cc9-9b55-70d711ba2a9e # ╟─a09de9e4-7ecc-4d23-9135-384077f0c03f +# ╟─5a691d10-44f7-4d44-a2c9-a7d4d720b7cc # ╟─f3d155c6-5384-481a-8373-582e753ea8d6 +# ╟─ee5c5e2e-e9f1-4f94-95c9-21d506281ae1 # ╟─b9aeab8a-f8ea-4310-8568-5d6bda0bb4d3 -# ╠═e1dc6ecf-4e62-415a-a620-0731953c5ab4 +# ╟─30a013a8-c02e-4816-af0d-9280473c916b +# ╟─2cc57795-717a-46f0-9bb5-67b601a766de +# ╟─59f6167d-796c-4844-89c0-c796fb59aa2e +# ╟─58c2e1f2-819d-40fc-8e92-03a1a3019a3d +# ╟─70690e72-c31e-4c91-b211-35c74d1d9973 # ╟─5f35a169-887f-477f-b010-167627f7ce4c # ╟─e860d92b-cc8f-479b-a0fc-e5f7a11ae1fd # ╟─bb4bfa72-bf69-41f5-b017-7cbf31653bae @@ -1239,4 +1538,3 @@ end # ╟─86ce1303-e77c-4b93-a2ed-dc0c54a1f191 # ╠═b857efd5-dba1-4872-b133-59e80d7cd489 # ╟─de4807ca-4e17-4020-9810-5f7c0fcae9a3 -# ╟─97994ed8-5606-46ef-bd30-c5343c1d99cf diff --git a/class01/rocket_physics.png b/class01/rocket_physics.png new file mode 100644 index 0000000..ab1c163 Binary files /dev/null and b/class01/rocket_physics.png differ