Hang Glider Flight#

Mathematical Formulation#

Problem Statement#

Find the optimal lift coefficient \(C_L(t)\) that maximizes the final range:

(1)#\[J = -x(t_f)\]

Subject to the atmospheric flight dynamics:

\[\frac{dx}{dt} = v_x\]
\[\frac{dy}{dt} = v_y\]
\[\frac{dv_x}{dt} = -\frac{1}{m}(L \sin\eta + D \cos\eta)\]
\[\frac{dv_y}{dt} = \frac{1}{m}(L \cos\eta - D \sin\eta) - g\]

Where the aerodynamic forces are:

\[L = \frac{1}{2} C_L \rho S v_r^2\]
\[D = \frac{1}{2} C_D \rho S v_r^2\]
\[C_D = C_0 + K C_L^2\]

And the flight path parameters:

\[v_r = \sqrt{v_x^2 + V_y^2}\]
\[V_y = v_y - u_a(x)\]
\[\sin\eta = \frac{V_y}{v_r}, \quad \cos\eta = \frac{v_x}{v_r}\]

Thermal Updraft Model#

\[u_a(x) = U_m \left(1 - X\right) e^{-X}\]

where \(X = (\frac{x}{R} - 2.5)^2\)

Boundary Conditions#

  • Initial conditions: \(x(0) = 0\), \(y(0) = 1000\), \(v_x(0) = 13.2275675\), \(v_y(0) = -1.28750052\)

  • Final conditions: \(x(t_f) = \text{free}\), \(y(t_f) = 900\), \(v_x(t_f) = 13.2275675\), \(v_y(t_f) = -1.28750052\)

  • Control bounds: \(0 \leq C_L \leq 1.4\)

  • State bounds: \(0 \leq x \leq 1500\) m, \(0 \leq y \leq 1100\) m, \(0 \leq v_x \leq 15\) m/s, \(-4 \leq v_y \leq 4\) m/s

Physical Parameters#

  • Mass: \(m = 100\) kg

  • Gravitational acceleration: \(g = 9.80665\) m/s²

  • Maximum updraft velocity: \(U_m = 2.5\) m/s

  • Thermal radius: \(R = 100\) m

  • Parasitic drag coefficient: \(C_0 = 0.034\)

  • Induced drag factor: \(K = 0.069662\)

  • Wing area: \(S = 14\) m²

  • Air density: \(\rho = 1.13\) kg/m³

State Variables#

  • \(x(t)\): Horizontal position (m)

  • \(y(t)\): Altitude (m)

  • \(v_x(t)\): Horizontal velocity (m/s)

  • \(v_y(t)\): Vertical velocity (m/s)

Control Variable#

  • \(C_L(t)\): Lift coefficient

Notes#

This problem demonstrates atmospheric flight optimization where the hang glider must exploit thermal updrafts to maximize range while satisfying altitude and velocity constraints.

Running This Example#

cd examples/hang_glider
python hang_glider.py

Code Implementation#

examples/hang_glider/hang_glider.py#
  1import casadi as ca
  2import numpy as np
  3
  4import maptor as mtor
  5
  6
  7# ============================================================================
  8# Physical Constants and Parameters
  9# ============================================================================
 10
 11# Physical constants
 12M = 100.0
 13G = 9.80665
 14U_M = 2.5
 15R = 100.0
 16C0 = 0.034
 17K = 0.069662
 18S = 14.0
 19RHO = 1.13
 20
 21
 22# ============================================================================
 23# Scaling Factors
 24# ============================================================================
 25
 26
 27X_SCALE = 1000.0  # position scaling (km)
 28V_SCALE = 10.0  # velocity scaling (10 m/s)
 29
 30
 31# ============================================================================
 32# Initial and Target Conditions
 33# ============================================================================
 34
 35# Scaled values
 36X0_SCALED = 0.0 / X_SCALE
 37Y0_SCALED = 1000.0 / X_SCALE
 38YF_SCALED = 900.0 / X_SCALE
 39VX0_SCALED = 13.2275675 / V_SCALE
 40VY0_SCALED = -1.28750052 / V_SCALE
 41VXF_SCALED = 13.2275675 / V_SCALE
 42VYF_SCALED = -1.28750052 / V_SCALE
 43
 44# Scaled bounds
 45X_MAX_SCALED = 1500.0 / X_SCALE
 46Y_MAX_SCALED = 1100.0 / X_SCALE
 47VX_MAX_SCALED = 15.0 / V_SCALE
 48VY_MIN_SCALED = -4.0 / V_SCALE
 49VY_MAX_SCALED = 4.0 / V_SCALE
 50
 51
 52# ============================================================================
 53# Problem Setup
 54# ============================================================================
 55
 56problem = mtor.Problem("Hang Glider")
 57phase = problem.set_phase(1)
 58
 59
 60# ============================================================================
 61# Variables
 62# ============================================================================
 63
 64t = phase.time(initial=0.0, final=(0.1, 200.0))
 65x_s = phase.state("x_scaled", initial=X0_SCALED, boundary=(0.0, X_MAX_SCALED))
 66y_s = phase.state("y_scaled", initial=Y0_SCALED, final=YF_SCALED, boundary=(0.0, Y_MAX_SCALED))
 67vx_s = phase.state("vx_scaled", initial=VX0_SCALED, final=VXF_SCALED, boundary=(0.0, VX_MAX_SCALED))
 68vy_s = phase.state(
 69    "vy_scaled", initial=VY0_SCALED, final=VYF_SCALED, boundary=(VY_MIN_SCALED, VY_MAX_SCALED)
 70)
 71CL = phase.control("CL", boundary=(0.0, 1.4))
 72
 73
 74# ============================================================================
 75# Dynamics
 76# ============================================================================
 77
 78
 79def _define_scaled_dynamics():
 80    # Convert scaled variables to physical units
 81    x_phys = x_s * X_SCALE
 82    y_s * X_SCALE
 83    vx_phys = vx_s * V_SCALE
 84    vy_phys = vy_s * V_SCALE
 85
 86    # Numerical safeguards
 87    vr_squared = vx_phys * vx_phys + vy_phys * vy_phys
 88    vr = ca.sqrt(ca.fmax(vr_squared, 1e-6))
 89
 90    # Aerodynamic model
 91    CD = C0 + K * CL * CL
 92    D = 0.5 * CD * RHO * S * vr * vr
 93    L = 0.5 * CL * RHO * S * vr * vr
 94
 95    # Updraft model with numerical bounds
 96    X_arg = x_phys / R - 2.5
 97    X = X_arg * X_arg
 98    X_safe = ca.fmin(X, 50.0)
 99    ua = U_M * (1.0 - X_safe) * ca.exp(-X_safe)
