Skip to content

Commit 43d6eb8

Browse files
andrewrosembergandrewrosemberg
authored
Merge branch 'master' into ar/dualparameter
2 parents f8042b9 + 23c27b6 commit 43d6eb8

23 files changed

+624
-119
lines changed

‎.github/workflows/format_check.yml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ jobs:
1818
shell: julia --color=yes {0}
1919
run: |
2020
using Pkg
21-
Pkg.add(PackageSpec(name="JuliaFormatter", version="1"))
21+
Pkg.add(PackageSpec(name="JuliaFormatter", version="2"))
2222
using JuliaFormatter
2323
format(".", verbose=true)
2424
out = String(read(Cmd(`git diff`)))

‎docs/src/examples/Mean_Variance_Portfolio_Example.jl

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
# # Mean Variance Portfolio Example
2+
3+
# 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}$:
4+
# ```math
5+
# \max_{x} \quad \mu^\top x
6+
# \quad\text{s.t.}\quad
7+
# x^\top \Sigma x \;\le\; v_{\max}, \quad
8+
# \mathbf{1}^\top x = 1,\quad
9+
# x \succeq 0,
10+
# ```
11+
# 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:
12+
13+
# ```math
14+
# \| \Sigma^{1/2} x \|_{2} \;\le\; \sigma_{\max}
15+
# ```
16+
# where $\Sigma^{1/2}$ is the Cholesky factorization of the covariance matrix and $\sigma_{\max}$ is the standard deviation limit.
17+
18+
# 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:
19+
# ```math
20+
# \frac{dL}{d\,\sigma_{\max}} \;=\;
21+
# \underbrace{\frac{\partial L}{\partial x}}_{\text{(1) decision impact}}\;
22+
# \cdot\;
23+
# \underbrace{\frac{\partial x^*}{\partial \sigma_{\max}}}_{\text{(2) from DiffOpt.jl}},
24+
# ```
25+
# where $x^*(\sigma_{\max})$ is the portfolio that solves the conic Markowitz problem under a given risk limit.
26+
27+
# ## Define and solve the Mean-Variance Portfolio Problem for a range of risk limits
28+
29+
# First, import the libraries.
30+
31+
using Test
32+
using JuMP
33+
import DiffOpt
34+
using LinearAlgebra
35+
import SCS
36+
using Plots
37+
using Plots.Measures
38+
39+
# Fixed data
40+
41+
# Training data (in-sample)
42+
Σ = [
43+
0.002 0.0005 0.001
44+
0.0005 0.003 0.0002
45+
0.001 0.0002 0.0025
46+
]
47+
μ_train = [0.05, 0.08, 0.12]
48+
49+
# Test data (out-of-sample)
50+
μ_test = [0.02, -0.3, 0.1] # simple forecast error example
51+
52+
# Sweep over σ_max
53+
σ_grid = 0.002:0.002:0.06
54+
N = length(σ_grid)
55+
56+
predicted_ret = zeros(N) # μ_train' * x*
57+
realised_ret = zeros(N) # μ_test' * x*
58+
loss = zeros(N) # L(x*)
59+
dL_dσ = zeros(N) # ∂L/∂σ_max from DiffOpt
60+
61+
for (k, σ_val) in enumerate(σ_grid)
62+
63+
## 1) differentiable conic model
64+
model = DiffOpt.conic_diff_model(SCS.Optimizer)
65+
set_silent(model)
66+
67+
## 2) parameter σ_max
68+
@variable(model, σ_max in Parameter(σ_val))
69+
70+
## 3) portfolio weights
71+
@variable(model, x[1:3] >= 0)
72+
@constraint(model, sum(x) <= 1)
73+
74+
## 4) objective: maximise expected return (training data)
75+
@objective(model, Max, dot(μ_train, x))
76+
77+
## 5) conic variance constraint ||L*x|| <= σ_max
78+
L_chol = cholesky(Symmetric(Σ)).L
79+
@variable(model, t >= 0)
80+
@constraint(model, [t; L_chol * x] in SecondOrderCone())
81+
@constraint(model, t <= σ_max)
82+
83+
optimize!(model)
84+
85+
x_opt = value.(x)
86+
println("Optimal portfolio weights: ", x_opt)
87+
88+
## store performance numbers
89+
predicted_ret[k] = dot(μ_train, x_opt)
90+
realised_ret[k] = dot(μ_test, x_opt)
91+
92+
## -------- reverse differentiation wrt σ_max --------
93+
DiffOpt.empty_input_sensitivities!(model)
94+
## ∂L/∂x (adjoint) = -μ_test
95+
DiffOpt.set_reverse_variable.(model, x, μ_test)
96+
DiffOpt.reverse_differentiate!(model)
97+
dL_dσ[k] = DiffOpt.get_reverse_parameter(model, σ_max)
98+
end
99+
100+
# ## Results with Plot graphs
101+
102+
default(;
103+
size = (1150, 350),
104+
legendfontsize = 8,
105+
guidefontsize = 9,
106+
tickfontsize = 7,
107+
)
108+
109+
# (a) predicted vs realised return
110+
plt_ret = plot(
111+
σ_grid,
112+
realised_ret;
113+
lw = 2,
114+
label = "Realised (test)",
115+
xlabel = "σ_max (risk limit)",
116+
ylabel = "Return",
117+
title = "Return vs risk limit",
118+
legend = :bottomright,
119+
);
120+
plot!(
121+
plt_ret,
122+
σ_grid,
123+
predicted_ret;
124+
lw = 2,
125+
ls = :dash,
126+
label = "Predicted (train)",
127+
);
128+
129+
# (b) out-of-sample loss and its gradient
130+
plt_loss = plot(
131+
σ_grid,
132+
dL_dσ;
133+
xlabel = "σ_max (risk limit)",
134+
ylabel = "∂L/∂σ_max",
135+
title = "Return Gradient",
136+
legend = false,
137+
);
138+
139+
plot_all = plot(
140+
plt_ret,
141+
plt_loss;
142+
layout = (1, 2),
143+
left_margin = 5Plots.Measures.mm,
144+
bottom_margin = 5Plots.Measures.mm,
145+
)
146+
147+
# Impact of the risk limit $\sigma_{\max}$ on Markowitz
148+
# portfolios. **Left:** predicted in-sample return versus
149+
# realized out-of-sample return. **Right:** the
150+
# out-of-sample loss $L(x)$ together with the absolute gradient
151+
# $|\partial L/\partial\sigma_{\max}|$ obtained from
152+
# `DiffOpt.jl`. The gradient tells the practitioner which
153+
# way—and how aggressively—to adjust $\sigma_{\max}$ to reduce
154+
# forecast error; its value is computed in one reverse-mode call
155+
# without re-solving the optimization for perturbed risk limits.

