add brachistochrone

ext/Descriptions/balanced_field.md

@@ -0,0 +1,48 @@
+The **Aircraft Balanced Field Length Calculation** is a classic aeronautical engineering problem used to determine the minimum runway length required for a safe takeoff or a safe stop in the event of an engine failure.  
+This implementation focuses on the **takeoff climb** phase with one engine out, starting from the critical engine failure speed $V_1$.
+
+### Mathematical formulation
+
+The problem is to minimise the final range $r(t_f)$ to reach the screen height (usually 35 ft).  
+The state vector is $x(t) = [r(t), v(t), h(t), \gamma(t)]^	op$ and the control is the angle of attack $\alpha(t)$.
+
+```math
+\begin{aligned}
+\min_{\alpha, t_f} \quad & r(t_f) \[0.5em]
+	ext{s.t.} \quad & \dot{r}(t) = v(t) \cos \gamma(t), \[0.5em]
+& \dot{v}(t) = \frac{T \cos \alpha(t) - D}{m} - g \sin \gamma(t), \[0.5em]
+& \dot{h}(t) = v(t) \sin \gamma(t), \[0.5em]
+& \dot{\gamma}(t) = \frac{T \sin \alpha(t) + L}{m v(t)} - \frac{g \cos \gamma(t)}{v(t)}, \[0.5em]
+& h(t_f) = 10.668 	ext{ m (35 ft)}, \[0.5em]
+& \gamma(t_f) = 5^\circ.
+\end{aligned}
+```
+
+The aerodynamic forces $L$ and $D$ depend on the angle of attack $\alpha$, velocity $v$, and altitude $h$ (ground effect).
+
+### Parameters
+
+The parameters are based on a mid-size jet aircraft (nominal values from Dymos):
+
+| Parameter | Symbol | Value |
+|-----------|--------|-------|
+| Mass | $m$ | 79015.8 kg |
+| Surface | $S$ | 124.7 m² |
+| Thrust (1 engine) | $T$ | 120102.0 N |
+| Screen height | $h_{tf}$ | 10.668 m |
+
+### Characteristics
+
+- Multi-state nonlinear dynamics,
+- Free final time $t_f$,
+- Ground effect inclusion in the drag model ($K$ coefficient),
+- Control bounds on the angle of attack $\alpha$,
+- Boundary conditions representing $V_1$ conditions and final climb requirements.
+
+### References
+
+- **Dymos Documentation**. [*Aircraft Balanced Field Length Calculation*.](https://openmdao.github.io/dymos/examples/balanced_field/balanced_field.html).
+  Illustrates the full multi-phase BFL calculation using Dymos and OpenMDAO.
+
+- **Betts, J. T. (2010)**. *Practical Methods for Optimal Control and Estimation Using Nonlinear Programming*. SIAM.  
+  Discusses takeoff and landing trajectory optimization in Chapter 6.

ext/Descriptions/brachistochrone.md

@@ -0,0 +1,55 @@
+The **Brachistochrone problem** is a classical benchmark in the history of calculus of variations and optimal control.  
+It seeks the shape of a curve (or wire) connecting two points $A$ and $B$ within a vertical plane, such that a bead sliding frictionlessly under the influence of uniform gravity travels from $A$ to $B$ in the shortest possible time.  
+Originating from the challenge posed by Johann Bernoulli in 1696 [Bernoulli 1696](https://doi.org/10.1002/andp.19163551307), it marks the birth of modern optimal control theory.  
+The problem involves nonlinear dynamics where the state includes position and velocity, and the control is the path's angle, with the final time $t_f$ being a decision variable to be minimised.
+
+### Mathematical formulation
+
+The problem can be stated as a free final time optimal control problem:
+
+```math
+\begin{aligned}
+\min_{u(\cdot), t_f} \quad & J = t_f = \int_0^{t_f} 1 \,\mathrm{d}t \\[1em]
+\text{s.t.} \quad & \dot{x}(t) = v(t) \sin u(t), \\[0.5em]
+& \dot{y}(t) = v(t) \cos u(t), \\[0.5em]
+& \dot{v}(t) = g \cos u(t), \\[0.5em]
+& x(0) = x_0, \; y(0) = y_0, \; v(0) = v_0, \\[0.5em]
+& x(t_f) = x_f, \; y(t_f) = y_f,
+\end{aligned}
+```
+
+where $u(t)$ represents the angle of the tangent to the curve with respect to the vertical axis, and $g$ is the gravitational acceleration.
+
+### Qualitative behaviour
+
+This problem is a classic example of **minimum time control** with nonlinear dynamics.  
+Unlike the shortest path problem (a straight line), the optimal trajectory balances the need to minimize path length with the need to maximize velocity.
+
+The analytical solution to this problem is a **cycloid**—the curve traced by a point on the rim of a circular wheel as the wheel rolls along a straight line without slipping.  
+Specifically:
+
+- **Energy Conservation**: Since the system is conservative and frictionless, the speed at any height $h$ is given by $v = \sqrt{2gh}$ (assuming start from rest). This implies that maximizing vertical drop early in the trajectory increases velocity for the remainder of the path.
+- **Concavity**: The optimal curve is concave up. It starts steeply (vertical tangent if $v_0=0$) to gain speed quickly, then flattens out near the bottom before potentially rising again to reach the target.
+- **Singularity**: At the initial point, if $v_0=0$, the control is theoretically singular as the bead has no velocity to direct. Numerical solvers often require a small non-zero initial velocity or a robust initial guess to handle this start.
+
+The control $u(t)$ is continuous and varies smoothly to trace the cycloidal arc.
+
+### Characteristics
+
+- **Nonlinear dynamics** involving trigonometric functions of the control.
+- **Free final time** problem (time-optimal).
+- Analytical solution is a **Cycloid**.
+- Historically significant as the first problem solved using techniques that evolved into the Calculus of Variations.
+
+### References
+
+- Bernoulli, J. (1696). *Problema novum ad cujus solutionem Mathematici invitantur*.  
+  Acta Eruditorum.  
+  The original publication posing the challenge to the mathematicians of Europe, solved by Newton, Leibniz, L'Hôpital, and the Bernoulli brothers.
+
+- Sussmann, H. J., & Willems, J. C. (1997). *300 years of optimal control: from the brachystochrone to the maximum principle*.  
+  IEEE Control Systems Magazine. [doi.org/10.1109/37.588108](https://doi.org/10.1109/37.588108)  
+  A comprehensive historical review linking the classical Brachistochrone problem to modern optimal control theory and Pontryagin's Maximum Principle.
+
+- Dymos Examples: Brachistochrone. [OpenMDAO/Dymos](https://openmdao.github.io/dymos/examples/brachistochrone/brachistochrone.html)  
+  The numerical implementation serving as the basis for the current benchmark formulation.

ext/Descriptions/bryson_denham.md

@@ -0,0 +1,56 @@
+The **Bryson–Denham problem** is a classic optimal control benchmark involving a state inequality constraint.  
+It models the movement of a unit mass on a frictionless plane, where the goal is to relocate the mass over a fixed time interval with minimum control effort, while ensuring the position never exceeds a certain limit.  
+Originally introduced in [Bryson et al. 1963](https://doi.org/10.2514/3.2107), it is a standard test case for evaluating the ability of numerical solvers to handle state-constrained trajectories.
+
+### Mathematical formulation
+
+The problem can be stated as
+
+```math
+$$
+\begin{aligned}
+\min_{x,u} \quad & J(x,u) = \int_0^1 \frac{1}{2} u^2(t) \,\mathrm{d}t \\[1em]
+\text{s.t.} \quad & \dot{x}_1(t) = x_2(t), \quad \dot{x}_2(t) = u(t), \\[0.5em]
+& x_1(0) = 0, \; x_1(1) = 0, \; x_2(0) = 1, \; x_2(1) = -1, \\[0.5em]
+& x_1(t) \le a.
+\end{aligned}
+$$
+```
+
+where $x_1$ represents position, $x_2$ velocity, $u$ acceleration (control), and $a$ is the maximum allowable displacement (e.g., $a=1/9$).
+
+### Qualitative behaviour
+
+The problem features a **second-order state constraint** because the control $u$ only appears in the second derivative of the constrained state $x_1$:
+
+```math
+$$
+\ddot{x}_1(t) = u(t).
+$$
+```
+
+When the trajectory is on a **boundary arc** ($x_1(t) = a$), the derivatives $\dot{x}_1$ and $\ddot{x}_1$ must be zero. This implies that while the constraint is active, the velocity $x_2$ is zero and the control $u$ is zero:
+
+```math
+$$
+u(t) = 0.
+$$
+```
+
+The solution structure changes based on the value of the constraint parameter $a$:
+
+* **Unconstrained Case**: If $a$ is sufficiently large ($a \ge 1/4$), the state constraint is never active.
+* **Touch Point**: For a critical value of $a$, the trajectory just touches the boundary $x_1 = a$ at a single point $t=0.5$.
+* **Boundary Arc**: For smaller values (e.g., $a=1/9$), the trajectory remains on the boundary for a finite time interval. During this interval, the control vanishes.
+
+### Characteristics
+
+* Linear–quadratic dynamics with a second-order state inequality constraint.  
+* Demonstrates a "bang-off-bang" control structure where the control is zero on the constraint boundary.  
+* Used to test the accuracy of state constraint handling and the detection of switching times.
+
+### References
+
+* 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)  
+* Dymos Documentation: Bryson-Denham Problem. [dymos/examples/bryson_denham](https://openmdao.github.io/dymos/examples/bryson_denham/bryson_denham.html)

ext/Descriptions/mountain_car.md

@@ -0,0 +1,57 @@
+The **Mountain Car problem** is a classic benchmark in control and reinforcement learning.  
+It involves an underpowered car that must climb a steep hill to reach a target position.  
+Because the car's engine is too weak to climb the hill directly, it must first drive away from the goal to gain momentum by oscillating in a valley.  
+The goal is to reach the target position in **minimum time**.
+
+### Mathematical formulation
+
+The problem can be stated as
+
+```math
+\begin{aligned}
+\min_{x,v,u,t_f} \quad & J = t_f \[0.5em]
+	ext{s.t.} \quad & \dot{x}(t) = v(t), \[0.5em]
+& \dot{v}(t) = a \, u(t) - b \, \cos(c \, x(t)), \[0.5em]
+& x(0) = -0.5, \quad v(0) = 0.0, \[0.5em]
+& x(t_f) = 0.5, \quad v(t_f) \ge 0.0, \[0.5em]
+& -1.2 \le x(t) \le 0.5, \[0.5em]
+& -0.07 \le v(t) \le 0.07, \[0.5em]
+& -1 \le u(t) \le 1, \[0.5em]
+& t_f \ge 0.
+\end{aligned}
+```
+
+### Parameters
+
+| Parameter | Symbol | Value |
+|-----------|--------|-------|
+| Power coefficient | $a$ | 0.001 |
+| Gravity coefficient | $b$ | 0.0025 |
+| Frequency coefficient | $c$ | 3.0 |
+| Initial position | $x_0$ | -0.5 |
+| Target position | $x_f$ | 0.5 |
+
+### Qualitative behaviour
+
+- The car must oscillate back and forth in the valley to build enough kinetic energy to climb the right hill.
+- The optimal strategy typically involves a sequence of full-throttle accelerations in alternating directions (**bang-bang control**).
+- The minimum time objective makes the problem sensitive to the initial guess for the final time.
+
+### Characteristics
+
+- Two-dimensional state space (position and velocity),
+- Scalar control (effort/acceleration),
+- **Free final time** with minimum time objective,
+- State constraints (position and velocity bounds),
+- Control constraints (bounds on acceleration).
+
+### References
+
+- **Dymos Documentation**. [*The Mountain Car Problem*](https://openmdao.github.io/dymos/examples/mountain_car/mountain_car.html).
+  Provides a detailed description and implementation of the Mountain Car problem using the Dymos library.
+
+- **OpenAI Gym**. [*MountainCar-v0*](https://gym.openai.com/envs/MountainCar-v0/).
+  A standard reinforcement learning environment for the Mountain Car problem.
+
+- Moore, A. W. (1990). *Efficient Memory-based Learning for Robot Control*. PhD thesis, University of Cambridge.  
+  Introduces the Mountain Car problem as a benchmark for reinforcement learning and control.

ext/JuMPModels/balanced_field.jl

@@ -0,0 +1,130 @@
+"""
+$(TYPEDSIGNATURES)
+
+Constructs and returns a JuMP model for the **Aircraft Balanced Field Length Calculation (Takeoff phase)**.  
+The model represents the dynamics of the aircraft during takeoff and seeks to minimise the final range.
+
+# Arguments
+
+- `::JuMPBackend`: Specifies the backend for building the JuMP model.
+- `grid_size::Int=200`: (Keyword) Number of discretisation steps for the time horizon.
+
+# Returns
+
+- `model::JuMP.Model`: A JuMP model representing the problem.
+"""
+function OptimalControlProblems.balanced_field(
+    ::JuMPBackend,
+    args...;
+    grid_size::Int=grid_size_data(:balanced_field),
+    parameters::Union{Nothing,NamedTuple}=nothing,
+    kwargs...,
+)
+
+    # parameters
+    params = parameters_data(:balanced_field, parameters)
+    t0 = params[:t0]
+    m = params[:m]
+    g = params[:g]
+    ρ = params[:ρ]
+    S = params[:S]
+    CD0 = params[:CD0]
+    AR = params[:AR]
+    e = params[:e]
+    CL0 = params[:CL0]
+    CL_max = params[:CL_max]
+    h_w = params[:h_w]
+    span = params[:span]
+    α_max = params[:α_max]
+    T = params[:T]
+
+    r_t0 = params[:r_t0]
+    v_t0 = params[:v_t0]
+    h_t0 = params[:h_t0]
+    γ_t0 = params[:γ_t0]
+    h_tf = params[:h_tf]
+    γ_tf = params[:γ_tf]
+
+    α_min = params[:α_min]
+    α_max_ctrl = params[:α_max_ctrl]
+
+    # Constants
+    b = span / 2.0
+    K_nom = 1.0 / (pi * AR * e)
+
+    # model
+    model = JuMP.Model(args...; kwargs...)
+
+    # metadata
+    model[:time_grid] = () -> range(t0, value(model[:tf]), grid_size + 1)
+    model[:state_components] = ["r", "v", "h", "γ"]
+    model[:costate_components] = ["∂r", "∂v", "∂h", "∂γ"]
+    model[:control_components] = ["α"]
+    model[:variable_components] = ["tf"]
+
+    @expression(model, N, grid_size)
+
+    # variables
+    @variables(
+        model,
+        begin
+            r[0:N], (start = r_t0)
+            v[0:N] ≥ 1.0, (start = v_t0 + 10.0)
+            h[0:N] ≥ 0.0, (start = h_tf / 2.0)
+            γ[0:N], (start = γ_tf)
+            0 ≤ α[0:N] ≤ α_max_ctrl, (start = 0.1)
+            tf ≥ 0.1, (start = 10.0)
+        end
+    )
+
+    # boundary constraints
+    @constraints(
+        model,
+        begin
+            r[0] == r_t0
+            v[0] == v_t0
+            h[0] == h_t0
+            γ[0] == γ_t0
+
+            h[N] == h_tf
+            γ[N] == γ_tf
+        end
+    )
+
+    # dynamics
+    @expressions(
+        model,
+        begin
+            Δt, (tf - t0) / N
+
+            q[i = 0:N], 0.5 * ρ * v[i]^2
+            CL[i = 0:N], CL0 + (α[i] / α_max) * (CL_max - CL0)
+            L[i = 0:N], q[i] * S * CL[i]
+
+            h_eff[i = 0:N], h[i] + h_w
+            term_h[i = 0:N], 33.0 * abs(h_eff[i] / b)^1.5
+            K[i = 0:N], K_nom * term_h[i] / (1.0 + term_h[i])
+            D[i = 0:N], q[i] * S * (CD0 + K[i] * CL[i]^2)
+
+            dr[i = 0:N], v[i] * cos(γ[i])
+            dv[i = 0:N], (T * cos(α[i]) - D[i]) / m - g * sin(γ[i])
+            dh[i = 0:N], v[i] * sin(γ[i])
+            dγ[i = 0:N], (T * sin(α[i]) + L[i]) / (m * v[i]) - (g * cos(γ[i])) / v[i]
+        end
+    )
+
+    @constraints(
+        model,
+        begin
+            ∂r[i = 1:N], r[i] == r[i - 1] + 0.5 * Δt * (dr[i] + dr[i - 1])
+            ∂v[i = 1:N], v[i] == v[i - 1] + 0.5 * Δt * (dv[i] + dv[i - 1])
+            ∂h[i = 1:N], h[i] == h[i - 1] + 0.5 * Δt * (dh[i] + dh[i - 1])
+            ∂γ[i = 1:N], γ[i] == γ[i - 1] + 0.5 * Δt * (dγ[i] + dγ[i - 1])
+        end
+    )
+
+    # objective
+    @objective(model, Min, r[N])
+
+    return model
+end