100
101    Vy = vy_phys - ua
102    sin_eta = Vy / vr
103    cos_eta = vx_phys / vr
104    W = M * G
105
106    # Dynamics in physical units
107    x_dot_phys = vx_phys
108    y_dot_phys = vy_phys
109    vx_dot_phys = (1.0 / M) * (-L * sin_eta - D * cos_eta)
110    vy_dot_phys = (1.0 / M) * (L * cos_eta - D * sin_eta - W)
111
112    # Manual scaling of derivatives (NO TIME SCALING)
113    x_dot_scaled = x_dot_phys / X_SCALE
114    y_dot_scaled = y_dot_phys / X_SCALE
115    vx_dot_scaled = vx_dot_phys / V_SCALE
116    vy_dot_scaled = vy_dot_phys / V_SCALE
117
118    return {
119        x_s: x_dot_scaled,
120        y_s: y_dot_scaled,
121        vx_s: vx_dot_scaled,
122        vy_s: vy_dot_scaled,
123    }
124
125
126phase.dynamics(_define_scaled_dynamics())
127
128
129# ============================================================================
130# Objective
131# ============================================================================
132
133# Objective: maximize x(tf)
134problem.minimize(-x_s.final)
135
136
137# ============================================================================
138# Mesh Configuration and Initial Guess
139# ============================================================================
140
141phase.mesh(
142    [3, 3, 3, 3, 3, 3, 3, 3, 3, 3],  # Low degree polynomials
143    np.linspace(-1.0, 1.0, 11),  # Fine, uniform spacing
144)
145
146# Initial guess
147states_guess = []
148controls_guess = []
149for N in [3, 3, 3, 3, 3, 3, 3, 3, 3, 3]:
150    tau = np.linspace(-1, 1, N + 1)
151    t_norm = (tau + 1) / 2
152
153    x_vals = (0.0 + 1500 * t_norm) / X_SCALE
154    y_vals = (1000.0 + (900.0 - 1000.0) * t_norm) / X_SCALE
155    vx_vals = (13.23 * np.ones(N + 1)) / V_SCALE
156    vy_vals = (-1.288 * np.ones(N + 1)) / V_SCALE
157
158    states_guess.append(np.vstack([x_vals, y_vals, vx_vals, vy_vals]))
159    controls_guess.append(np.ones((1, N)) * 0.7)
160
161phase.guess(
162    states=states_guess,
163    controls=controls_guess,
164    terminal_time=105,
165)
166
167
168# ============================================================================
169# Solve
170# ============================================================================
171
172solution = mtor.solve_adaptive(
173    problem,
174    error_tolerance=1e-6,
175    max_iterations=30,
176    min_polynomial_degree=3,
177    max_polynomial_degree=8,
178    ode_solver_tolerance=1e-5,
179    ode_method="DOP853",
180    nlp_options={
181        "ipopt.max_iter": 2000,
182        "ipopt.tol": 1e-6,
183        "ipopt.constr_viol_tol": 1e-6,
184        "ipopt.acceptable_tol": 1e-3,
185        "ipopt.linear_solver": "mumps",
186        "ipopt.print_level": 0,
187    },
188)
189
190
191# ============================================================================
192# Results Analysis
193# ============================================================================
194
195if solution.status["success"]:
196    final_x_scaled = solution[(1, "x_scaled")][-1]
197    final_x_physical = final_x_scaled * X_SCALE  #
198    flight_time = solution.phases[1]["times"]["final"]
199
200    print(f"MAPTOR result: {final_x_physical:.2f} m in {flight_time:.2f} s")
201    print("Literature ref: 1248.26 m in 98.47 s")
202    print(f"Range error: {abs(final_x_physical - 1248.26):.2f} m")
203    print(f"Time error: {abs(flight_time - 98.47):.2f} s")
204
205    # Final state verification in physical units
206    x_final = solution[(1, "x_scaled")][-1] * X_SCALE
207    y_final = solution[(1, "y_scaled")][-1] * X_SCALE
208    vx_final = solution[(1, "vx_scaled")][-1] * V_SCALE
209    vy_final = solution[(1, "vy_scaled")][-1] * V_SCALE
210
211    print("Final states (physical units):")
212    print(f"  x: {x_final:.2f} m")
213    print(f"  y: {y_final:.2f} m (target: 900.0)")
214    print(f"  vx: {vx_final:.6f} m/s (target: 13.2275675)")
215    print(f"  vy: {vy_final:.6f} m/s (target: -1.28750052)")
216
217    solution.plot()
218else:
219    print(f"Failed: {solution.status['message']}")