‎docs/src/examples/Planar_Arm_Example.jl

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
1+
# # Planar Arm Example
2+
3+
# 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:
4+
# ```math
5+
# f(\theta_1,\theta_2) = \bigl(\ell_1\cos(\theta_1) + \ell_2\cos(\theta_1+\theta_2),\,\,
6+
# \ell_1\sin(\theta_1) + \ell_2\sin(\theta_1+\theta_2)\bigr).
7+
# ```
8+
# We can solve an NLP:
9+
# ```math
10+
# \min_{\theta_1,\theta_2} \;\; (\theta_1^2 + \theta_2^2),
11+
# \quad\text{s.t.}\quad f(\theta_1,\theta_2) = (x_t,y_t).
12+
# ```
13+
# 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*.
14+
15+
# ## Define and solve the 2-link planar arm problem and build sensitivity map of joint angles to target
16+
17+
# First, import the libraries.
18+
19+
using Test
20+
using JuMP
21+
import DiffOpt
22+
using LinearAlgebra
23+
using Statistics
24+
import Ipopt
25+
using Plots
26+
using Plots.Measures
27+
28+
# Fixed data
29+
30+
# Arm geometry
31+
l1 = 1.0;
32+
l2 = 1.0;
33+
reach = l1 + l2 # 2.0
34+
tol = 1e-6 # numerical tolerance for feasibility
35+
# Sampling grid in workspace
36+
grid_res = 25 # grid resolution low for documentation compilation requirements
37+
xs = range(-reach, reach; length = grid_res)
38+
ys = range(-reach, reach; length = grid_res)
39+
40+
heat = fill(NaN, grid_res, grid_res) # store ‖J_inv‖₂
41+
feas = fill(0.0, grid_res, grid_res) # feasibility mask
42+
θ1mat = similar(heat)
43+
44+
function ik_angles(x, y; l1 = 1.0, l2 = 1.0, elbow_up = true)
45+
c2 = (x^2 + y^2 - l1^2 - l2^2) / (2 * l1 * l2)
46+
θ2 = elbow_up ? acos(clamp(c2, -1, 1)) : -acos(clamp(c2, -1, 1))
47+
k1 = l1 + l2 * cos(θ2)
48+
k2 = l2 * sin(θ2)
49+
θ1 = atan(y, x) - atan(k2, k1)
50+
return θ1, θ2
51+
end
52+
53+
for (i, x_t) in enumerate(xs), (j, y_t) in enumerate(ys)
54+
global θ1mat, heat
55+
## skip points outside the circular reach
56+
norm([x_t, y_t]) > reach && continue
57+
58+
## ---------- build differentiable NLP ----------
59+
model = DiffOpt.nonlinear_diff_model(Ipopt.Optimizer)
60+
set_silent(model)
61+
62+
@variable(model, xt in Parameter(x_t))
63+
@variable(model, yt in Parameter(y_t))
64+
@variable(model, θ1)
65+
@variable(model, θ2)
66+
67+
@objective(model, Min, θ1^2 + θ2^2)
68+
@constraint(model, l1 * cos(θ1) + l2 * cos(θ1 + θ2) == xt)
69+
@constraint(model, l1 * sin(θ1) + l2 * sin(θ1 + θ2) == yt)
70+
71+
## --- supply analytic start values ---
72+
θ1₀, θ2₀ = ik_angles(x_t, y_t; elbow_up = true)
73+
set_start_value(θ1, θ1₀)
74+
set_start_value(θ2, θ2₀)
75+
76+
optimize!(model)
77+
println("Solving for target (", x_t, ", ", y_t, ")")
78+
## check for optimality
79+
status = termination_status(model)
80+
println("Status: ", status)
81+
82+
status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || continue
83+
84+
θ1̂ = value(θ1)
85+
θ1mat[j, i] = θ1̂ # save pose
86+
87+
## ---- forward diff wrt xt (∂θ/∂x) ----
88+
DiffOpt.empty_input_sensitivities!(model)
89+
DiffOpt.set_forward_parameter(model, xt, 0.01)
90+
DiffOpt.forward_differentiate!(model)
91+
dθ1_dx = DiffOpt.get_forward_variable(model, θ1)
92+
dθ2_dx = DiffOpt.get_forward_variable(model, θ2)
93+
94+
## check first order approximation keeps solution close to target withing tolerance
95+
θ_approx = [θ1̂ + dθ1_dx, θ1̂ + dθ2_dx]
96+
x_approx = l1 * cos(θ_approx[1]) + l2 * cos(θ_approx[1] + θ_approx[2])
97+
y_approx = l1 * sin(θ_approx[1]) + l2 * sin(θ_approx[1] + θ_approx[2])
98+
_error = [x_approx - (x_t + 0.01), y_approx - y_t]
99+
println("Error in first order approximation: ", _error)
100+
feas[j, i] = norm(_error)
101+
102+
## ---- forward diff wrt yt (∂θ/∂y) ----
103+
DiffOpt.empty_input_sensitivities!(model)
104+
DiffOpt.set_forward_parameter(model, yt, 0.01)
105+
DiffOpt.forward_differentiate!(model)
106+
dθ1_dy = DiffOpt.get_forward_variable(model, θ1)
107+
dθ2_dy = DiffOpt.get_forward_variable(model, θ2)
108+
109+
## 2-norm of inverse Jacobian
110+
Jinv = [
111+
dθ1_dx dθ1_dy
112+
dθ2_dx dθ2_dy
113+
]
114+
heat[j, i] = opnorm(Jinv) # σ_max of Jinv
115+
end
116+
# Replace nans with 0.0
117+
heat = replace(heat, NaN => 0.0)
118+
119+
# ## Results with Plot graphs
120+
121+
default(;
122+
size = (1150, 350),
123+
legendfontsize = 8,
124+
guidefontsize = 9,
125+
tickfontsize = 7,
126+
)
127+
128+
plt = heatmap(
129+
xs,
130+
ys,
131+
heat;
132+
xlabel = "x target",
133+
ylabel = "y target",
134+
clims = (0, quantile(skipmissing(heat), 0.95)), # clip extremes
135+
colorbar_title = "‖∂θ/∂(x,y)‖₂",
136+
left_margin = 5Plots.Measures.mm,
137+
bottom_margin = 5Plots.Measures.mm,
138+
);
139+
140+
# Overlay workspace boundary
141+
θ = range(0, 2π; length = 200)
142+
plot!(plt, reach * cos.(θ), reach * sin.(θ); c = :white, lw = 1, lab = "reach");
143+
144+
plt_feas = heatmap(
145+
xs,
146+
ys,
147+
feas;
148+
xlabel = "x target",
149+
ylabel = "y target",
150+
clims = (0, 1),
151+
colorbar_title = "Precision Error",
152+
left_margin = 5Plots.Measures.mm,
153+
bottom_margin = 5Plots.Measures.mm,
154+
);
155+
156+
plot!(
157+
plt_feas,
158+
reach * cos.(θ),
159+
reach * sin.(θ);
160+
c = :white,
161+
lw = 1,
162+
lab = "reach",
163+
);
164+
165+
plt_all = plot(
166+
plt,
167+
plt_feas;
168+
layout = (1, 2),
169+
left_margin = 5Plots.Measures.mm,
170+
bottom_margin = 5Plots.Measures.mm,
171+
legend = :bottomright,
172+
)
173+
174+
# Left figure shows the spectral-norm heat-map
175+
# $\bigl\lVert\partial\boldsymbol{\theta}/\partial(x,y)\bigr\rVert_2$
176+
# 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.

