Merge pull request #758 from control-toolbox/state-constraint

Add comprehensive state constraint example with Bryson-Denham problem

Author: ocots

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.

docs/make.jl

@@ -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",

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"

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.