From 76e5a30b4a6b412614cc499b0f169add2683009b Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 16:57:41 +0200 Subject: [PATCH 01/17] Extract state constraint example to separate file - Create new autonomous example-state-constraint.md with state constraint variant - Remove state constraint section from example-double-integrator-energy.md - Add cross-references between the two examples - Add state constraint example to Basic Examples menu in make.jl --- docs/make.jl | 3 +- docs/src/example-double-integrator-energy.md | 147 +---------- docs/src/example-state-constraint.md | 244 +++++++++++++++++++ 3 files changed, 247 insertions(+), 147 deletions(-) create mode 100644 docs/src/example-state-constraint.md diff --git a/docs/make.jl b/docs/make.jl index b9dcd6a1c..3dd601fe3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -151,7 +151,7 @@ cp( Draft = false ``` =# -draft = false # Draft mode: if true, @example blocks in markdown are not executed +draft = true # Draft mode: if true, @example blocks in markdown are not executed # ═══════════════════════════════════════════════════════════════════════════════ # Load extensions @@ -207,6 +207,7 @@ with_api_reference(src_dir, ext_dir) do api_pages "Time mininimisation" => "example-double-integrator-time.md", "Control-free problems" => "example-control-free.md", "Singular control" => "example-singular-control.md", + "State constraint" => "example-state-constraint.md", ], "Manual" => [ "Define a problem" => "manual-abstract.md", diff --git a/docs/src/example-double-integrator-energy.md b/docs/src/example-double-integrator-energy.md index 54c88cbd5..59b67be91 100644 --- a/docs/src/example-double-integrator-energy.md +++ b/docs/src/example-double-integrator-energy.md @@ -166,149 +166,4 @@ plot(indirect_sol) - You can use [MINPACK.jl](@extref Tutorials Resolution-of-the-shooting-equation) instead of [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve). - For more details about the flow construction, visit the [Compute flows from optimal control problems](@ref manual-flow-ocp) page. - In this simple example, we have set an arbitrary initial guess. It can be helpful to use the solution of the direct method to initialise the shooting method. See the [Goddard tutorial](@extref Tutorials tutorial-goddard) for such a concrete application. - -## State constraint - -The following example illustrates both direct and indirect solution approaches for the energy minimization problem with a state constraint on the maximal velocity. The workflow demonstrates a practical strategy: a direct method on a coarse grid first identifies the problem structure and provides an initial guess for the indirect method, which then computes a precise solution via shooting based on Pontryagin's Maximum Principle. - -!!! note - - The direct solution can be refined using a finer discretization grid for higher accuracy. - -### Direct method: constrained case - -We add the path constraint - -```math - v(t) \le 1.2. -``` - -Let us model, solve and plot the optimal control problem with this constraint. - -```@example main -# the upper bound for v -v_max = 1.2 - -# the optimal control problem -ocp = @def begin - t ∈ [t0, tf], time - x = (q, v) ∈ R², state - u ∈ R, control - - v(t) ≤ v_max # state constraint - - x(t0) == x0 - x(tf) == xf - - ẋ(t) == [v(t), u(t)] - - 0.5∫( u(t)^2 ) → min -end - -# solve with a direct method -direct_sol = solve(ocp; grid_size=50) - -# plot the solution -plt = plot(direct_sol; label="Direct", size=(800, 600)) -``` - -The solution has three phases (unconstrained-constrained-unconstrained arcs), requiring definition of Hamiltonian flows for each phase and a shooting function to enforce boundary and switching conditions. - -### Indirect method: constrained case - -Under the normal case, the pseudo-Hamiltonian reads: - -```math -H(x, p, u, \mu) = p_1 v + p_2 u - \frac{u^2}{2} + \mu\, g(x), -``` - -where $g(x) = v_{\max} - v$. Along a boundary arc we have $g(x(t)) = 0$; differentiating gives: - -```math - \frac{\mathrm{d}}{\mathrm{d}t}g(x(t)) = -\dot{v}(t) = -u(t) = 0. -``` - -The zero control maximises the Hamiltonian, so $p_2(t) = 0$ along that arc. From the adjoint equation we then have - -```math - \dot{p}_2(t) = -p_1(t) + \mu(t) = 0 \quad \Rightarrow \mu(t) = p_1(t). -``` - -Because the adjoint vector is continuous at both the entry time $t_1$ and the exit time $t_2$, the unknowns are $p_0 \in \mathbb{R}^2$ together with $t_1$ and $t_2$. The target condition supplies two equations, $g(x(t_1)) = 0$ enforces the state constraint, and $p_2(t_1) = 0$ encodes the switching condition. - -```@example main -# flow for unconstrained extremals -f_interior = Flow(ocp, (x, p) -> p[2]) - -ub = 0 # boundary control -g(x) = v_max - x[2] # constraint: g(x) ≥ 0 -μ(p) = p[1] # dual variable - -# flow for boundary extremals -f_boundary = Flow(ocp, (x, p) -> ub, (x, u) -> g(x), (x, p) -> μ(p)) - -# shooting function -function shoot!(s, p0, t1, t2) - x_t0, p_t0 = x0, p0 - x_t1, p_t1 = f_interior(t0, x_t0, p_t0, t1) - x_t2, p_t2 = f_boundary(t1, x_t1, p_t1, t2) - x_tf, p_tf = f_interior(t2, x_t2, p_t2, tf) - s[1:2] = x_tf - xf - s[3] = g(x_t1) - s[4] = p_t1[2] -end -nothing # hide -``` - -We can derive an initial guess for the costate and the entry/exit times from the direct solution: - -```@example main -t = time_grid(direct_sol) # the time grid as a vector -x = state(direct_sol) # the state as a function of time -p = costate(direct_sol) # the costate as a function of time - -# initial costate -p0 = p(t0) - -# times where constraint is active -t12 = t[ 0 .≤ (g ∘ x).(t) .≤ 1e-3 ] - -# entry and exit times -t1 = minimum(t12) # entry time -t2 = maximum(t12) # exit time -nothing # hide -``` - -We can now solve the shooting equations. - -```@example main -# auxiliary in-place NLE function -nle!(s, ξ, _) = shoot!(s, ξ[1:2], ξ[3], ξ[4]) - -# initial guess for the Newton solver -ξ_guess = [p0..., t1, t2] - -# NLE problem with initial guess -prob = NonlinearProblem(nle!, ξ_guess) - -# resolution of the shooting equations -shooting_sol = solve(prob; show_trace=Val(true)) -p0, t1, t2 = shooting_sol.u[1:2], shooting_sol.u[3], shooting_sol.u[4] - -# print the costate solution and the entry and exit times -println("\np0 = ", p0, "\nt1 = ", t1, "\nt2 = ", t2) -``` - -To reconstruct the constrained trajectory, concatenate the flows as follows: an unconstrained arc until $t_1$, a boundary arc from $t_1$ to $t_2$, and a final unconstrained arc from $t_2$ to $t_f$. -This composition yields the full solution (state, costate, and control), which we then plot alongside the direct method for comparison. - -```@example main -# concatenation of the flows -φ = f_interior * (t1, f_boundary) * (t2, f_interior) - -# compute the solution: state, costate, control... -indirect_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 100)) - -# plot the solution on the previous plot -plot!(plt, indirect_sol; label="Indirect", color=2, linestyle=:dash) -``` + - For a version with a state constraint on the velocity, see the [State constraint](@ref example-state-constraint) example. diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md new file mode 100644 index 000000000..f6ed5d605 --- /dev/null +++ b/docs/src/example-state-constraint.md @@ -0,0 +1,244 @@ +# [State constraint](@id example-state-constraint) + +```@meta +Draft = false +``` + +This example illustrates the energy minimization problem for the double integrator with a state constraint on the maximal velocity. It demonstrates both direct and indirect solution approaches, showing how to handle state constraints using practical strategies. + +Let us consider a wagon moving along a rail, whose acceleration can be controlled by a force $u$. +We denote by $x = (q, v)$ the state of the wagon, where $q$ is the position and $v$ the velocity. + +```@raw html + +``` + +We assume that the mass is constant and equal to one, and that there is no friction. The dynamics are given by + +```math + \dot q(t) = v(t), \quad \dot v(t) = u(t),\quad u(t) \in \R, +``` + +which is simply the [double integrator](https://en.wikipedia.org/w/index.php?title=Double_integrator&oldid=1071399674) system. Let us consider a transfer starting at time $t_0 = 0$ and ending at time $t_f = 1$, for which we want to minimise the transfer energy + +```math + \frac{1}{2}\int_{0}^{1} u^2(t) \, \mathrm{d}t +``` + +starting from $x(0) = (-1, 0)$ and aiming to reach the target $x(1) = (0, 0)$. + +First, we need to import the [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) package to define the optimal control problem, [NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl) to solve it, and [Plots.jl](https://docs.juliaplots.org) to visualise the solution. + +```@example main +using OptimalControl +using NLPModelsIpopt +using Plots +``` + +## Optimal control problem + +Let us define the problem with the [`@def`](@ref) macro: + +```@raw html +
+
+``` + +```@example main +t0 = 0; tf = 1; x0 = [-1, 0]; xf = [0, 0] + +ocp = @def begin + t ∈ [t0, tf], time + x = (q, v) ∈ R², state + u ∈ R, control + + x(t0) == x0 + x(tf) == xf + + ẋ(t) == [v(t), u(t)] + + 0.5∫( u(t)^2 ) → min +end +nothing # hide +``` + +```@raw html +
+
+``` + +### Mathematical formulation + +```math + \begin{aligned} + & \text{Minimise} && \frac{1}{2}\int_0^1 u^2(t) \,\mathrm{d}t \\ + & \text{subject to} \\ + & && \dot{x}(t) = [v(t), u(t)], \\[1.0em] + & && x(0) = (-1,0), \\[0.5em] + & && x(1) = (0,0). + \end{aligned} +``` + +```@raw html +
+
+``` + +!!! note "Nota bene" + + For a comprehensive introduction to the syntax used above to define the optimal control problem, see [this abstract syntax tutorial](@ref manual-abstract-syntax). In particular, non-Unicode alternatives are available for derivatives, integrals, *etc.* + +## Adding a state constraint + +We now add a path constraint on the maximal velocity: + +```math + v(t) \le 1.2. +``` + +The workflow demonstrates a practical strategy: a direct method on a coarse grid first identifies the problem structure and provides an initial guess for the indirect method, which then computes a precise solution via shooting based on Pontryagin's Maximum Principle. + +!!! note + + The direct solution can be refined using a finer discretization grid for higher accuracy. + +### Direct method: constrained case + +Let us model, solve and plot the optimal control problem with this constraint. + +```@example main +# the upper bound for v +v_max = 1.2 + +# the optimal control problem +ocp = @def begin + t ∈ [t0, tf], time + x = (q, v) ∈ R², state + u ∈ R, control + + v(t) ≤ v_max # state constraint + + x(t0) == x0 + x(tf) == xf + + ẋ(t) == [v(t), u(t)] + + 0.5∫( u(t)^2 ) → min +end + +# solve with a direct method +direct_sol = solve(ocp; grid_size=50) + +# plot the solution +plt = plot(direct_sol; label="Direct", size=(800, 600)) +``` + +The solution has three phases (unconstrained-constrained-unconstrained arcs), requiring definition of Hamiltonian flows for each phase and a shooting function to enforce boundary and switching conditions. + +### Indirect method: constrained case + +Under the normal case, the pseudo-Hamiltonian reads: + +```math +H(x, p, u, \mu) = p_1 v + p_2 u - \frac{u^2}{2} + \mu\, g(x), +``` + +where $g(x) = v_{\max} - v$. Along a boundary arc we have $g(x(t)) = 0$; differentiating gives: + +```math + \frac{\mathrm{d}}{\mathrm{d}t}g(x(t)) = -\dot{v}(t) = -u(t) = 0. +``` + +The zero control maximises the Hamiltonian, so $p_2(t) = 0$ along that arc. From the adjoint equation we then have + +```math + \dot{p}_2(t) = -p_1(t) + \mu(t) = 0 \quad \Rightarrow \mu(t) = p_1(t). +``` + +Because the adjoint vector is continuous at both the entry time $t_1$ and the exit time $t_2$, the unknowns are $p_0 \in \mathbb{R}^2$ together with $t_1$ and $t_2$. The target condition supplies two equations, $g(x(t_1)) = 0$ enforces the state constraint, and $p_2(t_1) = 0$ encodes the switching condition. + +```@example main +using OrdinaryDiffEq # Ordinary Differential Equations (ODE) solver +using NonlinearSolve # Nonlinear Equations (NLE) solver + +# flow for unconstrained extremals +f_interior = Flow(ocp, (x, p) -> p[2]) + +ub = 0 # boundary control +g(x) = v_max - x[2] # constraint: g(x) ≥ 0 +μ(p) = p[1] # dual variable + +# flow for boundary extremals +f_boundary = Flow(ocp, (x, p) -> ub, (x, u) -> g(x), (x, p) -> μ(p)) + +# shooting function +function shoot!(s, p0, t1, t2) + x_t0, p_t0 = x0, p0 + x_t1, p_t1 = f_interior(t0, x_t0, p_t0, t1) + x_t2, p_t2 = f_boundary(t1, x_t1, p_t1, t2) + x_tf, p_tf = f_interior(t2, x_t2, p_t2, tf) + s[1:2] = x_tf - xf + s[3] = g(x_t1) + s[4] = p_t1[2] +end +nothing # hide +``` + +We can derive an initial guess for the costate and the entry/exit times from the direct solution: + +```@example main +t = time_grid(direct_sol) # the time grid as a vector +x = state(direct_sol) # the state as a function of time +p = costate(direct_sol) # the costate as a function of time + +# initial costate +p0 = p(t0) + +# times where constraint is active +t12 = t[ 0 .≤ (g ∘ x).(t) .≤ 1e-3 ] + +# entry and exit times +t1 = minimum(t12) # entry time +t2 = maximum(t12) # exit time +nothing # hide +``` + +We can now solve the shooting equations. + +```@example main +# auxiliary in-place NLE function +nle!(s, ξ, _) = shoot!(s, ξ[1:2], ξ[3], ξ[4]) + +# initial guess for the Newton solver +ξ_guess = [p0..., t1, t2] + +# NLE problem with initial guess +prob = NonlinearProblem(nle!, ξ_guess) + +# resolution of the shooting equations +shooting_sol = solve(prob; show_trace=Val(true)) +p0, t1, t2 = shooting_sol.u[1:2], shooting_sol.u[3], shooting_sol.u[4] + +# print the costate solution and the entry and exit times +println("\np0 = ", p0, "\nt1 = ", t1, "\nt2 = ", t2) +``` + +To reconstruct the constrained trajectory, concatenate the flows as follows: an unconstrained arc until $t_1$, a boundary arc from $t_1$ to $t_2$, and a final unconstrained arc from $t_2$ to $t_f$. +This composition yields the full solution (state, costate, and control), which we then plot alongside the direct method for comparison. + +```@example main +# concatenation of the flows +φ = f_interior * (t1, f_boundary) * (t2, f_interior) + +# compute the solution: state, costate, control... +indirect_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 100)) + +# plot the solution on the previous plot +plot!(plt, indirect_sol; label="Indirect", color=2, linestyle=:dash) +``` + +!!! note + + - You can use [MINPACK.jl](@extref Tutorials Resolution-of-the-shooting-equation) instead of [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve). + - For more details about the flow construction, visit the [Compute flows from optimal control problems](@ref manual-flow-ocp) page. + - For the unconstrained version of this problem, see the [Energy minimisation](@ref example-double-integrator-energy) example. From d1305b996049794f3bf26405daaa760872816f44 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 18:17:48 +0200 Subject: [PATCH 02/17] Add second-order state constraint (Bryson-Denham) to state constraint example - Rename 'Adding a state constraint' -> 'First-order state constraint' - Add explanation of constraint order for order-1 case - Add new 'Second-order state constraint' section (Bryson-Denham problem) - Explain order-2 via two differentiations of g(x) - Show two cases: touch point (a=1/4) and boundary arc (a=1/9) - Direct method only, solutions superimposed on same plot - Add Bryson et al. 1963 reference --- docs/src/example-state-constraint.md | 102 ++++++++++++++++++++++++++- 1 file changed, 100 insertions(+), 2 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index f6ed5d605..735a9d3c8 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -4,7 +4,7 @@ Draft = false ``` -This example illustrates the energy minimization problem for the double integrator with a state constraint on the maximal velocity. It demonstrates both direct and indirect solution approaches, showing how to handle state constraints using practical strategies. +This example illustrates how state constraints of different orders affect the structure of optimal solutions for the double integrator energy minimization problem. It demonstrates both direct and indirect solution approaches. Let us consider a wagon moving along a rail, whose acceleration can be controlled by a force $u$. We denote by $x = (q, v)$ the state of the wagon, where $q$ is the position and $v$ the velocity. @@ -88,7 +88,7 @@ nothing # hide For a comprehensive introduction to the syntax used above to define the optimal control problem, see [this abstract syntax tutorial](@ref manual-abstract-syntax). In particular, non-Unicode alternatives are available for derivatives, integrals, *etc.* -## Adding a state constraint +## First-order state constraint We now add a path constraint on the maximal velocity: @@ -96,6 +96,14 @@ We now add a path constraint on the maximal velocity: v(t) \le 1.2. ``` +This is a **first-order state constraint**: differentiating $g(x) = v_{\max} - v$ once already makes the control appear, + +```math + \frac{\mathrm{d}}{\mathrm{d}t}g(x(t)) = -\dot{v}(t) = -u(t), +``` + +which fixes $u = 0$ on the boundary arc. + The workflow demonstrates a practical strategy: a direct method on a coarse grid first identifies the problem structure and provides an initial guess for the indirect method, which then computes a precise solution via shooting based on Pontryagin's Maximum Principle. !!! note @@ -242,3 +250,93 @@ plot!(plt, indirect_sol; label="Indirect", color=2, linestyle=:dash) - You can use [MINPACK.jl](@extref Tutorials Resolution-of-the-shooting-equation) instead of [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve). - For more details about the flow construction, visit the [Compute flows from optimal control problems](@ref manual-flow-ocp) page. - For the unconstrained version of this problem, see the [Energy minimisation](@ref example-double-integrator-energy) example. + +## Second-order state constraint + +We now consider the same double integrator with different boundary conditions and a constraint on the **position** $x_1 = q$:[^2] + +```math + q(t) \le a. +``` + +The boundary conditions are $x(0) = (0, 1)$ and $x(1) = (0, -1)$. + +This is a **second-order state constraint**: the control $u$ appears only after differentiating $g(x) = a - q$ twice, + +```math + \frac{\mathrm{d}}{\mathrm{d}t}g(x(t)) = -\dot{q}(t) = -v(t) \quad \text{(no control)}, +``` + +```math + \frac{\mathrm{d}^2}{\mathrm{d}t^2}g(x(t)) = -\dot{v}(t) = -u(t) \quad \text{(control appears)}. +``` + +On a boundary arc where $g(x(t)) = 0$, both derivatives must vanish, forcing $v(t) = 0$ and $u(t) = 0$. + +### Solution structure + +The unconstrained optimal trajectory for these boundary conditions is $q(t) = t - t^2$, which reaches its maximum $1/4$ at $t = 1/2$. The solution structure depends on $a$: + +- $a \ge 1/4$: constraint never active (unconstrained solution) +- $a = 1/4$: **touch point** — the trajectory touches $q = a$ at a single instant $t = 1/2$ +- $a < 1/4$: **boundary arc** — the trajectory remains on $q = a$ for a finite time interval + +### Direct method + +We compare the two constrained cases using the direct method. + +```@example main +# new boundary conditions +x0_bd = [0.0, 1.0]; xf_bd = [0.0, -1.0] + +# Case 1: touch point (a = 1/4) +a_touch = 1/4 + +ocp_touch = @def begin + t ∈ [t0, tf], time + x = (q, v) ∈ R², state + u ∈ R, control + + q(t) ≤ a_touch + + x(t0) == x0_bd + x(tf) == xf_bd + + ẋ(t) == [v(t), u(t)] + + 0.5∫( u(t)^2 ) → min +end + +sol_touch = solve(ocp_touch; grid_size=100) +nothing # hide +``` + +```@example main +# Case 2: boundary arc (a = 1/9) +a_arc = 1/9 + +ocp_arc = @def begin + t ∈ [t0, tf], time + x = (q, v) ∈ R², state + u ∈ R, control + + q(t) ≤ a_arc + + x(t0) == x0_bd + x(tf) == xf_bd + + ẋ(t) == [v(t), u(t)] + + 0.5∫( u(t)^2 ) → min +end + +sol_arc = solve(ocp_arc; grid_size=100) +nothing # hide +``` + +```@example main +plt_bd = plot(sol_touch; label="Touch point (a = 1/4)", size=(800, 600)) +plot!(plt_bd, sol_arc; label="Boundary arc (a = 1/9)", color=2, linestyle=:dash) +``` + +[^2]: Bryson, A.E., Denham, W.F., & Dreyfus, S.E. (1963). *Optimal programming problems with inequality constraints I: necessary conditions for extremal solutions*. AIAA Journal. [doi.org/10.2514/3.2107](https://doi.org/10.2514/3.2107) From 15e91a35e491a70c5a354a7fa869f1a259ececc2 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 18:21:49 +0200 Subject: [PATCH 03/17] Add references to state constraint theory in intro (Bryson 1963, Jacobson 1971) --- docs/src/example-state-constraint.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 735a9d3c8..4dade0a57 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -4,7 +4,7 @@ Draft = false ``` -This example illustrates how state constraints of different orders affect the structure of optimal solutions for the double integrator energy minimization problem. It demonstrates both direct and indirect solution approaches. +This example illustrates how state constraints of different orders affect the structure of optimal solutions for the double integrator energy minimization problem. It demonstrates both direct and indirect solution approaches. Some examples with state constraints of different orders are solved analytically in Bryson et al.[^1] and Jacobson et al.[^2]. Let us consider a wagon moving along a rail, whose acceleration can be controlled by a force $u$. We denote by $x = (q, v)$ the state of the wagon, where $q$ is the position and $v$ the velocity. @@ -253,7 +253,7 @@ plot!(plt, indirect_sol; label="Indirect", color=2, linestyle=:dash) ## Second-order state constraint -We now consider the same double integrator with different boundary conditions and a constraint on the **position** $x_1 = q$:[^2] +We now consider the same double integrator with different boundary conditions and a constraint on the **position** $x_1 = q$:[^1] ```math q(t) \le a. @@ -339,4 +339,6 @@ plt_bd = plot(sol_touch; label="Touch point (a = 1/4)", size=(800, 600)) plot!(plt_bd, sol_arc; label="Boundary arc (a = 1/9)", color=2, linestyle=:dash) ``` -[^2]: Bryson, A.E., Denham, W.F., & Dreyfus, S.E. (1963). *Optimal programming problems with inequality constraints I: necessary conditions for extremal solutions*. AIAA Journal. [doi.org/10.2514/3.2107](https://doi.org/10.2514/3.2107) +[^1]: Bryson, A.E., Denham, W.F., & Dreyfus, S.E. (1963). *Optimal programming problems with inequality constraints I: necessary conditions for extremal solutions*. AIAA Journal, 1(11), 2544–2550. [doi.org/10.2514/3.2107](https://doi.org/10.2514/3.2107) + +[^2]: Jacobson, D.H., Lele, M.M., & Speyer, J.L. (1971). *New necessary conditions of optimality for control problems with state-variable inequality constraints*. Journal of Mathematical Analysis and Applications, 35, 255–284. From 438887122638ffc834c0ff6eddd36d611d0170d2 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 18:25:32 +0200 Subject: [PATCH 04/17] Fix list items to not start with math expressions --- docs/src/example-state-constraint.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 4dade0a57..1a524b553 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -277,9 +277,9 @@ On a boundary arc where $g(x(t)) = 0$, both derivatives must vanish, forcing $v( The unconstrained optimal trajectory for these boundary conditions is $q(t) = t - t^2$, which reaches its maximum $1/4$ at $t = 1/2$. The solution structure depends on $a$: -- $a \ge 1/4$: constraint never active (unconstrained solution) -- $a = 1/4$: **touch point** — the trajectory touches $q = a$ at a single instant $t = 1/2$ -- $a < 1/4$: **boundary arc** — the trajectory remains on $q = a$ for a finite time interval +- For large values ($a \ge 1/4$): the constraint is never active (unconstrained solution) +- **Touch point** ($a = 1/4$): the trajectory touches $q = a$ at a single instant $t = 1/2$ +- **Boundary arc** ($a < 1/4$): the trajectory remains on $q = a$ for a finite time interval ### Direct method From d30620c2ec5b64f26ecf4b68ac83d3d1bc7fde90 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 18:41:14 +0200 Subject: [PATCH 05/17] Fix order-2 solution structure and add Bryson & Ho 1975 reference - Correct three-regime description: unconstrained (a>=1/4), touch point (1/6<=a<=1/4), boundary arc (a<1/6) - Use a=1/5 (touch point) and a=1/7 (boundary arc) as example values - Add Bryson & Ho (1975) Applied Optimal Control as [^3] --- docs/src/example-state-constraint.md | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 1a524b553..0d057833b 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -275,22 +275,22 @@ On a boundary arc where $g(x(t)) = 0$, both derivatives must vanish, forcing $v( ### Solution structure -The unconstrained optimal trajectory for these boundary conditions is $q(t) = t - t^2$, which reaches its maximum $1/4$ at $t = 1/2$. The solution structure depends on $a$: +The unconstrained optimal trajectory for these boundary conditions is $q(t) = t - t^2$, which reaches its maximum $1/4$ at $t = 1/2$. A characteristic feature of second-order state constraints is the existence of an intermediate regime between the unconstrained and boundary-arc cases[^3]. The solution structure depends on $a$: -- For large values ($a \ge 1/4$): the constraint is never active (unconstrained solution) -- **Touch point** ($a = 1/4$): the trajectory touches $q = a$ at a single instant $t = 1/2$ -- **Boundary arc** ($a < 1/4$): the trajectory remains on $q = a$ for a finite time interval +- **Unconstrained** ($a \ge 1/4$): the constraint is never active +- **Touch point** ($1/6 \le a \le 1/4$): the trajectory touches $q = a$ at a single instant, without sliding along the boundary +- **Boundary arc** ($a < 1/6$): the trajectory remains on $q = a$ for a finite time interval, during which $v(t) = 0$ and $u(t) = 0$ ### Direct method -We compare the two constrained cases using the direct method. +We compare the two constrained cases using the direct method, taking $a = 1/5$ (touch point) and $a = 1/7$ (boundary arc). ```@example main # new boundary conditions x0_bd = [0.0, 1.0]; xf_bd = [0.0, -1.0] -# Case 1: touch point (a = 1/4) -a_touch = 1/4 +# Case 1: touch point (a = 1/5, in range [1/6, 1/4]) +a_touch = 1/5 ocp_touch = @def begin t ∈ [t0, tf], time @@ -307,13 +307,13 @@ ocp_touch = @def begin 0.5∫( u(t)^2 ) → min end -sol_touch = solve(ocp_touch; grid_size=100) +sol_touch = solve(ocp_touch; grid_size=100, display=false) nothing # hide ``` ```@example main -# Case 2: boundary arc (a = 1/9) -a_arc = 1/9 +# Case 2: boundary arc (a = 1/7, below 1/6) +a_arc = 1/7 ocp_arc = @def begin t ∈ [t0, tf], time @@ -330,15 +330,17 @@ ocp_arc = @def begin 0.5∫( u(t)^2 ) → min end -sol_arc = solve(ocp_arc; grid_size=100) +sol_arc = solve(ocp_arc; grid_size=100, display=false) nothing # hide ``` ```@example main -plt_bd = plot(sol_touch; label="Touch point (a = 1/4)", size=(800, 600)) -plot!(plt_bd, sol_arc; label="Boundary arc (a = 1/9)", color=2, linestyle=:dash) +plt_bd = plot(sol_touch; label="Touch point (a = 1/5)", size=(800, 600)) +plot!(plt_bd, sol_arc; label="Boundary arc (a = 1/7)", color=2, linestyle=:dash) ``` [^1]: Bryson, A.E., Denham, W.F., & Dreyfus, S.E. (1963). *Optimal programming problems with inequality constraints I: necessary conditions for extremal solutions*. AIAA Journal, 1(11), 2544–2550. [doi.org/10.2514/3.2107](https://doi.org/10.2514/3.2107) [^2]: Jacobson, D.H., Lele, M.M., & Speyer, J.L. (1971). *New necessary conditions of optimality for control problems with state-variable inequality constraints*. Journal of Mathematical Analysis and Applications, 35, 255–284. + +[^3]: Bryson, A.E. & Ho, Y.-C. (1975). *Applied Optimal Control: Optimization, Estimation and Control*. CRC Press. From 5b65abf84ddafdf0360dc36f259f75565a3417fb Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 18:44:40 +0200 Subject: [PATCH 06/17] Refactor order-2 section: use make_ocp(a) function to avoid duplication --- docs/src/example-state-constraint.md | 50 ++++++++-------------------- 1 file changed, 14 insertions(+), 36 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 0d057833b..12d0c8650 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -289,52 +289,30 @@ We compare the two constrained cases using the direct method, taking $a = 1/5$ ( # new boundary conditions x0_bd = [0.0, 1.0]; xf_bd = [0.0, -1.0] -# Case 1: touch point (a = 1/5, in range [1/6, 1/4]) -a_touch = 1/5 +# parametric OCP: double integrator with position constraint q(t) ≤ a +function make_ocp(a) + @def begin + t ∈ [t0, tf], time + x = (q, v) ∈ R², state + u ∈ R, control -ocp_touch = @def begin - t ∈ [t0, tf], time - x = (q, v) ∈ R², state - u ∈ R, control + q(t) ≤ a - q(t) ≤ a_touch + x(t0) == x0_bd + x(tf) == xf_bd - x(t0) == x0_bd - x(tf) == xf_bd + ẋ(t) == [v(t), u(t)] - ẋ(t) == [v(t), u(t)] - - 0.5∫( u(t)^2 ) → min + 0.5∫( u(t)^2 ) → min + end end - -sol_touch = solve(ocp_touch; grid_size=100, display=false) nothing # hide ``` ```@example main -# Case 2: boundary arc (a = 1/7, below 1/6) -a_arc = 1/7 - -ocp_arc = @def begin - t ∈ [t0, tf], time - x = (q, v) ∈ R², state - u ∈ R, control +sol_touch = solve(make_ocp(1/5); grid_size=100, display=false) # touch point +sol_arc = solve(make_ocp(1/7); grid_size=100, display=false) # boundary arc - q(t) ≤ a_arc - - x(t0) == x0_bd - x(tf) == xf_bd - - ẋ(t) == [v(t), u(t)] - - 0.5∫( u(t)^2 ) → min -end - -sol_arc = solve(ocp_arc; grid_size=100, display=false) -nothing # hide -``` - -```@example main plt_bd = plot(sol_touch; label="Touch point (a = 1/5)", size=(800, 600)) plot!(plt_bd, sol_arc; label="Boundary arc (a = 1/7)", color=2, linestyle=:dash) ``` From 9ccd68bdd30f66cdfef889ae20f868518178b491 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 19:16:25 +0200 Subject: [PATCH 07/17] Add indirect method for touch point case (a=0.2) in order-2 section MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Two unconstrained arcs joined at t1 with costate jump [Δpq, 0] - Shooting unknowns: p0 (2D), t1, Δpq (4 equations) - Initial guess extracted from direct solution sol_touch - Flow concatenation: fs_bd * (t1, [Δpq, 0], fs_bd) - Analytical verification: pq = ±4.8, jump = 9.6, cost = 2.24 --- docs/src/example-state-constraint.md | 105 +++++++++++++++++++++++++-- 1 file changed, 100 insertions(+), 5 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 12d0c8650..0236d816f 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -283,7 +283,7 @@ The unconstrained optimal trajectory for these boundary conditions is $q(t) = t ### Direct method -We compare the two constrained cases using the direct method, taking $a = 1/5$ (touch point) and $a = 1/7$ (boundary arc). +We compare the two constrained cases using the direct method, taking $a = 0.2$ (touch point) and $a = 0.1$ (boundary arc). ```@example main # new boundary conditions @@ -310,11 +310,106 @@ nothing # hide ``` ```@example main -sol_touch = solve(make_ocp(1/5); grid_size=100, display=false) # touch point -sol_arc = solve(make_ocp(1/7); grid_size=100, display=false) # boundary arc +sol_touch = solve(make_ocp(0.2); grid_size=100, display=false) # touch point +sol_arc = solve(make_ocp(0.1); grid_size=100, display=false) # boundary arc + +state_style = (legend=false, ) +costate_style = (legend=false, ) +plt_bd = plot( + sol_touch; + label="a = 0.2", + size=(800, 600), + state_style=state_style, + costate_style=costate_style, +) +plot!( + plt_bd, + sol_arc; + label="a = 0.1", + color=2, + linestyle=:dash, + state_style=state_style, + costate_style=costate_style, +) +``` + +### Indirect method: touch point case + +For the touch point case ($a = 0.2$), the optimal solution consists of two unconstrained arcs on $[t_0, t_1]$ and $[t_1, t_f]$, joined at the contact instant $t_1$ where $q(t_1) = a$ and $v(t_1) = 0$. The costate is discontinuous at $t_1$: the first component $p_q$ undergoes a jump $\Delta p_q$ while $p_v$ remains continuous. + +The shooting unknowns are therefore the initial costate $p_0 \in \mathbb{R}^2$, the contact time $t_1$, and the costate jump $\Delta p_q$. The four shooting conditions are: + +```math +x(t_f) = x_f, \quad q(t_1) = a, \quad v(t_1) = 0. +``` + +```@example main +a_touch = 0.2 + +# interior (unconstrained) flow +fs_bd = Flow(make_ocp(a_touch), (x, p) -> p[2]) + +# constraint: g(x) = a - q ≥ 0 +g_bd(x) = a_touch - x[1] + +# shooting function: unknowns p0 (2D), t1 (contact time), Δpq (costate jump) +function shoot_touch!(s, p0, t1, Δpq) + x_t1, p_t1 = fs_bd(t0, x0_bd, p0, t1) # arc 1: t0 → t1 + p_t1_plus = [p_t1[1] + Δpq, p_t1[2]] # costate jump at t1 + x_tf, p_tf = fs_bd(t1, x_t1, p_t1_plus, tf) # arc 2: t1 → tf + s[1:2] = x_tf - xf_bd # reach target + s[3] = g_bd(x_t1) # touch: q(t1) = a + s[4] = x_t1[2] # tangency: v(t1) = 0 +end +nothing # hide +``` + +We extract the initial guess from the direct solution `sol_touch`. + +```@example main +t_grid = time_grid(sol_touch) +x_sol = state(sol_touch) +p_sol = costate(sol_touch) + +p0_guess = p_sol(t0) + +# t1: time where q(t) is closest to the constraint bound a +t1_guess = t_grid[argmin(abs.(g_bd.(x_sol.(t_grid))))] + +# Δpq: estimated costate jump around t1 +ε = 0.05 * (tf - t0) +Δpq_guess = p_sol(t1_guess + ε)[1] - p_sol(t1_guess - ε)[1] + +println("p0 guess = ", p0_guess) +println("t1 guess = ", t1_guess) +println("Δpq guess = ", Δpq_guess) +nothing # hide +``` + +```@example main +nle_touch!(s, ξ, _) = shoot_touch!(s, ξ[1:2], ξ[3], ξ[4]) + +ξ_guess = [p0_guess..., t1_guess, Δpq_guess] +sol_shoot_touch = solve(NonlinearProblem(nle_touch!, ξ_guess); show_trace=Val(true)) + +p0_touch = sol_shoot_touch.u[1:2] +t1_touch = sol_shoot_touch.u[3] +Δpq_touch = sol_shoot_touch.u[4] + +println("\np0 = ", p0_touch, "\nt1 = ", t1_touch, "\nΔpq = ", Δpq_touch) +``` + +The analytical solution gives $p_q = -4.8$ on $[t_0, t_1)$, $p_q = +4.8$ on $(t_1, t_f]$, with a jump of $9.6$ and an optimal cost of $2.24$. + +```@example main +# concatenate: arc 1 → costate jump → arc 2 +f_touch = fs_bd * (t1_touch, [Δpq_touch, 0.0], fs_bd) + +# reconstruct the indirect solution +indirect_touch = f_touch((t0, tf), x0_bd, p0_touch; saveat=range(t0, tf, 100)) -plt_bd = plot(sol_touch; label="Touch point (a = 1/5)", size=(800, 600)) -plot!(plt_bd, sol_arc; label="Boundary arc (a = 1/7)", color=2, linestyle=:dash) +plot(indirect_touch; label="Indirect", size=(800, 600), + state_style=(legend=false,), costate_style=(legend=false,)) ``` [^1]: Bryson, A.E., Denham, W.F., & Dreyfus, S.E. (1963). *Optimal programming problems with inequality constraints I: necessary conditions for extremal solutions*. AIAA Journal, 1(11), 2544–2550. [doi.org/10.2514/3.2107](https://doi.org/10.2514/3.2107) From b1e3cc70b6c25d0c2cb22bc510b67d7ed900ca33 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 19:21:12 +0200 Subject: [PATCH 08/17] Add t1 = 1/2 to analytical solution and minor code cleanup --- docs/src/example-state-constraint.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 0236d816f..23c15e021 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -354,12 +354,12 @@ g_bd(x) = a_touch - x[1] # shooting function: unknowns p0 (2D), t1 (contact time), Δpq (costate jump) function shoot_touch!(s, p0, t1, Δpq) - x_t1, p_t1 = fs_bd(t0, x0_bd, p0, t1) # arc 1: t0 → t1 + x_t1, p_t1 = fs_bd(t0, x0_bd, p0, t1) # arc 1: t0 → t1 p_t1_plus = [p_t1[1] + Δpq, p_t1[2]] # costate jump at t1 - x_tf, p_tf = fs_bd(t1, x_t1, p_t1_plus, tf) # arc 2: t1 → tf - s[1:2] = x_tf - xf_bd # reach target - s[3] = g_bd(x_t1) # touch: q(t1) = a - s[4] = x_t1[2] # tangency: v(t1) = 0 + x_tf, _ = fs_bd(t1, x_t1, p_t1_plus, tf) # arc 2: t1 → tf + s[1:2] = x_tf - xf_bd # reach target + s[3] = g_bd(x_t1) # touch: q(t1) = a + s[4] = x_t1[2] # tangency: v(t1) = 0 end nothing # hide ``` @@ -399,7 +399,7 @@ t1_touch = sol_shoot_touch.u[3] println("\np0 = ", p0_touch, "\nt1 = ", t1_touch, "\nΔpq = ", Δpq_touch) ``` -The analytical solution gives $p_q = -4.8$ on $[t_0, t_1)$, $p_q = +4.8$ on $(t_1, t_f]$, with a jump of $9.6$ and an optimal cost of $2.24$. +The analytical solution gives $t_1 = 1/2$, $p_q = -4.8$ on $[t_0, t_1)$, $p_q = +4.8$ on $(t_1, t_f]$, with a jump of $9.6$ and an optimal cost of $2.24$. ```@example main # concatenate: arc 1 → costate jump → arc 2 From 767771b0b2096db21ec6687298e50cfeaf02f794 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 20:37:57 +0200 Subject: [PATCH 09/17] Add indirect method for boundary arc case (a=0.1) in order-2 section MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Three arcs: unconstrained → boundary arc → unconstrained - Flow: fs_arc * (t1, [Δpq1,0], fc_bd) * (t2, [Δpq2,0], fs_arc) - 6 unknowns: p0 (2D), t1, t2, Δpq1, Δpq2 - 6 conditions: x(tf)=xf (×2), q(t1)=a, v(t1)=0, pv(t1+)=0, pq(t1+)=0 - Initial guess extracted from sol_arc - Symmetry check: Δpq1 = Δpq2, t2 = 1 - t1 --- docs/src/example-state-constraint.md | 98 +++++++++++++++++++++++++++- 1 file changed, 96 insertions(+), 2 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 23c15e021..513c95be9 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -317,7 +317,7 @@ state_style = (legend=false, ) costate_style = (legend=false, ) plt_bd = plot( sol_touch; - label="a = 0.2", + label="a = 0.2", size=(800, 600), state_style=state_style, costate_style=costate_style, @@ -408,7 +408,101 @@ f_touch = fs_bd * (t1_touch, [Δpq_touch, 0.0], fs_bd) # reconstruct the indirect solution indirect_touch = f_touch((t0, tf), x0_bd, p0_touch; saveat=range(t0, tf, 100)) -plot(indirect_touch; label="Indirect", size=(800, 600), +plot(indirect_touch; label="Indirect (a = 0.2)", size=(800, 600), + state_style=(legend=false,), costate_style=(legend=false,)) +``` + +### Indirect method: boundary arc case + +For the boundary arc case ($a = 0.1$), the optimal solution consists of three arcs: two unconstrained arcs on $[t_0, t_1]$ and $[t_2, t_f]$, separated by a boundary arc on $[t_1, t_2]$ where $q(t) = a$ and $v(t) = 0$. Along the boundary arc, the control is $u = 0$, the costate satisfies $p_q = 0$ and $p_v = 0$, and the Lagrange multiplier $\mu = 0$. The costate has jumps $[\Delta p_q^1, 0]$ and $[\Delta p_q^2, 0]$ at $t_1$ and $t_2$ respectively. + +The six shooting unknowns are the initial costate $p_0 \in \mathbb{R}^2$, the entry and exit times $t_1$ and $t_2$, and the two jumps $\Delta p_q^1$ and $\Delta p_q^2$. The shooting conditions are: + +```math +x(t_f) = x_f, \quad q(t_1) = a, \quad v(t_1) = 0, \quad p_v(t_1^+) = 0, \quad p_q(t_1^+) = 0. +``` + +```@example main +a_arc = 0.1 + +# interior (unconstrained) flow +fs_arc = Flow(make_ocp(a_arc), (x, p) -> p[2]) + +# boundary arc flow: u = 0, constraint g(x) = a - q ≥ 0, multiplier μ = 0 +fc_bd = Flow(make_ocp(a_arc), (x, p) -> 0, (x, u) -> a_arc - x[1], (x, p) -> 0) + +# constraint function +g_arc(x) = a_arc - x[1] + +# shooting function: unknowns p0 (2D), t1, t2, Δpq1, Δpq2 +function shoot_arc!(s, p0, t1, t2, Δpq1, Δpq2) + x_t1, p_t1 = fs_arc(t0, x0_bd, p0, t1) # arc 1: t0 → t1 + p_t1_plus = [p_t1[1] + Δpq1, p_t1[2]] # costate jump at t1 + x_t2, p_t2 = fc_bd(t1, x_t1, p_t1_plus, t2) # arc 2: t1 → t2 (boundary) + p_t2_plus = [p_t2[1] + Δpq2, p_t2[2]] # costate jump at t2 + x_tf, _ = fs_arc(t2, x_t2, p_t2_plus, tf) # arc 3: t2 → tf + s[1:2] = x_tf - xf_bd # reach target + s[3] = g_arc(x_t1) # touch: q(t1) = a + s[4] = x_t1[2] # tangency: v(t1) = 0 + s[5] = p_t1_plus[2] # switching: pv(t1+) = 0 + s[6] = p_t1_plus[1] # arc condition: pq(t1+) = 0 +end +nothing # hide +``` + +We extract the initial guess from the direct solution `sol_arc`. + +```@example main +t_grid_arc = time_grid(sol_arc) +x_sol_arc = state(sol_arc) +p_sol_arc = costate(sol_arc) + +p0_guess_arc = p_sol_arc(t0) + +# t1, t2: entry and exit of the boundary arc (q ≈ a) +active = findall(t -> g_arc(x_sol_arc(t)) ≤ 1e-3, t_grid_arc) +t1_guess_arc = t_grid_arc[first(active)] +t2_guess_arc = t_grid_arc[last(active)] + +# jumps: costate difference around t1 and t2 +ε_arc = 0.05 * (tf - t0) +Δpq1_guess = p_sol_arc(t1_guess_arc + ε_arc)[1] - p_sol_arc(t1_guess_arc - ε_arc)[1] +Δpq2_guess = p_sol_arc(t2_guess_arc + ε_arc)[1] - p_sol_arc(t2_guess_arc - ε_arc)[1] + +println("p0 guess = ", p0_guess_arc) +println("t1 guess = ", t1_guess_arc) +println("t2 guess = ", t2_guess_arc) +println("Δpq1 guess = ", Δpq1_guess) +println("Δpq2 guess = ", Δpq2_guess) +nothing # hide +``` + +```@example main +nle_arc!(s, ξ, _) = shoot_arc!(s, ξ[1:2], ξ[3], ξ[4], ξ[5], ξ[6]) + +ξ_guess_arc = [p0_guess_arc..., t1_guess_arc, t2_guess_arc, Δpq1_guess, Δpq2_guess] +sol_shoot_arc = solve(NonlinearProblem(nle_arc!, ξ_guess_arc); show_trace=Val(true)) + +p0_arc = sol_shoot_arc.u[1:2] +t1_arc = sol_shoot_arc.u[3] +t2_arc = sol_shoot_arc.u[4] +Δpq1 = sol_shoot_arc.u[5] +Δpq2 = sol_shoot_arc.u[6] + +println("\np0 = ", p0_arc) +println("t1 = ", t1_arc, " (expect ", 3a_arc, ")") +println("t2 = ", t2_arc, " (expect ", 1 - 3a_arc, ")") +println("Δpq1 = ", Δpq1, " Δpq2 = ", Δpq2, " (expect equal by symmetry)") +``` + +```@example main +# concatenate: arc 1 → jump → boundary arc → jump → arc 3 +f_arc = fs_arc * (t1_arc, [Δpq1, 0.0], fc_bd) * (t2_arc, [Δpq2, 0.0], fs_arc) + +# reconstruct the indirect solution +indirect_arc = f_arc((t0, tf), x0_bd, p0_arc; saveat=range(t0, tf, 100)) + +plot(indirect_arc; label="Indirect (a = 0.1)", size=(800, 600), state_style=(legend=false,), costate_style=(legend=false,)) ``` From 7deee6b5e3aeb792925262604011c78eaae399e6 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 20:39:58 +0200 Subject: [PATCH 10/17] Expand boundary arc description with adjoint chain explanation --- docs/src/example-state-constraint.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 513c95be9..e80dec955 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -414,7 +414,7 @@ plot(indirect_touch; label="Indirect (a = 0.2)", size=(800, 600), ### Indirect method: boundary arc case -For the boundary arc case ($a = 0.1$), the optimal solution consists of three arcs: two unconstrained arcs on $[t_0, t_1]$ and $[t_2, t_f]$, separated by a boundary arc on $[t_1, t_2]$ where $q(t) = a$ and $v(t) = 0$. Along the boundary arc, the control is $u = 0$, the costate satisfies $p_q = 0$ and $p_v = 0$, and the Lagrange multiplier $\mu = 0$. The costate has jumps $[\Delta p_q^1, 0]$ and $[\Delta p_q^2, 0]$ at $t_1$ and $t_2$ respectively. +For the boundary arc case ($a = 0.1$), the optimal solution consists of three arcs: two unconstrained arcs on $[t_0, t_1]$ and $[t_2, t_f]$, separated by a boundary arc on $[t_1, t_2]$ where $q(t) = a$ and $v(t) = 0$. Along the boundary arc, the control is $u = 0$, since differentiating $g(x) = a - q \geq 0$ twice gives $\ddot{q} = u = 0$. From the maximisation condition, $p_v(t) = 0$ along the arc. Differentiating the adjoint equation $\dot{p}_v = -p_q$ and using $p_v = 0$ yields $p_q = 0$. Differentiating further gives $\mu = \dot{p}_q = 0$. The costate has jumps $[\Delta p_q^1, 0]$ and $[\Delta p_q^2, 0]$ at $t_1$ and $t_2$ respectively. The six shooting unknowns are the initial costate $p_0 \in \mathbb{R}^2$, the entry and exit times $t_1$ and $t_2$, and the two jumps $\Delta p_q^1$ and $\Delta p_q^2$. The shooting conditions are: From 8b26b488691dca3c9c654cb8e597fcd38dec561c Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 20:45:57 +0200 Subject: [PATCH 11/17] Add Hamiltonian expression to boundary arc description with adapted notation --- docs/src/example-state-constraint.md | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index e80dec955..597e02e64 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -414,7 +414,13 @@ plot(indirect_touch; label="Indirect (a = 0.2)", size=(800, 600), ### Indirect method: boundary arc case -For the boundary arc case ($a = 0.1$), the optimal solution consists of three arcs: two unconstrained arcs on $[t_0, t_1]$ and $[t_2, t_f]$, separated by a boundary arc on $[t_1, t_2]$ where $q(t) = a$ and $v(t) = 0$. Along the boundary arc, the control is $u = 0$, since differentiating $g(x) = a - q \geq 0$ twice gives $\ddot{q} = u = 0$. From the maximisation condition, $p_v(t) = 0$ along the arc. Differentiating the adjoint equation $\dot{p}_v = -p_q$ and using $p_v = 0$ yields $p_q = 0$. Differentiating further gives $\mu = \dot{p}_q = 0$. The costate has jumps $[\Delta p_q^1, 0]$ and $[\Delta p_q^2, 0]$ at $t_1$ and $t_2$ respectively. +For the boundary arc case ($a = 0.1$), the optimal solution consists of three arcs: two unconstrained arcs on $[t_0, t_1]$ and $[t_2, t_f]$, separated by a boundary arc on $[t_1, t_2]$ where $q(t) = a$ and $v(t) = 0$. The pseudo-Hamiltonian is + +```math +H(x, p, u, \mu) = p_q\, v + p_v\, u + 0.5\, p^0 u^2 + \mu\, g(x), +``` + +where $p^0 = -1$ in the normal case and $g(x) = a - q \geq 0$ is the constraint. Along the boundary arc, the control is $u = 0$, since differentiating $g(x) = a - q \geq 0$ twice gives $\ddot{q} = u = 0$. From the maximisation condition, $p_v(t) = 0$ along the arc. Differentiating the adjoint equation $\dot{p}_v = -p_q$ and using $p_v = 0$ yields $p_q = 0$. Differentiating further gives $\mu = \dot{p}_q = 0$. The costate has jumps $[\Delta p_q^1, 0]$ and $[\Delta p_q^2, 0]$ at $t_1$ and $t_2$ respectively. The six shooting unknowns are the initial costate $p_0 \in \mathbb{R}^2$, the entry and exit times $t_1$ and $t_2$, and the two jumps $\Delta p_q^1$ and $\Delta p_q^2$. The shooting conditions are: From ec05944c53518afa3a83d9a195daaaa9f36db3f7 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 20:54:44 +0200 Subject: [PATCH 12/17] Superpose indirect arc solution on touch point plot using plot --- docs/src/example-state-constraint.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 597e02e64..63ebe786d 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -408,8 +408,8 @@ f_touch = fs_bd * (t1_touch, [Δpq_touch, 0.0], fs_bd) # reconstruct the indirect solution indirect_touch = f_touch((t0, tf), x0_bd, p0_touch; saveat=range(t0, tf, 100)) -plot(indirect_touch; label="Indirect (a = 0.2)", size=(800, 600), - state_style=(legend=false,), costate_style=(legend=false,)) +plt_indirect = plot(indirect_touch; label="Indirect (a = 0.2)", size=(800, 600), + state_style=(legend=false,), costate_style=(legend=false,)) ``` ### Indirect method: boundary arc case @@ -420,7 +420,7 @@ For the boundary arc case ($a = 0.1$), the optimal solution consists of three ar H(x, p, u, \mu) = p_q\, v + p_v\, u + 0.5\, p^0 u^2 + \mu\, g(x), ``` -where $p^0 = -1$ in the normal case and $g(x) = a - q \geq 0$ is the constraint. Along the boundary arc, the control is $u = 0$, since differentiating $g(x) = a - q \geq 0$ twice gives $\ddot{q} = u = 0$. From the maximisation condition, $p_v(t) = 0$ along the arc. Differentiating the adjoint equation $\dot{p}_v = -p_q$ and using $p_v = 0$ yields $p_q = 0$. Differentiating further gives $\mu = \dot{p}_q = 0$. The costate has jumps $[\Delta p_q^1, 0]$ and $[\Delta p_q^2, 0]$ at $t_1$ and $t_2$ respectively. +where $p^0 = -1$ in the normal case and $g(x) = a - q \geq 0$ is the constraint. Along the boundary arc, the control is $u = 0$, since differentiating $g(x) = a - q \geq 0$ twice gives $\ddot{q} = u = 0$. From the maximisation condition, $p_v(t) = 0$ along the arc. Differentiating the adjoint equation $\dot{p}_v = -p_q$ and using $p_v = 0$ yields $p_q = 0$. Differentiating further gives $\mu = \dot{p}_q = 0$. The costate has jumps $(\Delta p_q^1, 0)$ and $(\Delta p_q^2, 0)$ at $t_1$ and $t_2$ respectively. The six shooting unknowns are the initial costate $p_0 \in \mathbb{R}^2$, the entry and exit times $t_1$ and $t_2$, and the two jumps $\Delta p_q^1$ and $\Delta p_q^2$. The shooting conditions are: @@ -442,16 +442,16 @@ g_arc(x) = a_arc - x[1] # shooting function: unknowns p0 (2D), t1, t2, Δpq1, Δpq2 function shoot_arc!(s, p0, t1, t2, Δpq1, Δpq2) - x_t1, p_t1 = fs_arc(t0, x0_bd, p0, t1) # arc 1: t0 → t1 + x_t1, p_t1 = fs_arc(t0, x0_bd, p0, t1) # arc 1: t0 → t1 p_t1_plus = [p_t1[1] + Δpq1, p_t1[2]] # costate jump at t1 - x_t2, p_t2 = fc_bd(t1, x_t1, p_t1_plus, t2) # arc 2: t1 → t2 (boundary) + x_t2, p_t2 = fc_bd(t1, x_t1, p_t1_plus, t2) # arc 2: t1 → t2 (boundary) p_t2_plus = [p_t2[1] + Δpq2, p_t2[2]] # costate jump at t2 x_tf, _ = fs_arc(t2, x_t2, p_t2_plus, tf) # arc 3: t2 → tf - s[1:2] = x_tf - xf_bd # reach target - s[3] = g_arc(x_t1) # touch: q(t1) = a - s[4] = x_t1[2] # tangency: v(t1) = 0 - s[5] = p_t1_plus[2] # switching: pv(t1+) = 0 - s[6] = p_t1_plus[1] # arc condition: pq(t1+) = 0 + s[1:2] = x_tf - xf_bd # reach target + s[3] = g_arc(x_t1) # touch: q(t1) = a + s[4] = x_t1[2] # tangency: v(t1) = 0 + s[5] = p_t1_plus[2] # switching: pv(t1+) = 0 + s[6] = p_t1_plus[1] # arc condition: pq(t1+) = 0 end nothing # hide ``` @@ -508,8 +508,8 @@ f_arc = fs_arc * (t1_arc, [Δpq1, 0.0], fc_bd) * (t2_arc, [Δpq2, 0.0], fs_arc) # reconstruct the indirect solution indirect_arc = f_arc((t0, tf), x0_bd, p0_arc; saveat=range(t0, tf, 100)) -plot(indirect_arc; label="Indirect (a = 0.1)", size=(800, 600), - state_style=(legend=false,), costate_style=(legend=false,)) +plot!(plt_indirect, indirect_arc; label="Indirect (a = 0.1)", color=2, linestyle=:dash, + state_style=(legend=false,), costate_style=(legend=false,)) ``` [^1]: Bryson, A.E., Denham, W.F., & Dreyfus, S.E. (1963). *Optimal programming problems with inequality constraints I: necessary conditions for extremal solutions*. AIAA Journal, 1(11), 2544–2550. [doi.org/10.2514/3.2107](https://doi.org/10.2514/3.2107) From a4d58ae792eb2245f7f329bc87726e4fed99a378 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 20:55:17 +0200 Subject: [PATCH 13/17] Add transition text between make_ocp definition and solve cells --- docs/src/example-state-constraint.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 63ebe786d..a11ab279f 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -309,6 +309,8 @@ end nothing # hide ``` +We now solve both cases using this parametric OCP definition. + ```@example main sol_touch = solve(make_ocp(0.2); grid_size=100, display=false) # touch point sol_arc = solve(make_ocp(0.1); grid_size=100, display=false) # boundary arc From d69cdc0992ed4386dd376c16a9eadc830d3e933f Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 20:57:24 +0200 Subject: [PATCH 14/17] Use findall pattern for entry/exit time extraction in order-1 case --- docs/src/example-state-constraint.md | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index a11ab279f..6a094b836 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -188,6 +188,7 @@ function shoot!(s, p0, t1, t2) s[1:2] = x_tf - xf s[3] = g(x_t1) s[4] = p_t1[2] + return nothing end nothing # hide ``` @@ -202,12 +203,10 @@ p = costate(direct_sol) # the costate as a function of time # initial costate p0 = p(t0) -# times where constraint is active -t12 = t[ 0 .≤ (g ∘ x).(t) .≤ 1e-3 ] - -# entry and exit times -t1 = minimum(t12) # entry time -t2 = maximum(t12) # exit time +# t1, t2: entry and exit of the constrained arc (v ≈ v_max) +active = findall(t -> 0 ≤ g(x(t)) ≤ 1e-3, t) +t1 = t[first(active)] # entry time +t2 = t[last(active)] # exit time nothing # hide ``` @@ -362,6 +361,7 @@ function shoot_touch!(s, p0, t1, Δpq) s[1:2] = x_tf - xf_bd # reach target s[3] = g_bd(x_t1) # touch: q(t1) = a s[4] = x_t1[2] # tangency: v(t1) = 0 + return nothing end nothing # hide ``` @@ -454,6 +454,7 @@ function shoot_arc!(s, p0, t1, t2, Δpq1, Δpq2) s[4] = x_t1[2] # tangency: v(t1) = 0 s[5] = p_t1_plus[2] # switching: pv(t1+) = 0 s[6] = p_t1_plus[1] # arc condition: pq(t1+) = 0 + return nothing end nothing # hide ``` @@ -468,7 +469,7 @@ p_sol_arc = costate(sol_arc) p0_guess_arc = p_sol_arc(t0) # t1, t2: entry and exit of the boundary arc (q ≈ a) -active = findall(t -> g_arc(x_sol_arc(t)) ≤ 1e-3, t_grid_arc) +active = findall(t -> 0 ≤ g_arc(x_sol_arc(t)) ≤ 1e-3, t_grid_arc) t1_guess_arc = t_grid_arc[first(active)] t2_guess_arc = t_grid_arc[last(active)] From f0404bf23280a66413d0846b2427d755a19493f7 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 21:00:38 +0200 Subject: [PATCH 15/17] Final changes: adjust epsilon for arc guess and update docs Project.toml --- docs/src/assets/Project.toml | 2 +- docs/src/example-state-constraint.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml index 519de2ab4..22c049a82 100644 --- a/docs/src/assets/Project.toml +++ b/docs/src/assets/Project.toml @@ -35,7 +35,7 @@ CTFlows = "0.8" CTModels = "0.9" CTParser = "0.8" CTSolvers = "0.4" -CUDA = "5, 6" +CUDA = "5" CommonSolve = "0.2" DataFrames = "1" Documenter = "1" diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index 6a094b836..cc88d53b4 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -474,7 +474,7 @@ t1_guess_arc = t_grid_arc[first(active)] t2_guess_arc = t_grid_arc[last(active)] # jumps: costate difference around t1 and t2 -ε_arc = 0.05 * (tf - t0) +ε_arc = 0.1 * (tf - t0) Δpq1_guess = p_sol_arc(t1_guess_arc + ε_arc)[1] - p_sol_arc(t1_guess_arc - ε_arc)[1] Δpq2_guess = p_sol_arc(t2_guess_arc + ε_arc)[1] - p_sol_arc(t2_guess_arc - ε_arc)[1] From 8f97c032f48d6ada9ba7d78abd50cf43f63f6493 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 21:06:08 +0200 Subject: [PATCH 16/17] Update CHANGELOG for v2.0.2 with state constraint example documentation --- CHANGELOG.md | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9df7ab59c..8f721c031 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,28 @@ Versions follow [Semantic Versioning](https://semver.org/spec/v2.0.0.html). --- +## [2.0.2] — 2026-04-14 + +### Added + +- **Documentation examples**: + - New comprehensive state constraint example (`example-state-constraint.md`) demonstrating first-order and second-order (Bryson-Denham) state constraints + - Direct method implementation with parametric OCP for both touch point and boundary arc cases + - Indirect methods for touch point (2-arc) and boundary arc (3-arc) cases with shooting functions + - Theoretical references: Bryson et al. (1963), Jacobson et al. (1971), Bryson & Ho (1975) + - Hamiltonian-based adjoint chain explanations for boundary arc dynamics + +### Changed + +- **Documentation organization**: + - Extracted state constraint section from `example-double-integrator-energy.md` into dedicated example file + - Added cross-references between energy minimization and state constraint examples + +- **Dependencies**: + - Updated UnoSolver from v0.2 to v0.3 + +--- + ## [2.0.1] — 2026-04-13 ### Changed @@ -22,15 +44,6 @@ Versions follow [Semantic Versioning](https://semver.org/spec/v2.0.0.html). --- -## [2.0.2] — 2026-04-14 - -### Changed - -- **Dependencies**: - - Updated UnoSolver from v0.2 to v0.3 - ---- - ## [2.0.0] — 2026-04-03 **Major version release** with complete solve architecture redesign. This release introduces breaking changes from v1.1.6 (last stable release). See [BREAKING.md](BREAKING.md) for detailed migration guide. From 796ab288688af5aaab9870db938c48967a6f6e29 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 14 Apr 2026 21:06:42 +0200 Subject: [PATCH 17/17] Enable draft mode for state constraint example --- docs/make.jl | 2 +- docs/src/example-state-constraint.md | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 3dd601fe3..093efa034 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -151,7 +151,7 @@ cp( Draft = false ``` =# -draft = true # Draft mode: if true, @example blocks in markdown are not executed +draft = false # Draft mode: if true, @example blocks in markdown are not executed # ═══════════════════════════════════════════════════════════════════════════════ # Load extensions diff --git a/docs/src/example-state-constraint.md b/docs/src/example-state-constraint.md index cc88d53b4..5de30d409 100644 --- a/docs/src/example-state-constraint.md +++ b/docs/src/example-state-constraint.md @@ -1,9 +1,5 @@ # [State constraint](@id example-state-constraint) -```@meta -Draft = false -``` - This example illustrates how state constraints of different orders affect the structure of optimal solutions for the double integrator energy minimization problem. It demonstrates both direct and indirect solution approaches. Some examples with state constraints of different orders are solved analytically in Bryson et al.[^1] and Jacobson et al.[^2]. Let us consider a wagon moving along a rail, whose acceleration can be controlled by a force $u$.