‎docs/src/examples/Thermal_Generation_Dispatch_Example.jl

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -102,7 +102,7 @@ function diff_reverse(model::Model, ϵ::Float64 = 1.0)
102102
perturbation = zeros(I + 1)
103103

104104
## Loop for each primal variable
105-
for i in 1:I+1
105+
for i in 1:(I+1)
106106
## Set the perturbation in the Primal Variables and set the context to Backward
107107
perturbation[i] = ϵ
108108
MOI.set.(model, DiffOpt.ReverseVariablePrimal(), vect_ref, perturbation)
@@ -156,7 +156,7 @@ end
156156
# Result Primal Values:
157157
Plots.plot(
158158
d,
159-
data_results[1, :, 1:I+1];
159+
data_results[1, :, 1:(I+1)];
160160
title = "Generation by Demand",
161161
label = ["Thermal Generation 1" "Thermal Generation 2" "Thermal Generation 3" "Generation Deficit"],
162162
xlabel = "Demand [unit]",
@@ -166,7 +166,7 @@ Plots.plot(
166166
# Result Sensitivity Analysis:
167167
Plots.plot(
168168
d,
169-
data_results[1, :, I+2:2*(I+1)];
169+
data_results[1, :, (I+2):(2*(I+1))];
170170
title = "Sensitivity of Generation by Demand",
171171
label = ["T. Gen. 1 Sensitivity" "T. Gen. 2 Sensitivity" "T. Gen. 3 Sensitivity" "Gen. Deficit Sensitivity"],
172172
xlabel = "Demand [unit]",
@@ -177,7 +177,7 @@ Plots.plot(
177177
# Result Primal Values:
178178
Plots.plot(
179179
d,
180-
data_results[2, :, 1:I+1];
180+
data_results[2, :, 1:(I+1)];
181181
title = "Generation by Demand",
182182
label = ["Thermal Generation 1" "Thermal Generation 2" "Thermal Generation 3" "Generation Deficit"],
183183
xlabel = "Demand [unit]",
@@ -187,7 +187,7 @@ Plots.plot(
187187
# Result Sensitivity Analysis:
188188
Plots.plot(
189189
d,
190-
data_results[2, :, I+2:2*(I+1)];
190+
data_results[2, :, (I+2):(2*(I+1))];
191191
title = "Sensitivity of Generation by Demand",
192192
label = ["T. Gen. 1 Sensitivity" "T. Gen. 2 Sensitivity" "T. Gen. 3 Sensitivity" "Gen. Deficit Sensitivity"],
193193
xlabel = "Demand [unit]",

0 commit comments

Comments
 (0)