add brachistochrone

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/JuMPModels/brachistochrone.jl

@@ -0,0 +1,128 @@
+"""
+$(TYPEDSIGNATURES)
+
+Constructs a JuMP model representing the **Brachistochrone optimal control problem**.
+The objective is to minimize the final time `tf` to travel between two points under gravity.
+The problem uses a direct transcription (trapezoidal rule) manually implemented in JuMP.
+
+# Arguments
+
+- `::JuMPBackend`: Placeholder type to specify the JuMP backend or solver interface.
+- `grid_size::Int=50`: (Keyword) Number of discretisation steps for the time grid.
+
+# Returns
+
+- `model::JuMP.Model`: A JuMP model containing the decision variables (including final time), dynamics constraints, and boundary conditions.
+
+# Example
+
+```julia-repl
+julia> using OptimalControlProblems
+julia> using JuMP
+julia> model = OptimalControlProblems.brachistochrone(JuMPBackend(); N=100)
+```
+
+# References
+
+- Dymos Brachistochrone: [https://openmdao.github.io/dymos/examples/brachistochrone/brachistochrone.html](https://openmdao.github.io/dymos/examples/brachistochrone/brachistochrone.html)
+"""
+function OptimalControlProblems.brachistochrone(
+    ::JuMPBackend,
+    args...;
+    grid_size::Int=grid_size_data(:brachistochrone),
+    parameters::Union{Nothing,NamedTuple}=nothing,
+    kwargs...,
+)
+
+    # parameters
+    params = parameters_data(:brachistochrone, parameters)
+    g  = params[:g]
+    t0 = params[:t0]
+    x0 = params[:x0]
+    y0 = params[:y0]
+    v0 = params[:v0]
+    xf = params[:xf]
+    yf = params[:yf]
+
+    # model
+    model = JuMP.Model(args...; kwargs...)
+
+    # N = grid_size
+    @expression(model, N, grid_size)
+
+    @variables(
+        model,
+        begin
+            0.1 <= tf <= 20.0, (start = 2.0)
+            px[0:N]
+            py[0:N]
+            v[0:N]
+            u[0:N]
+        end
+    )
+
+
+    model[:time_grid] = () -> range(t0, value(tf), N+1)
+    model[:state_components] = ["px", "py", "v"]
+    model[:costate_components] = ["∂px", "∂py", "∂v"]
+    model[:control_components] = ["u"]
+    model[:variable_components] = ["tf"]
+
+    for i in 0:N
+        alpha = i / N
+        set_start_value(px[i], x0 + alpha * (xf - x0))
+        set_start_value(py[i], y0 + alpha * (yf - y0))
+        set_start_value(v[i], v0 + alpha * 10.0) # Estimate speed
+        set_start_value(u[i], 1.57) # ~90 degrees
+    end
+
+    # boundary constraints
+    @constraints(
+        model,
+        begin
+            # Start
+            px[0] == x0
+            py[0] == y0
+            v[0] == v0
+
+            # End
+            px[N] == xf
+            py[N] == yf
+            # v[N] is free
+        end
+    )
+
+    # dynamics
+    @expressions(
+        model,
+        begin
+            # Time step is variable: dt = (tf - t0) / N
+            dt, (tf - t0) / N
+
+            # Dynamics expressions (Dymos formulation)
+            # dpx/dt = v * sin(u)
+            dpx[i = 0:N], v[i] * sin(u[i])
+
+            # dpy/dt = -v * cos(u)
+            dpy[i = 0:N], -v[i] * cos(u[i])
+
+            # dv/dt = g * cos(u)
+            dv[i = 0:N], g * cos(u[i])
+        end
+    )
+
+    # Trapezoidal rule integration (Implicit)
+    @constraints(
+        model,
+        begin
+            ∂px[i = 1:N], px[i] == px[i - 1] + 0.5 * dt * (dpx[i] + dpx[i - 1])
+            ∂py[i = 1:N], py[i] == py[i - 1] + 0.5 * dt * (dpy[i] + dpy[i - 1])
+            ∂v[i = 1:N], v[i] == v[i - 1] + 0.5 * dt * (dv[i] + dv[i - 1])
+        end
+    )
+
+    # objective: Minimize final time
+    @objective(model, Min, tf)
+
+    return model
+end

ext/MetaData/brachistochrone.jl

@@ -0,0 +1,12 @@
+brachistochrone_meta = OrderedDict(
+    :grid_size => 500,                 
+    :parameters => (
+        g  = 9.80665,                 
+        t0 = 0.0,                      
+        x0 = 0.0,                      
+        y0 = 10.0,                     
+        v0 = 0.0,                    
+        xf = 10.0,                    
+        yf = 5.0,                     
+    ),
+)

ext/OptimalControlModels/brachistochrone.jl

