Add Examples 2.0

docs/src/examples/Mean_Variance_Portfolio_Example.jl

@@ -0,0 +1,155 @@
+# # Thermal Generation Dispatch Sweep Example
+
+# Consider the Markowitz portfolio selection problem, which allocates weights $x \in \mathbb{R}^n$ to $n$ assets so as to maximize returns subject to a variance limit $v_{\max}$:
+# ```math
+# \max_{x} \quad \mu^\top x
+# \quad\text{s.t.}\quad
+# x^\top \Sigma x \;\le\; v_{\max}, \quad
+# \mathbf{1}^\top x = 1,\quad
+# x \succeq 0,
+# ```
+# where $\mu$ is the vector of expected returns, $\Sigma$ is the covariance matrix, and $x$ must sum to 1 (fully invest the budget). An efficient conic version of this problem casts the variance limit as a second order cone constraint:
+
+# ```math
+# \| \Sigma^{1/2} x \|_{2} \;\le\; \sigma_{\max}
+# ```
+# where $\Sigma^{1/2}$ is the Cholesky factorization of the covariance matrix and $\sigma_{\max}$ is the standard deviation limit.
+
+# Practitioners often care about an \emph{out-of-sample performance metric} $L(x)$ evaluated on test data or scenarios that differ from those used to form $\mu$ and $\Sigma$. To assess the impact of the risk profile in the performance evaluation, one can compute:
+# ```math
+# \frac{dL}{d\,\sigma_{\max}} \;=\;
+# \underbrace{\frac{\partial L}{\partial x}}_{\text{(1) decision impact}}\;
+# \cdot\;
+# \underbrace{\frac{\partial x^*}{\partial \sigma_{\max}}}_{\text{(2) from DiffOpt.jl}},
+# ```
+# where $x^*(\sigma_{\max})$ is the portfolio that solves the conic Markowitz problem under a given risk limit.
+
+# ## Define and solve the Mean-Variance Portfolio Problem for a range of risk limits
+
+# First, import the libraries.
+
+using Test
+using JuMP
+import DiffOpt
+using LinearAlgebra
+import SCS
+using Plots
+using Plots.Measures
+
+# Fixed data
+
+# Training data (in-sample)
+Σ = [
+    0.002 0.0005 0.001
+    0.0005 0.003 0.0002
+    0.001 0.0002 0.0025
+]
+μ_train = [0.05, 0.08, 0.12]
+
+# Test data (out-of-sample)
+μ_test = [0.02, -0.3, 0.1]             # simple forecast error example
+
+# Sweep over σ_max
+σ_grid = 0.002:0.002:0.06
+N = length(σ_grid)
+
+predicted_ret = zeros(N)                 # μ_train' * x*
+realised_ret = zeros(N)                 # μ_test'  * x*
+loss = zeros(N)                 # L(x*)
+dL_dσ = zeros(N)                 # ∂L/∂σ_max  from DiffOpt
+
+for (k, σ_val) in enumerate(σ_grid)
+
+    ## 1) differentiable conic model
+    model = DiffOpt.conic_diff_model(SCS.Optimizer)
+    set_silent(model)
+
+    ## 2) parameter σ_max
+    @variable(model, σ_max in Parameter(σ_val))
+
+    ## 3) portfolio weights
+    @variable(model, x[1:3] >= 0)
+    @constraint(model, sum(x) <= 1)
+
+    ## 4) objective: maximise expected return (training data)
+    @objective(model, Max, dot(μ_train, x))
+
+    ## 5) conic variance constraint  ||L*x|| <= σ_max
+    L_chol = cholesky(Symmetric(Σ)).L
+    @variable(model, t >= 0)
+    @constraint(model, [t; L_chol * x] in SecondOrderCone())
+    @constraint(model, t <= σ_max)
+
+    optimize!(model)
+
+    x_opt = value.(x)
+    println("Optimal portfolio weights: ", x_opt)
+
+    ## store performance numbers
+    predicted_ret[k] = dot(μ_train, x_opt)
+    realised_ret[k] = dot(μ_test, x_opt)
+
+    ## -------- reverse differentiation wrt σ_max --------
+    DiffOpt.empty_input_sensitivities!(model)
+    ## ∂L/∂x   (adjoint)  =  -μ_test
+    DiffOpt.set_reverse_variable.(model, x, μ_test)
+    DiffOpt.reverse_differentiate!(model)
+    dL_dσ[k] = DiffOpt.get_reverse_parameter(model, σ_max)
+end
+
+# ## Results with Plot graphs
+
+default(;
+    size = (1150, 350),
+    legendfontsize = 8,
+    guidefontsize = 9,
+    tickfontsize = 7,
+)
+
+# (a) predicted vs realised return
+plt_ret = plot(
+    σ_grid,
+    realised_ret;
+    lw = 2,
+    label = "Realised (test)",
+    xlabel = "σ_max (risk limit)",
+    ylabel = "Return",
+    title = "Return vs risk limit",
+    legend = :bottomright,
+);
+plot!(
+    plt_ret,
+    σ_grid,
+    predicted_ret;
+    lw = 2,
+    ls = :dash,
+    label = "Predicted (train)",
+);
+
+# (b) out-of-sample loss and its gradient
+plt_loss = plot(
+    σ_grid,
+    dL_dσ;
+    xlabel = "σ_max (risk limit)",
+    ylabel = "∂L/∂σ_max",
+    title = "Return Gradient",
+    legend = false,
+);
+
+plot_all = plot(
+    plt_ret,
+    plt_loss;
+    layout = (1, 2),
+    left_margin = 5Plots.Measures.mm,
+    bottom_margin = 5Plots.Measures.mm,
+)
+
+# Impact of the risk limit $\sigma_{\max}$ on Markowitz
+# portfolios.  **Left:** predicted in-sample return versus
+# realized out-of-sample return.  **Right:** the
+# out-of-sample loss $L(x)$ together with the absolute gradient
+# $|\partial L/\partial\sigma_{\max}|$ obtained from
+# `DiffOpt.jl`.  The gradient tells the practitioner which
+# way—and how aggressively—to adjust $\sigma_{\max}$ to reduce
+# forecast error; its value is computed in one reverse-mode call
+# without re-solving the optimization for perturbed risk limits.

