Multiphase Vehicle Launch#
Mathematical Formulation#
Problem Statement#
Find the optimal thrust direction vectors \(\mathbf{u}_i(t) = [u_{i,1}, u_{i,2}, u_{i,3}]^T\) for each phase that maximize the final payload mass:
(1)#\[J = -m_4(t_f)\]
Subject to the atmospheric flight dynamics with Earth rotation:
\[\frac{d\mathbf{r}}{dt} = \mathbf{v}\]
\[\frac{d\mathbf{v}}{dt} = \mathbf{a}_{\text{thrust}} + \mathbf{a}_{\text{drag}} + \mathbf{a}_{\text{grav}}\]
\[\frac{dm}{dt} = -\dot{m}_{\text{prop}}\]
Where:
\(\mathbf{a}_{\text{thrust}} = \frac{T_{\text{total}}}{m} \mathbf{u}\)
\(\mathbf{a}_{\text{drag}} = -\frac{\rho |\mathbf{v}_{\text{rel}}|}{2m} S C_D \mathbf{v}_{\text{rel}}\)
\(\mathbf{a}_{\text{grav}} = -\frac{\mu}{|\mathbf{r}|^3} \mathbf{r}\)
\(\mathbf{v}_{\text{rel}} = \mathbf{v} - \boldsymbol{\Omega} \times \mathbf{r}\)
Phase Configuration#
Phase 1 (\(t \in [0, 75.2]\) s): 6 SRBs + First Stage
Thrust: \(T_1 = 6 \times 628500 + 1083100 = 4854100\) N
Mass flow: \(\dot{m}_1 = 6 \times 17010/75.2 + 95550/261 = 1719.7\) kg/s
Phase 2 (\(t \in [75.2, 150.4]\) s): 3 SRBs + First Stage
Thrust: \(T_2 = 3 \times 628500 + 1083100 = 2968600\) N
Mass flow: \(\dot{m}_2 = 3 \times 17010/75.2 + 95550/261 = 1043.7\) kg/s
Phase 3 (\(t \in [150.4, 261.0]\) s): First Stage Only
Thrust: \(T_3 = 1083100\) N
Mass flow: \(\dot{m}_3 = 95550/261 = 366.1\) kg/s
Phase 4 (\(t \in [261.0, 961.0]\) s): Second Stage Only
Thrust: \(T_4 = 110094\) N
Mass flow: \(\dot{m}_4 = 16820/700 = 24.0\) kg/s
Terminal Orbital Constraints#
The final state must satisfy the target orbital elements:
Semi-major axis: \(a = 24361140\) m
Eccentricity: \(e = 0.7308\)
Inclination: \(i = 28.5°\)
Right ascension of ascending node: \(\Omega = 269.8°\)
Argument of periapsis: \(\omega = 130.5°\)
Boundary Conditions#
Initial Conditions:
Position: \(\mathbf{r}(0) = [5505524, 0, 3050952]^T\) m
Velocity: \(\mathbf{v}(0) = [0, 403.2, 0]^T\) m/s (Earth rotation)
Mass: \(m(0) = 631464\) kg (total vehicle)
Final Conditions:
Orbital elements as specified above
Mass: \(m(t_f) = \text{maximized}\)
Control Constraints:
\(|\mathbf{u}_i| = 1\) for all phases (unit thrust vector)
Physical Parameters#
Earth gravitational parameter: \(\mu = 3.986012 \times 10^{14}\) m³/s²
Earth radius: \(R_E = 6378145\) m
Earth rotation rate: \(\Omega = 7.29211585 \times 10^{-5}\) rad/s
Atmospheric density: \(\rho = 1.225 e^{-h/7200}\) kg/m³
Vehicle parameters: \(S = 4\pi\) m², \(C_D = 0.5\)
Vehicle Mass Breakdown#
Solid Rocket Boosters: 9 × 19290 kg total, 9 × 17010 kg propellant
First Stage: 104380 kg total, 95550 kg propellant
Second Stage: 19300 kg total, 16820 kg propellant
Payload: 4164 kg
State Variables#
\(\mathbf{r}(t) = [x, y, z]^T\): Position vector (m)
\(\mathbf{v}(t) = [v_x, v_y, v_z]^T\): Velocity vector (m/s)
\(m(t)\): Vehicle mass (kg)
Control Variables#
\(\mathbf{u}_i(t) = [u_{i,1}, u_{i,2}, u_{i,3}]^T\): Unit thrust direction vector for phase \(i\)
Notes#
This problem demonstrates multiphase rocket trajectory optimization with realistic propulsion systems, atmospheric effects, Earth rotation, and orbital mechanics constraints. Each phase has different thrust and mass flow characteristics corresponding to stage separations.
Running This Example#
cd examples/multiphase_vehicle_launch
python multiphase_vehicle_launch.py
Code Implementation#
examples/multiphase_vehicle_launch/multiphase_vehicle_launch.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# Earth and gravitational parameters
12OMEGA_EARTH = 7.29211585e-5 # Earth rotation rate (rad/s)
13MU_EARTH = 3.986012e14 # Gravitational parameter (m^3/s^2)
14RE_EARTH = 6378145.0 # Earth radius (m)
15G0 = 9.80665 # Sea level gravity (m/s^2)
16
17# Atmospheric parameters
18RHO0 = 1.225 # Sea level density (kg/m^3)
19H_SCALE_HEIGHT = 7200.0 # Density scale height (m)
20CD = 0.5 # Drag coefficient
21SA = 4 * np.pi # Surface area (m^2)
22
23# Vehicle mass parameters (kg)
24M_TOT_SRB = 19290.0
25M_PROP_SRB = 17010.0
26M_DRY_SRB = M_TOT_SRB - M_PROP_SRB
27M_TOT_FIRST = 104380.0
28M_PROP_FIRST = 95550.0
29M_DRY_FIRST = M_TOT_FIRST - M_PROP_FIRST
30M_TOT_SECOND = 19300.0
31M_PROP_SECOND = 16820.0
32M_DRY_SECOND = M_TOT_SECOND - M_PROP_SECOND
33M_PAYLOAD = 4164.0
34
35# Thrust parameters (N)
36THRUST_SRB = 628500.0
37THRUST_FIRST = 1083100.0
38THRUST_SECOND = 110094.0
39
40# Mission timeline (s)
41T_PHASE_1_END = 75.2
42T_PHASE_2_END = 150.4
43T_PHASE_3_END = 261.0
44T_PHASE_4_END = 961.0
45
46# Derived propulsion parameters
47MDOT_SRB = M_PROP_SRB / T_PHASE_1_END
48MDOT_FIRST = M_PROP_FIRST / T_PHASE_3_END
49MDOT_SECOND = M_PROP_SECOND / 700.0
50ISP_SRB = THRUST_SRB / (G0 * MDOT_SRB)
51ISP_FIRST = THRUST_FIRST / (G0 * MDOT_FIRST)
52ISP_SECOND = THRUST_SECOND / (G0 * MDOT_SECOND)
53
54# Numerical scaling factors
55R_SCALE = 1e6 # Position scaling (m)
56V_SCALE = 1e3 # Velocity scaling (m/s)
57M_SCALE = 1e4 # Mass scaling (kg)
58
59# Earth rotation matrix
60OMEGA_MATRIX = ca.MX(3, 3)
61OMEGA_MATRIX[0, 0] = 0.0
62OMEGA_MATRIX[0, 1] = -OMEGA_EARTH
63OMEGA_MATRIX[0, 2] = 0.0
64OMEGA_MATRIX[1, 0] = OMEGA_EARTH
65OMEGA_MATRIX[1, 1] = 0.0
66OMEGA_MATRIX[1, 2] = 0.0
67OMEGA_MATRIX[2, 0] = 0.0
68OMEGA_MATRIX[2, 1] = 0.0
69OMEGA_MATRIX[2, 2] = 0.0
70
71
72# ============================================================================
73# Initial and Target Conditions
74# ============================================================================
75
76LAT_LAUNCH = 28.5 * np.pi / 180.0
77X_LAUNCH = RE_EARTH * np.cos(LAT_LAUNCH)
78Y_LAUNCH = 0.0
79Z_LAUNCH = RE_EARTH * np.sin(LAT_LAUNCH)
80R_INITIAL_SCALED = [X_LAUNCH / R_SCALE, Y_LAUNCH / R_SCALE, Z_LAUNCH / R_SCALE]
81
82# Initial velocity due to Earth rotation
83omega_matrix_np = np.array([[0.0, -OMEGA_EARTH, 0.0], [OMEGA_EARTH, 0.0, 0.0], [0.0, 0.0, 0.0]])
84v_initial_phys = omega_matrix_np @ np.array([X_LAUNCH, Y_LAUNCH, Z_LAUNCH])
85V_INITIAL_SCALED = [
86 v_initial_phys[0] / V_SCALE,
87 v_initial_phys[1] / V_SCALE,
88 v_initial_phys[2] / V_SCALE,
89]
90
91# Target orbital elements
92A_TARGET = 24361140.0 # Semi-major axis (m)
93E_TARGET = 0.7308 # Eccentricity
94INC_TARGET = 28.5 * np.pi / 180.0 # Inclination (rad)
95RAAN_TARGET = 269.8 * np.pi / 180.0 # Right ascension of ascending node (rad)
96AOP_TARGET = 130.5 * np.pi / 180.0 # Argument of periapsis (rad)
97
98# Initial masses for each phase (scaled)
99M_INITIAL_P1 = (M_PAYLOAD + M_TOT_SECOND + M_TOT_FIRST + 9 * M_TOT_SRB) / M_SCALE
100M_FINAL_P1 = (M_INITIAL_P1 * M_SCALE - (6 * MDOT_SRB + MDOT_FIRST) * T_PHASE_1_END) / M_SCALE
101M_INITIAL_P2 = (M_FINAL_P1 * M_SCALE - 6 * M_DRY_SRB) / M_SCALE
102M_FINAL_P2 = (
103 M_INITIAL_P2 * M_SCALE - (3 * MDOT_SRB + MDOT_FIRST) * (T_PHASE_2_END - T_PHASE_1_END)
104) / M_SCALE
105M_INITIAL_P3 = (M_FINAL_P2 * M_SCALE - 3 * M_DRY_SRB) / M_SCALE
106M_FINAL_P3 = (M_INITIAL_P3 * M_SCALE - MDOT_FIRST * (T_PHASE_3_END - T_PHASE_2_END)) / M_SCALE
107M_INITIAL_P4 = (M_FINAL_P3 * M_SCALE - M_DRY_FIRST) / M_SCALE
108M_FINAL_P4 = M_PAYLOAD / M_SCALE
109
110# State bounds (scaled)
111R_MIN = -2 * RE_EARTH / R_SCALE
112R_MAX = 2 * RE_EARTH / R_SCALE
113V_MIN = -10000.0 / V_SCALE
114V_MAX = 10000.0 / V_SCALE
115
116
117# ============================================================================
118# Helper Functions
119# ============================================================================
120
121
122def _oe2rv(oe, mu):
123 # Convert orbital elements to position and velocity vectors.
124 a, e, i, Om, om, nu = oe
125 p = a * (1 - e * e)
126 r = p / (1 + e * np.cos(nu))
127
128 # Position in perifocal frame
129 rv_pf = np.array([r * np.cos(nu), r * np.sin(nu), 0.0])
130
131 # Velocity in perifocal frame
132 vv_pf = np.array([-np.sin(nu), e + np.cos(nu), 0.0]) * np.sqrt(mu / p)
133
134 # Rotation matrix from perifocal to inertial frame
135 cO, sO = np.cos(Om), np.sin(Om)
136 co, so = np.cos(om), np.sin(om)
137 ci, si = np.cos(i), np.sin(i)
138
139 R = np.array(
140 [
141 [cO * co - sO * so * ci, -cO * so - sO * co * ci, sO * si],
142 [sO * co + cO * so * ci, -sO * so + cO * co * ci, -cO * si],
143 [so * si, co * si, ci],
144 ]
145 )
146
147 ri = R @ rv_pf
148 vi = R @ vv_pf
149
150 return ri, vi
151
152
153def _cross_product(a, b):
154 # Cross product of two 3D vectors
155 return ca.vertcat(
156 a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]
157 )
158
159
160def _dot_product(a, b):
161 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
162
163
164def _smooth_heaviside(x, a_eps=0.1):
165 return 0.5 * (1 + ca.tanh(x / a_eps))
166
167
168def _rv2oe(rv, vv, mu):
169 # Convert position and velocity to orbital elements with smooth transitions
170 eps = 1e-12
171
172 K = ca.vertcat(0.0, 0.0, 1.0)
173 hv = _cross_product(rv, vv)
174 nv = _cross_product(K, hv)
175
176 n = ca.sqrt(ca.fmax(_dot_product(nv, nv), eps))
177 h2 = ca.fmax(_dot_product(hv, hv), eps)
178 v2 = ca.fmax(_dot_product(vv, vv), eps)
179 r = ca.sqrt(ca.fmax(_dot_product(rv, rv), eps))
180
181 rv_dot_vv = _dot_product(rv, vv)
182
183 ev = ca.vertcat(
184 (1 / mu) * ((v2 - mu / r) * rv[0] - rv_dot_vv * vv[0]),
185 (1 / mu) * ((v2 - mu / r) * rv[1] - rv_dot_vv * vv[1]),
186 (1 / mu) * ((v2 - mu / r) * rv[2] - rv_dot_vv * vv[2]),
187 )
188
189 p = h2 / mu
190 e = ca.sqrt(ca.fmax(_dot_product(ev, ev), eps))
191 a = p / (1 - e * e)
192 i = ca.acos(ca.fmax(ca.fmin(hv[2] / ca.sqrt(h2), 1.0 - eps), -1.0 + eps))
193
194 # Smooth transitions for angular elements
195 a_eps = 0.1
196 nv_dot_ev = _dot_product(nv, ev)
197
198 Om = _smooth_heaviside(nv[1] + eps, a_eps) * ca.acos(
199 ca.fmax(ca.fmin(nv[0] / n, 1.0 - eps), -1.0 + eps)
200 ) + _smooth_heaviside(-(nv[1] + eps), a_eps) * (
201 2 * np.pi - ca.acos(ca.fmax(ca.fmin(nv[0] / n, 1.0 - eps), -1.0 + eps))
202 )
203
204 om = _smooth_heaviside(ev[2], a_eps) * ca.acos(
205 ca.fmax(ca.fmin(nv_dot_ev / (n * e), 1.0 - eps), -1.0 + eps)
206 ) + _smooth_heaviside(-ev[2], a_eps) * (
207 2 * np.pi - ca.acos(ca.fmax(ca.fmin(nv_dot_ev / (n * e), 1.0 - eps), -1.0 + eps))
208 )
209
210 ev_dot_rv = _dot_product(ev, rv)
211 nu = _smooth_heaviside(rv_dot_vv, a_eps) * ca.acos(
212 ca.fmax(ca.fmin(ev_dot_rv / (e * r), 1.0 - eps), -1.0 + eps)
213 ) + _smooth_heaviside(-rv_dot_vv, a_eps) * (
214 2 * np.pi - ca.acos(ca.fmax(ca.fmin(ev_dot_rv / (e * r), 1.0 - eps), -1.0 + eps))
215 )
216
217 return ca.vertcat(a, e, i, Om, om, nu)
218
219
220def _define_phase_dynamics(phase, r_vars_s, v_vars_s, m_var_s, u_vars, phase_num):
221 # Convert scaled variables to physical units
222 r_vec_s = ca.vertcat(*r_vars_s)
223 v_vec_s = ca.vertcat(*v_vars_s)
224 u_vec = ca.vertcat(*u_vars)
225
226 r_vec = r_vec_s * R_SCALE
227 v_vec = v_vec_s * V_SCALE
228 m_phys = m_var_s * M_SCALE
229
230 rad = ca.sqrt(ca.fmax(_dot_product(r_vec, r_vec), 1e-12))
231
232 # Relative velocity accounting for Earth rotation
233 vrel = v_vec - OMEGA_MATRIX @ r_vec
234 speedrel = ca.sqrt(ca.fmax(_dot_product(vrel, vrel), 1e-12))
235
236 # Atmospheric drag
237 altitude = rad - RE_EARTH
238 rho = RHO0 * ca.exp(-altitude / H_SCALE_HEIGHT)
239 bc = (rho / (2 * m_phys)) * SA * CD
240 bcspeed = bc * speedrel
241 drag = -vrel * bcspeed
242
243 # Gravitational acceleration
244 muoverradcubed = MU_EARTH / (rad**3)
245 grav = -muoverradcubed * r_vec
246
247 # Phase-specific thrust and mass flow
248 if phase_num == 1:
249 T_tot = 6 * THRUST_SRB + THRUST_FIRST
250 mdot = -(6 * THRUST_SRB / (G0 * ISP_SRB) + THRUST_FIRST / (G0 * ISP_FIRST))
251 elif phase_num == 2:
252 T_tot = 3 * THRUST_SRB + THRUST_FIRST
253 mdot = -(3 * THRUST_SRB / (G0 * ISP_SRB) + THRUST_FIRST / (G0 * ISP_FIRST))
254 elif phase_num == 3:
255 T_tot = THRUST_FIRST
256 mdot = -THRUST_FIRST / (G0 * ISP_FIRST)
257 elif phase_num == 4:
258 T_tot = THRUST_SECOND
259 mdot = -THRUST_SECOND / (G0 * ISP_SECOND)
260
261 # Thrust acceleration
262 Toverm = T_tot / m_phys
263 thrust = Toverm * u_vec
264
265 # Total acceleration
266 acceleration = thrust + drag + grav
267
268 # Scaled derivatives
269 r_dot_scaled = v_vec / R_SCALE
270 v_dot_scaled = acceleration / V_SCALE
271 m_dot_scaled = mdot / M_SCALE
272
273 # Define dynamics
274 dynamics_dict = {}
275 for i in range(3):
276 dynamics_dict[r_vars_s[i]] = r_dot_scaled[i]
277 dynamics_dict[v_vars_s[i]] = v_dot_scaled[i]
278 dynamics_dict[m_var_s] = m_dot_scaled
279
280 phase.dynamics(dynamics_dict)
281
282 # Unit thrust vector constraint
283 thrust_magnitude_squared = _dot_product(u_vec, u_vec)
284 phase.path_constraints(thrust_magnitude_squared == 1.0)
285
286
287def _generate_phase_guess(N_list, r_init_s, v_init_s, m_init_s, m_final_s):
288 states_guess = []
289 controls_guess = []
290
291 for N in N_list:
292 # State guess (7 states × N+1 points)
293 states = np.zeros((7, N + 1))
294 for i in range(3):
295 states[i, :] = r_init_s[i] # Constant position
296 states[i + 3, :] = v_init_s[i] # Constant velocity
297 states[6, :] = np.linspace(m_init_s, m_final_s, N + 1) # Linear mass variation
298 states_guess.append(states)
299
300 # Control guess (3 controls × N points)
301 controls = np.zeros((3, N))
302 controls[0, :] = 1.0 # [1, 0, 0] unit vector
303 controls_guess.append(controls)
304
305 return states_guess, controls_guess
306
307
308# ============================================================================
309# Problem Setup
310# ============================================================================
311
312problem = mtor.Problem("Multiphase Vehicle Launch")
313
314
315# ============================================================================
316# Phase 1: Six Solid Rocket Boosters + First Stage
317# ============================================================================
318
319phase1 = problem.set_phase(1)
320t_1 = phase1.time(initial=0.0, final=T_PHASE_1_END)
321r1_s = [
322 phase1.state(f"r{i}_scaled", initial=R_INITIAL_SCALED[i], boundary=(R_MIN, R_MAX))
323 for i in range(3)
324]
325v1_s = [
326 phase1.state(f"v{i}_scaled", initial=V_INITIAL_SCALED[i], boundary=(V_MIN, V_MAX))
327 for i in range(3)
328]
329m1_s = phase1.state("mass_scaled", initial=M_INITIAL_P1, boundary=(M_FINAL_P1, M_INITIAL_P1))
330u1 = [phase1.control(f"u{i}", boundary=(-1.0, 1.0)) for i in range(3)]
331
332_define_phase_dynamics(phase1, r1_s, v1_s, m1_s, u1, 1)
333
334
335# ============================================================================
336# Phase 2: Three Solid Rocket Boosters + First Stage
337# ============================================================================
338
339phase2 = problem.set_phase(2)
340t_2 = phase2.time(initial=T_PHASE_1_END, final=T_PHASE_2_END)
341r2_s = [
342 phase2.state(f"r{i}_scaled", initial=r1_s[i].final, boundary=(R_MIN, R_MAX)) for i in range(3)
343]
344v2_s = [
345 phase2.state(f"v{i}_scaled", initial=v1_s[i].final, boundary=(V_MIN, V_MAX)) for i in range(3)
346]
347m2_s = phase2.state(
348 "mass_scaled", initial=m1_s.final - 6 * M_DRY_SRB / M_SCALE, boundary=(M_FINAL_P2, M_INITIAL_P2)
349)
350u2 = [phase2.control(f"u{i}", boundary=(-1.0, 1.0)) for i in range(3)]
351
352_define_phase_dynamics(phase2, r2_s, v2_s, m2_s, u2, 2)
353
354
355# ============================================================================
356# Phase 3: First Stage Only
357# ============================================================================
358
359phase3 = problem.set_phase(3)
360t_3 = phase3.time(initial=T_PHASE_2_END, final=T_PHASE_3_END)
361r3_s = [
362 phase3.state(f"r{i}_scaled", initial=r2_s[i].final, boundary=(R_MIN, R_MAX)) for i in range(3)
363]
364v3_s = [
365 phase3.state(f"v{i}_scaled", initial=v2_s[i].final, boundary=(V_MIN, V_MAX)) for i in range(3)
366]
367m3_s = phase3.state(
368 "mass_scaled", initial=m2_s.final - 3 * M_DRY_SRB / M_SCALE, boundary=(M_FINAL_P3, M_INITIAL_P3)
369)
370u3 = [phase3.control(f"u{i}", boundary=(-1.0, 1.0)) for i in range(3)]
371
372_define_phase_dynamics(phase3, r3_s, v3_s, m3_s, u3, 3)
373
374
375# ============================================================================
376# Phase 4: Second Stage Only
377# ============================================================================
378
379phase4 = problem.set_phase(4)
380t_4 = phase4.time(initial=T_PHASE_3_END, final=(T_PHASE_3_END, T_PHASE_4_END))
381r4_s = [
382 phase4.state(f"r{i}_scaled", initial=r3_s[i].final, boundary=(R_MIN, R_MAX)) for i in range(3)
383]
384v4_s = [
385 phase4.state(f"v{i}_scaled", initial=v3_s[i].final, boundary=(V_MIN, V_MAX)) for i in range(3)
386]
387m4_s = phase4.state(
388 "mass_scaled", initial=m3_s.final - M_DRY_FIRST / M_SCALE, boundary=(M_FINAL_P4, M_INITIAL_P4)
389)
390u4 = [phase4.control(f"u{i}", boundary=(-1.0, 1.0)) for i in range(3)]
391
392_define_phase_dynamics(phase4, r4_s, v4_s, m4_s, u4, 4)
393
394
395# ============================================================================
396# Target Orbital Elements Constraint
397# ============================================================================
398
399# Convert target orbital elements to position/velocity
400target_oe = [A_TARGET, E_TARGET, INC_TARGET, RAAN_TARGET, AOP_TARGET, 0.0]
401rout_phys, vout_phys = _oe2rv(target_oe, MU_EARTH)
402
403# Convert scaled final state to physical units for constraint evaluation
404r_final_s = ca.vertcat(*[r4_s[i].final for i in range(3)])
405v_final_s = ca.vertcat(*[v4_s[i].final for i in range(3)])
406r_final = r_final_s * R_SCALE
407v_final = v_final_s * V_SCALE
408
409# Compute final orbital elements
410oe_final = _rv2oe(r_final, v_final, MU_EARTH)
411
412# Constrain five orbital elements (true anomaly is free)
413phase4.event_constraints(
414 oe_final[0] == A_TARGET, # Semi-major axis
415 oe_final[1] == E_TARGET, # Eccentricity
416 oe_final[2] == INC_TARGET, # Inclination
417 oe_final[3] == RAAN_TARGET, # Right ascension of ascending node
418 oe_final[4] == AOP_TARGET, # Argument of periapsis
419)
420
421
422# ============================================================================
423# Objective Function
424# ============================================================================
425
426# Maximize final mass (minimize negative mass)
427problem.minimize(-m4_s.final)
428
429
430# ============================================================================
431# Mesh Configuration and Initial Guess
432# ============================================================================
433
434# Mesh configuration
435phase1.mesh([4, 4], [-1.0, 0.0, 1.0])
436phase2.mesh([4, 4], [-1.0, 0.0, 1.0])
437phase3.mesh([4, 4], [-1.0, 0.0, 1.0])
438phase4.mesh([4, 4], [-1.0, 0.0, 1.0])
439
440# Generate initial guesses for each phase
441states_p1, controls_p1 = _generate_phase_guess(
442 [4, 4], R_INITIAL_SCALED, V_INITIAL_SCALED, M_INITIAL_P1, M_FINAL_P1
443)
444states_p2, controls_p2 = _generate_phase_guess(
445 [4, 4], R_INITIAL_SCALED, V_INITIAL_SCALED, M_INITIAL_P2, M_FINAL_P2
446)
447states_p3, controls_p3 = _generate_phase_guess(
448 [4, 4], R_INITIAL_SCALED, V_INITIAL_SCALED, M_INITIAL_P3, M_FINAL_P3
449)
450
451# Target position and velocity for Phase 4 guess
452rout_scaled = rout_phys / R_SCALE
453vout_scaled = vout_phys / V_SCALE
454states_p4, controls_p4 = _generate_phase_guess(
455 [4, 4], rout_scaled, vout_scaled, M_INITIAL_P4, M_FINAL_P4
456)
457
458phase1.guess(states=states_p1, controls=controls_p1)
459phase2.guess(states=states_p2, controls=controls_p2)
460phase3.guess(states=states_p3, controls=controls_p3)
461phase4.guess(states=states_p4, controls=controls_p4)
462
463
464# ============================================================================
465# Solve
466# ============================================================================
467
468solution = mtor.solve_adaptive(
469 problem,
470 error_tolerance=1e-6,
471 max_iterations=20,
472 min_polynomial_degree=3,
473 max_polynomial_degree=8,
474 nlp_options={
475 "ipopt.max_iter": 1000,
476 "ipopt.tol": 1e-6,
477 "ipopt.constr_viol_tol": 1e-6,
478 "ipopt.linear_solver": "mumps",
479 "ipopt.print_level": 5,
480 },
481)
482
483
484# ============================================================================
485# Results Analysis
486# ============================================================================
487
488if solution.status["success"]:
489 final_mass_scaled = -solution.status["objective"]
490 final_mass = final_mass_scaled * M_SCALE
491 mission_time = solution.status["total_mission_time"]
492
493 print("Multiphase Vehicle Launch Solution:")
494 print(f"Final mass: {final_mass:.1f} kg")
495 print(f"Mission time: {mission_time:.1f} seconds")
496 print(f"Payload fraction: {final_mass / (M_INITIAL_P1 * M_SCALE) * 100:.2f}%")
497
498 # Final state verification
499 r_final_scaled = (
500 solution[(4, "r0_scaled")][-1],
501 solution[(4, "r1_scaled")][-1],
502 solution[(4, "r2_scaled")][-1],
503 )
504 v_final_scaled = (
505 solution[(4, "v0_scaled")][-1],
506 solution[(4, "v1_scaled")][-1],
507 solution[(4, "v2_scaled")][-1],
508 )
509
510 print("\nFinal state verification:")
511 print(
512 f"Final position: [{r_final_scaled[0]:.6f}, {r_final_scaled[1]:.6f}, {r_final_scaled[2]:.6f}] (scaled)"
513 )
514 print(
515 f"Final velocity: [{v_final_scaled[0]:.6f}, {v_final_scaled[1]:.6f}, {v_final_scaled[2]:.6f}] (scaled)"
516 )
517
518 # Phase information
519 print("\nPhase Information:")
520 for phase_id, phase_data in solution.phases.items():
521 times = phase_data["times"]
522 print(f" Phase {phase_id}: {times['duration']:.1f}s duration")
523
524 solution.plot()
525else:
526 print(f"Failed: {solution.status['message']}")