@@ -0,0 +1,93 @@
+"""
+$(TYPEDSIGNATURES)
+
+Constructs an **OptimalControl problem** representing the Brachistochrone problem using the OptimalControl backend.
+The function sets up the state and control variables, boundary conditions, dynamics, and the objective functional.
+It then performs direct transcription to generate a discrete optimal control problem (DOCP).
+
+# Arguments
+
+- `::OptimalControlBackend`: Placeholder type to specify the OptimalControl backend or solver interface.
+- `grid_size::Int=grid_size_data(:brachistochrone)`: (Keyword) Number of discretisation points for the direct transcription grid.
+- `parameters::Union{Nothing,NamedTuple}=nothing`: (Keyword) Custom parameters to override defaults.
+
+# Returns
+
+- `docp`: The direct optimal control problem object, representing the discretised problem.
+
+# Example
+```julia-repl
+julia> using OptimalControlProblems
+
+julia> docp = OptimalControlProblems.brachistochrone(OptimalControlBackend(); N=100);
+```
+
+# References
+
+- Dymos Brachistochrone: https://openmdao.github.io/dymos/examples/brachistochrone/brachistochrone.html
+"""
+function OptimalControlProblems.brachistochrone(
+    ::OptimalControlBackend,
+    description::Symbol...;
+    grid_size::Int=grid_size_data(:brachistochrone),
+    parameters::Union{Nothing,NamedTuple}=nothing,
+    kwargs...,
+)
+
+    params = parameters_data(:brachistochrone, parameters)
+    g  = params[:g]
+    t0 = params[:t0]
+    x0 = params[:x0]
+    y0 = params[:y0]
+    v0 = params[:v0]
+    xf = params[:xf]
+    yf = params[:yf]
+
+    ocp = @def begin
+        tf ∈ R, variable
+        t ∈ [t0, tf], time
+
+        z = (px, py, v) ∈ R³, state
+        u ∈ R, control
+
+        0.1 ≤ tf ≤ 20.0 
+
+        px(t0) == x0
+        py(t0) == y0
+        v(t0)  == v0
+
+        px(tf) == xf
+        py(tf) == yf
+
+        ∂(px)(t) == v(t) * sin(u(t))
+        ∂(py)(t) == -v(t) * cos(u(t))
+        ∂(v)(t)  == g * cos(u(t))
+
+        # Objectif
+        tf → min
+    end
+
+    # initial guess: linear interpolation to match JuMP
+    tf_guess = 2.0
+    init = (
+        state = t -> [
+            x0 + (t - t0) / (tf_guess - t0) * (xf - x0),
+            y0 + (t - t0) / (tf_guess - t0) * (yf - y0),
+            v0 + (t - t0) / (tf_guess - t0) * 10.0
+        ],
+        control = 1.57,
+        variable = tf_guess
+    )
+
+    docp = direct_transcription(
+        ocp,
+        description...;
+        lagrange_to_mayer=false,
+        init = init,
+        grid_size = grid_size,
+        disc_method=:trapeze,
+        kwargs...
+    )
+
+    return docp
+end

ext/OptimalControlModels_s/brachistochrone_s.jl

@@ -0,0 +1,96 @@
+"""
+$(TYPEDSIGNATURES)
+
+Constructs an **OptimalControl problem** representing the Brachistochrone problem using the OptimalControl backend.
+The function sets up the state and control variables, boundary conditions, dynamics, and the objective functional.
+It then performs direct transcription to generate a discrete optimal control problem (DOCP).
+
+# Arguments
+
+- `::OptimalControlBackend`: Placeholder type to specify the OptimalControl backend or solver interface.
+- `grid_size::Int=grid_size_data(:brachistochrone)`: (Keyword) Number of discretisation points for the direct transcription grid.
+- `parameters::Union{Nothing,NamedTuple}=nothing`: (Keyword) Custom parameters to override defaults.
+
+# Returns
+
+- `docp`: The direct optimal control problem object, representing the discretised problem.
+
+# Example
+```julia-repl
+julia> using OptimalControlProblems
+
+julia> docp = OptimalControlProblems.brachistochrone(OptimalControlBackend(); N=100);
+```
+
+# References
+
+- Dymos Brachistochrone: https://openmdao.github.io/dymos/examples/brachistochrone/brachistochrone.html
+"""
+function OptimalControlProblems.brachistochrone_s(
+    ::OptimalControlBackend,
+    description::Symbol...;
+    grid_size::Int=grid_size_data(:brachistochrone),
+    parameters::Union{Nothing,NamedTuple}=nothing,
+    kwargs...,
+)
+
+    params = parameters_data(:brachistochrone, parameters)
+    g  = params[:g]
+    t0 = params[:t0]
+    x0 = params[:x0]
+    y0 = params[:y0]
+    v0 = params[:v0]
+    xf = params[:xf]
+    yf = params[:yf]
+
+    # model
+    ocp = @def begin
+
+        tf ∈ R, variable
+        t ∈ [t0, tf], time
+
+
+        z = (px, py, v) ∈ R³, state
+        u ∈ R, control
+
+
+        z(t0) == [x0, y0, v0]
+
+
+        px(tf) == xf
+        py(tf) == yf
+
+        0.1 ≤ tf ≤ 20.0
+
+        ∂(px)(t) == v(t) * sin(u(t))
+        ∂(py)(t) == -v(t) * cos(u(t))
+        ∂(v)(t)  == g * cos(u(t))
+
+        tf → min
+    end
+
+    # initial guess: linear interpolation to match JuMP
+    tf_guess = 2.0
+    init = (
+        state = t -> [
+            x0 + (t - t0) / (tf_guess - t0) * (xf - x0),
+            y0 + (t - t0) / (tf_guess - t0) * (yf - y0),
+            v0 + (t - t0) / (tf_guess - t0) * 10.0
+        ],
+        control = 1.57, 
+        variable = tf_guess
+    )
+
+    # discretise the optimal control problem
+    docp = direct_transcription(
+        ocp,
+        description...;
+        lagrange_to_mayer=false,
+        init=init,
+        grid_size=grid_size,
+        disc_method=:trapeze,
+        kwargs...,
+    )
+
+    return docp
+end