docs/src/examples/Planar_Arm_Example.jl

@@ -0,0 +1,176 @@
+# # Planar Arm Example
+
+# Inverse Kinematics (IK) computes joint angles that place a robot’s end-effector at a desired target $(x_t,y_t)$. For a 2-link planar arm with joint angles $\theta_1,\theta_2$, the end-effector position is:
+# ```math
+# f(\theta_1,\theta_2) = \bigl(\ell_1\cos(\theta_1) + \ell_2\cos(\theta_1+\theta_2),\,\,
+# \ell_1\sin(\theta_1) + \ell_2\sin(\theta_1+\theta_2)\bigr).
+# ```
+# We can solve an NLP:
+# ```math
+# \min_{\theta_1,\theta_2} \;\; (\theta_1^2 + \theta_2^2),
+# \quad\text{s.t.}\quad f(\theta_1,\theta_2) = (x_t,y_t).
+# ```
+# Treat $(x_t,y_t)$ as parameters. Once solved, we differentiate w.r.t. $(x_t,y_t)$ to find how small changes in the target location alter the optimal angles - the *differential kinematics*.
+
+# ## Define and solve the 2-link planar arm problem and build sensitivity map of joint angles to target
+
+# First, import the libraries.
+
+using Test
+using JuMP
+import DiffOpt
+using LinearAlgebra
+using Statistics
+import Ipopt
+using Plots
+using Plots.Measures
+
+# Fixed data
+
+# Arm geometry
+l1 = 1.0;
+l2 = 1.0;
+reach = l1 + l2          # 2.0
+tol = 1e-6               # numerical tolerance for feasibility
+# Sampling grid in workspace
+grid_res = 25 # grid resolution low for documentation compilation requirements
+xs = range(-reach, reach; length = grid_res)
+ys = range(-reach, reach; length = grid_res)
+
+heat = fill(NaN, grid_res, grid_res)   # store ‖J_inv‖₂
+feas = fill(0.0, grid_res, grid_res)  # feasibility mask
+θ1mat = similar(heat)
+
+function ik_angles(x, y; l1 = 1.0, l2 = 1.0, elbow_up = true)
+    c2 = (x^2 + y^2 - l1^2 - l2^2) / (2 * l1 * l2)
+    θ2 = elbow_up ? acos(clamp(c2, -1, 1)) : -acos(clamp(c2, -1, 1))
+    k1 = l1 + l2 * cos(θ2)
+    k2 = l2 * sin(θ2)
+    θ1 = atan(y, x) - atan(k2, k1)
+    return θ1, θ2
+end
+
+for (i, x_t) in enumerate(xs), (j, y_t) in enumerate(ys)
+    global θ1mat, heat
+    ## skip points outside the circular reach
+    norm([x_t, y_t]) > reach && continue
+
+    ## ---------- build differentiable NLP ----------
+    model = DiffOpt.nonlinear_diff_model(Ipopt.Optimizer)
+    set_silent(model)
+
+    @variable(model, xt in Parameter(x_t))
+    @variable(model, yt in Parameter(y_t))
+    @variable(model, θ1)
+    @variable(model, θ2)
+
+    @objective(model, Min, θ1^2 + θ2^2)
+    @constraint(model, l1 * cos(θ1) + l2 * cos(θ1 + θ2) == xt)
+    @constraint(model, l1 * sin(θ1) + l2 * sin(θ1 + θ2) == yt)
+
+    ## --- supply analytic start values ---
+    θ1₀, θ2₀ = ik_angles(x_t, y_t; elbow_up = true)
+    set_start_value(θ1, θ1₀)
+    set_start_value(θ2, θ2₀)
+
+    optimize!(model)
+    println("Solving for target (", x_t, ", ", y_t, ")")
+    ## check for optimality
+    status = termination_status(model)
+    println("Status: ", status)
+
+    status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || continue
+
+    θ1̂ = value(θ1)
+    θ1mat[j, i] = θ1̂   # save pose
+
+    ## ---- forward diff wrt  xt  (∂θ/∂x) ----
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_forward_parameter(model, xt, 0.01)
+    DiffOpt.forward_differentiate!(model)
+    dθ1_dx = DiffOpt.get_forward_variable(model, θ1)
+    dθ2_dx = DiffOpt.get_forward_variable(model, θ2)
+
+    ## check first order approximation keeps solution close to target withing tolerance
+    θ_approx = [θ1̂ + dθ1_dx, θ1̂ + dθ2_dx]
+    x_approx = l1 * cos(θ_approx[1]) + l2 * cos(θ_approx[1] + θ_approx[2])
+    y_approx = l1 * sin(θ_approx[1]) + l2 * sin(θ_approx[1] + θ_approx[2])
+    _error = [x_approx - (x_t + 0.01), y_approx - y_t]
+    println("Error in first order approximation: ", _error)
+    feas[j, i] = norm(_error)
+
+    ## ---- forward diff wrt  yt  (∂θ/∂y) ----
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_forward_parameter(model, yt, 0.01)
+    DiffOpt.forward_differentiate!(model)
+    dθ1_dy = DiffOpt.get_forward_variable(model, θ1)
+    dθ2_dy = DiffOpt.get_forward_variable(model, θ2)
+
+    ## 2-norm of inverse Jacobian
+    Jinv = [
+        dθ1_dx dθ1_dy
+        dθ2_dx dθ2_dy
+    ]
+    heat[j, i] = opnorm(Jinv)            # σ_max  of Jinv
+end
+# Replace nans with 0.0
+heat = replace(heat, NaN => 0.0)
+
+# ## Results with Plot graphs
+
+default(;
+    size = (1150, 350),
+    legendfontsize = 8,
+    guidefontsize = 9,
+    tickfontsize = 7,
+)
+
+plt = heatmap(
+    xs,
+    ys,
+    heat;
+    xlabel = "x target",
+    ylabel = "y target",
+    clims = (0, quantile(skipmissing(heat), 0.95)),   # clip extremes
+    colorbar_title = "‖∂θ/∂(x,y)‖₂",
+    left_margin = 5Plots.Measures.mm,
+    bottom_margin = 5Plots.Measures.mm,
+);
+
+# Overlay workspace boundary
+θ = range(0, 2π; length = 200)
+plot!(plt, reach * cos.(θ), reach * sin.(θ); c = :white, lw = 1, lab = "reach");
+
+plt_feas = heatmap(
+    xs,
+    ys,
+    feas;
+    xlabel = "x target",
+    ylabel = "y target",
+    clims = (0, 1),
+    colorbar_title = "Precision Error",
+    left_margin = 5Plots.Measures.mm,
+    bottom_margin = 5Plots.Measures.mm,
+);
+
+plot!(
+    plt_feas,
+    reach * cos.(θ),
+    reach * sin.(θ);
+    c = :white,
+    lw = 1,
+    lab = "reach",
+);
+
+plt_all = plot(
+    plt,
+    plt_feas;
+    layout = (1, 2),
+    left_margin = 5Plots.Measures.mm,
+    bottom_margin = 5Plots.Measures.mm,
+    legend = :bottomright,
+)
+
+# Left figure shows the spectral-norm heat-map
+# $\bigl\lVert\partial\boldsymbol{\theta}/\partial(x,y)\bigr\rVert_2$
+# for a two-link arm - Bright rings mark near-singular poses. Right figure shows the normalized precision error of the first order approximation derived from calculated sensitivities.