2025 Class Project: Minimum time to climb problem¶
The following page contains information about the course project.
The course project must be done by yourself. It is not a group project. You may discuss issues or problems you have, but you cannot share code.
This problem is taken from Betts, Practical Methods for Optimal Control Using Nonlinear Programming, Third Edition. https://epubs.siam.org/doi/book/10.1137/1.9781611976199.
The equations of motion of the aircraft are given as follows:
\begin{equation*} \begin{aligned} \dot{h} & = v \sin \gamma \\ \dot{r} &= v \cos \gamma \\ \dot{v} &= \dfrac{T \cos \alpha - D}{m} - g \sin \gamma \\ \dot{\gamma} &= \dfrac{T}{mv} \sin \alpha + \dfrac{L}{mv} - \dfrac{g \cos \gamma}{v} \\ \dot{m} &= - \dfrac{T}{g I_{sp}} \\ \end{aligned} \end{equation*}
Here \(h\) is the altitude, \(r\) is the down range distance, \(v\) is the velocity, \(\gamma\) is the flight path angle and \(m\) is the aircraft mass.
The control variable is the angle of attack \(\alpha\) which must satisfy the bounds \(-20^{\circ} \le \alpha \le 20^{\circ}\).
To close the equations of motion, we need the trust \(T\), the drag \(D\) and the lift \(L\), as well as the specific impulse, \(I_{sp}\), and the acceleration due to gravity \(g\).
In this problem, we will use a constant specific impulse \(I_{sp} = 1600\, \text{s}\) and a constant acceleration due to gravity of \(g = 32 \text{ft/s}\).
The lift, \(L\), and drag, \(D\), are computed as a function of the velocity, \(v\), altitude, \(h\), angle of attack, \(\alpha\), and Mach number \(M\), as:
\begin{equation*} \begin{aligned} L &= \dfrac{1}{2} \rho(h) v^2 S_{\text{ref}} C_{L}(\alpha, M) \\ D &= \dfrac{1}{2} \rho(h) v^2 S_{\text{ref}} C_{D}(\alpha, M) \end{aligned} \end{equation*}
The drag coefficient \(C_{D}\) and the lift coefficient \(C_{L}\) are computed based on the formula:
\begin{equation*} \begin{aligned} C_{D} = C_{D0}(M) + \eta(M) C_{L\alpha}(M) \alpha^2 \\ C_{L} = C_{L\alpha}(M) \alpha \end{aligned} \end{equation*}
The coefficients \(C_{D0}(M)\), \(C_{L\alpha}(M)\) and \(\eta(M)\) are computed from the data:
\(M\) |
0 |
0.4 |
0.8 |
0.9 |
1.0 |
1.2 |
1.4 |
1.6 |
1.8 |
|---|---|---|---|---|---|---|---|---|---|
\(C_{L\alpha}\) |
3.44 |
3.44 |
3.44 |
3.58 |
4.44 |
3.44 |
3.01 |
2.86 |
2.44 |
\(C_{D0}\) |
0.013 |
0.013 |
0.013 |
0.014 |
0.031 |
0.041 |
0.039 |
0.036 |
0.035 |
\(\eta\) |
0.54 |
0.54 |
0.54 |
0.75 |
0.79 |
0.78 |
0.89 |
0.93 |
0.93 |
The trust \(T\) is computed as a function of mach number and altitude \(h\). Here the table lists thousands of lb thrust vs thousands of ft altitude.
\(M\) vs \(h\) |
0 |
5 |
10 |
15 |
20 |
25 |
30 |
40 |
50 |
70 |
|---|---|---|---|---|---|---|---|---|---|---|
0 |
24.2 |
|||||||||
0.2 |
28.0 |
24.6 |
21.1 |
18.1 |
15.2 |
12.8 |
10.7 |
|||
0.4 |
28.3 |
25.2 |
21.9 |
18.7 |
15.9 |
13.4 |
11.2 |
|||
0.6 |
30.8 |
27.2 |
23.8 |
20.5 |
17.3 |
14.7 |
12.3 |
8.1 |
4.9 |
|
0.8 |
34.5 |
30.3 |
26.6 |
23.2 |
19.8 |
16.8 |
14.1 |
9.4 |
5.6 |
1.1 |
1.0 |
37.9 |
34.3 |
30.4 |
26.8 |
23.3 |
19.8 |
16.8 |
11.2 |
6.8 |
1.4 |
1.2 |
36.1 |
38.0 |
34.9 |
31.3 |
27.3 |
23.6 |
20.1 |
13.4 |
8.3 |
1.7 |
1.4 |
36.6 |
38.5 |
36.1 |
31.6 |
28.1 |
24.2 |
16.2 |
10.0 |
2.2 |
|
1.6 |
38.7 |
35.7 |
32.0 |
28.1 |
19.3 |
11.9 |
2.9 |
|||
1.8 |
34.6 |
31.1 |
21.7 |
13.3 |
3.1 |
Summary of important constants
Data |
Value |
|---|---|
\(g\) |
\(32\, \text{ft/s}\) |
\(S_{\text{ref}}\) |
\(530\, \text{ft}^2\) |
\(I_{sp}\) |
\(1600\, \text{s}\) |
To close the system of equations, you must use data from the 1976 Standard Atmosphere to interpolate for \(\rho(h)\), the atmospheric density, \(a(h)\), the speed of sound and compute the Mach number as \(M = v / a(h)\). For the atmospheric data, you are welcome to use either of the following resources, as long as you cite in your report what you chose to utilize.
Digital Dutch Atmospheric Calculator: https://www.digitaldutch.com/atmoscalc/
OpenMDAO 1976 Atmosphere Component: https://github.com/OpenMDAO/dymos/blob/master/dymos/models/atmosphere/atmos_1976.py
With the dynamics fully specified, the initial conditions are:
\begin{equation*} \begin{aligned} h &= 0 \\ r &= 0 \\ v &= 380\, \text{ft/s} \\ \gamma &= 1.7^{\circ} \\ m &= 1304\, \text{slug} \\ \end{aligned} \end{equation*}
The final conditions are given by:
\begin{equation*} \begin{aligned} h &= 65617 \, \text{ft} \\ v &= 986.5 \, \text{ft/s} \\ \gamma &= 0 \\ \end{aligned} \end{equation*}
The objective is to minimize the time to reach the final conditions.
Assignment¶
Implement the equations of motion in descriptor form. (You may use the course notes code as a starting point, but you may not use Dymos.) Use the trapezoid rule for implicit time integration described in the course notes.
Develop surrogate models for the thrust \(T(M, h)\), the coefficients \(C_{D0}(M)\), \(C_{L\alpha}(M)\) and \(\eta(M)\) and the atmospheric density \(\rho(h)\) and speed of sound \(a(h)\). Validate these models and make sure that the units of the inputs and outputs are consistent.
Implement a reduced space method for the problem, treating the angle of attack \(\alpha(t)\) and the final time \(t_{f}\), as the design variables. The remaining variables \(h\), \(r\), \(v\), \(\gamma\) and \(m\) are the state variables. The constraints consist of the discrepancy between the desired and achieved final conditions.
Implement an adjoint method for the final condition constraints. Verify your implementation of the adjoint method using the complex-step method.
Implement a full space method for the problem, treating all variables as design variables controlled by the optimizer. The constraints consist of the initial conditions, the governing equations at each point in time and the final conditions.
Verify your implementation of the constraint Jacobian using complex step. Make sure that all components of the Jacobian matrix are correct!
Solve the optimization problem with the reduced space formulation.
Solve the optimization problem with the full space formulation.
Report requirements¶
Write a five page report with a 300 word abstract. You may use the AIAA paper format if you wish. Keep your descriptions concise and make your plots clear with labels and concise captions.
Write out a description of your optimization problem for the full space and reduced space formulations.
Describe your approach for developing the required surrogate models. Plot each of the 6 surrogate models you used.
Describe the reduced space formulation and the adjoint method you implemented. Plot the accuracy of the adjoint method vs step size. What accuracy did you achieve compared to the complex-step method?
Plot the accuracy of your full space Jacobian vs step size. What accuracy does your full space Jacobian achieve?
Describe issues you found with your optimization results. What was successful. What modifications did you make? How did you overcome any challenges.
Describe key observations you found from the optimization results. How does your problem change as the inputs change?