diff --git a/ext/Descriptions/brachistochrone.md b/ext/Descriptions/brachistochrone.md new file mode 100644 index 00000000..da388cf9 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/ext/JuMPModels/brachistochrone.jl b/ext/JuMPModels/brachistochrone.jl new file mode 100644 index 00000000..71319d34 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/ext/MetaData/brachistochrone.jl b/ext/MetaData/brachistochrone.jl new file mode 100644 index 00000000..17d700d0 --- /dev/null +++ b/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, + ), +) \ No newline at end of file diff --git a/ext/OptimalControlModels/brachistochrone.jl b/ext/OptimalControlModels/brachistochrone.jl new file mode 100644 index 00000000..860ea7eb --- /dev/null +++ b/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 \ No newline at end of file diff --git a/ext/OptimalControlModels_s/brachistochrone_s.jl b/ext/OptimalControlModels_s/brachistochrone_s.jl new file mode 100644 index 00000000..8d46e45f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/OptimalControlProblems.jl b/src/OptimalControlProblems.jl index 9729407d..2a27ad5c 100644 --- a/src/OptimalControlProblems.jl +++ b/src/OptimalControlProblems.jl @@ -72,6 +72,17 @@ function make_list_of_problems() :robot, :space_shuttle, :steering, + :beam, + :cart_pendulum, + :chain, + :dielectrophoretic_particle, + :double_oscillator, + :electric_vehicle, + :glider, + :jackson, + :robbins, + :rocket, + :vanderpol, ] list_of_problems = setdiff(list_of_problems, problems_to_exclude)