test implement balanced field

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/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], (start = v_t0 + 10.0)
+            h[0:N], (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 * (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

ext/MetaData/balanced_field.jl

@@ -0,0 +1,29 @@
+balanced_field_meta = OrderedDict(
+    :grid_size => 200,
+    :parameters => (
+        t0 = 0.0,
+        m = 79015.8,      # mass (kg)
+        g = 9.80665,      # gravity (m/s^2)
+        ρ = 1.225,        # air density (kg/m^3)
+        S = 124.7,        # surface (m^2)
+        CD0 = 0.03,       # zero-lift drag coefficient
+        AR = 9.45,        # aspect ratio
+        e = 0.801,        # Oswald efficiency
+        CL0 = 0.5,        # zero-alpha lift coefficient
+        CL_max = 2.0,     # max lift coefficient
+        h_w = 1.0,        # height of wing above CoG (m)
+        span = 35.7,      # wingspan (m)
+        α_max = 0.1745,   # alpha at CL_max (10 degrees in rad)
+        T = 120102.0,     # thrust (N) - one engine out
+        r_t0 = 600.0,     # initial range (m) at V1
+        v_t0 = 70.0,      # initial velocity (m/s) at V1
+        h_t0 = 0.0,       # initial altitude (m)
+        γ_t0 = 0.0,       # initial flight path angle (rad)
+        h_tf = 10.668,    # final altitude (35 ft in m)
+        γ_tf = 0.0873,    # final flight path angle (5 degrees in rad)
+        v_min = 0.0,
+        v_max = 100.0,
+        α_min = -0.1745,  # -10 degrees
+        α_max_ctrl = 0.2618, # 15 degrees
+    ),
+)

ext/OptimalControlModels/balanced_field.jl

@@ -0,0 +1,126 @@
+"""
+$(TYPEDSIGNATURES)
+
+Constructs an **OptimalControl problem** for the Aircraft Balanced Field Length Calculation (Takeoff phase).  
+The objective is to minimise the final range to reach a specified altitude (35 ft) with one engine out.  
+The problem formulation is based on the Dymos balanced field example.
+
+# Arguments
+
+- `::OptimalControlBackend`: Placeholder type specifying the OptimalControl backend or solver interface.
+- `grid_size::Int=200`: (Keyword) Number of discretisation points for the direct transcription grid.
+
+# Returns
+
+- `docp`: The direct optimal control problem object representing the problem.
+- `nlp`: The corresponding nonlinear programming model.
+
+# Example
+
+```julia-repl
+julia> using OptimalControlProblems
+julia> docp = OptimalControlProblems.balanced_field(OptimalControlBackend());
+```
+"""
+function OptimalControlProblems.balanced_field(
+    ::OptimalControlBackend,
+    description::Symbol...;
+    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
+    ocp = @def begin
+        tf ∈ R, variable
+        t ∈ [t0, tf], time
+        x = (r, v, h, γ) ∈ R⁴, state
+        α ∈ R, control
+
+        r(t0) == r_t0
+        v(t0) == v_t0
+        h(t0) == h_t0
+        γ(t0) == γ_t0
+
+        h(tf) == h_tf
+        γ(tf) == γ_tf
+
+        0 ≤ α(t) ≤ α_max_ctrl # alpha
+        tf ≥ 0.1
+
+        ẋ(t) == dynamics(x(t), α(t), m, g, ρ, S, CD0, CL0, CL_max, α_max, T, b, K_nom, h_w)
+
+        r(tf) → min
+    end
+
+    function dynamics(x, α, m, g, ρ, S, CD0, CL0, CL_max, α_max, T, b, K_nom, h_w)
+        r, v, h, γ = x
+
+        q = 0.5 * ρ * v^2
+        CL = CL0 + (α / α_max) * (CL_max - CL0)
+        L = q * S * CL
+
+        h_eff = h + h_w
+        term_h = 33.0 * (h_eff / b)^1.5
+        K = K_nom * term_h / (1.0 + term_h)
+        D = q * S * (CD0 + K * CL^2)
+
+        rdot = v * cos(γ)
+        vdot = (T * cos(α) - D) / m - g * sin(γ)
+        hdot = v * sin(γ)
+        γdot = (T * sin(α) + L) / (m * v) - (g * cos(γ)) / v
+
+        return [rdot, vdot, hdot, γdot]
+    end
+
+    # initial guess
+    tf_guess = 10.0
+    xinit = [r_t0, v_t0 + 10.0, h_tf / 2.0, γ_tf]
+    uinit = [0.1]
+    vinit = [tf_guess]
+    init = (state=xinit, control=uinit, variable=vinit)
+
+    # discretise
+    docp = direct_transcription(
+        ocp,
+        description...;
+        lagrange_to_mayer=false,
+        init=init,
+        grid_size=grid_size,
+        disc_method=:trapeze,
+        kwargs...,
+    )
+
+    return docp
+end