Multi-Phase Cannonball#
Maximizing the range of a cannonball in a vacuum is a typical introductory problem for optimal control. In this example we are going to demonstrate a more multidisciplinary take on the problem. We will assume a density of the metal from which the cannonball is constructed, and a cannon that can fire any diameter cannonball but is limited to a maximum muzzle energy. If we make the cannonball large it will be heavy and the cannon will not be capable of propelling it very far. If we make the cannonball too small, it will have a low ballistic coefficient and not be able to sustain its momentum in the presence of atmospheric drag. Somewhere between these two extremes is the cannonball radius which allows for maximum range flight.
The presence of atmospheric drag also means that we typically want to launch the cannonball with more horizontal velocity, and thus use a launch angle less than 45 degrees.
The goal of our optimization is to find the optimal design for the cannonball (its radius) and the optimal flight profile (its launch angle) simultaneously.
Using two phases to capture an intermediate boundary constraint#
This problem demonstrates the use of two phases to capture the state of the system at an event in the trajectory. Here, we have the first phase (ascent) terminate when the flight path angle reaches zero (apogee). The descent phase follows until the cannonball impacts the ground.
The dynamics are given by
(70)#\[\begin{align}
\frac{dv}{dt} &= \frac{D}{m} - g \sin \gamma \\
\frac{d\gamma}{dt} &= - \frac{g \cos \gamma}{v} \\
\frac{dh}{dt} &= v \sin \gamma \\
\frac{dr}{dt} &= v \cos \gamma \\
\end{align}\]
The initial conditions are
(71)#\[\begin{align}
r_0 &= 0 \rm{\,m} \\
h_0 &= 100 \rm{\,m} \\
v_0 &= \rm{free} \\
\gamma_0 &= \rm{free}
\end{align}\]
and the final conditions are
(72)#\[\begin{align}
h_f &= 0 \rm{\,m}
\end{align}\]
Designing a cannonball for maximum range#
This problem demonstrates a very simple vehicle design capability that is run before the trajectory.
We assume our cannon can shoot a cannonball with some fixed kinetic energy and that our cannonball is made of solid iron. The volume (and mass) of the cannonball is proportional to its radius cubed, while the cross-sectional area is proportional to its radius squared. If we increase the size of the cannonball, the ballistic coefficient
(73)#\[\begin{align}
BC &= \frac{m}{C_D A}
\end{align}\]
will increase, meaning the cannonball overcome air resistance more easily and thus carry more distance.
However, making the cannonball larger also increases its mass. Our cannon can impart the cannonball with, at most, 400 kJ of kinetic energy. So making the cannonball larger will decrease the initial velocity, and thus negatively impact its range.
We therefore have a design that affects the objective in competing ways. We cannot make the cannonball too large, as it will be too heavy to shoot. We also cannot make the cannonball too small, as it will be more susceptible to air resistance. Somewhere in between is the sweet spot that provides the maximum range cannonball.
The cannonball sizing component#
This compoennt computes the area (needed to compute drag) and mass (needed to compute energy) of a cannonball of a given radius and density.
This component sits upstream of the trajectory model and feeds its outputs to the trajectory as parameters.
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.models.atmosphere.atmos_1976 import USatm1976Data
#############################################
# Component for the design part of the model
#############################################
class CannonballSizeComp(om.ExplicitComponent):
"""
Compute the area and mass of a cannonball with a given radius and density.
Notes
-----
This component is not vectorized with 'num_nodes' as is the usual way
with Dymos, but is instead intended to compute a scalar mass and reference
area from scalar radius and density inputs. This component does not reside
in the ODE but instead its outputs are connected to the trajectory via
input design parameters.
"""
def setup(self):
self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')
self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')
self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')
self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')
self.declare_partials(of='mass', wrt='dens')
self.declare_partials(of='mass', wrt='radius')
self.declare_partials(of='S', wrt='radius')
def compute(self, inputs, outputs):
radius = inputs['radius']
dens = inputs['dens']
outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3
outputs['S'] = np.pi * radius ** 2
def compute_partials(self, inputs, partials):
radius = inputs['radius']
dens = inputs['dens']
partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3
partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2
partials['S', 'radius'] = 2 * np.pi * radius
The cannonball ODE component#
This component computes the state rates and the kinetic energy of the cannonball.
By calling the declare_coloring method wrt all inputs and using method 'cs', we’re telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, and to automatically compute those outputs using complex-step approximation.
class CannonballODE(om.ExplicitComponent):
"""
Cannonball ODE assuming flat earth and accounting for air resistance
"""
def initialize(self):
self.options.declare('num_nodes', types=int)
def setup(self):
nn = self.options['num_nodes']
# static parameters
self.add_input('m', units='kg')
self.add_input('S', units='m**2')
# 0.5 good assumption for a sphere
self.add_input('CD', 0.5)
# time varying inputs
self.add_input('h', units='m', shape=nn)
self.add_input('v', units='m/s', shape=nn)
self.add_input('gam', units='rad', shape=nn)
# state rates
self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])
self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])
self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])
self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])
self.add_output('ke', shape=nn, units='J')
# Ask OpenMDAO to compute the partial derivatives using complex-step
# with a partial coloring algorithm for improved performance, and use
# a graph coloring algorithm to automatically detect the sparsity pattern.
self.declare_coloring(wrt='*', method='cs')
alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]
rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
self.rho_interp = interp1d(np.array(alt_data, dtype=complex),
np.array(rho_data, dtype=complex),
kind='linear')
def compute(self, inputs, outputs):
gam = inputs['gam']
v = inputs['v']
h = inputs['h']
m = inputs['m']
S = inputs['S']
CD = inputs['CD']
GRAVITY = 9.80665 # m/s**2
# handle complex-step gracefully from the interpolant
if np.iscomplexobj(h):
rho = self.rho_interp(inputs['h'])
else:
rho = self.rho_interp(inputs['h']).real
q = 0.5*rho*inputs['v']**2
qS = q * S
D = qS * CD
cgam = np.cos(gam)
sgam = np.sin(gam)
outputs['v_dot'] = - D/m-GRAVITY*sgam
outputs['gam_dot'] = -(GRAVITY/v)*cgam
outputs['h_dot'] = v*sgam
outputs['r_dot'] = v*cgam
outputs['ke'] = 0.5*m*v**2
Building and running the problem#
The following code defines the components for the physical cannonball calculations and ODE problem, sets up trajectory using two phases, and links them accordingly. The initial flight path angle is free, since 45 degrees is not necessarily optimal once air resistance is taken into account.
p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = 'SLSQP'
p.driver.declare_coloring()
p.model.add_subsystem('size_comp', CannonballSizeComp(),
promotes_inputs=['radius', 'dens'])
p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
p.model.add_design_var('radius', lower=0.01, upper=0.10,
ref0=0.01, ref=0.10, units='m')
traj = p.model.add_subsystem('traj', dm.Trajectory())
transcription = dm.Radau(num_segments=5, order=3, compressed=True)
ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
ascent = traj.add_phase('ascent', ascent)
# All initial states except flight path angle are fixed
# Final flight path angle is fixed (we will set it to zero
# so that the phase ends at apogee).
# The output of the ODE which provides the rate source for each state
# is obtained from the tags used on those outputs in the ODE.
# The units of the states are automatically inferred by multiplying the units
# of those rates by the time units.
ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
duration_ref=100, units='s')
ascent.set_state_options('r', fix_initial=True, fix_final=False)
ascent.set_state_options('h', fix_initial=True, fix_final=False)
ascent.set_state_options('gam', fix_initial=False, fix_final=True)
ascent.set_state_options('v', fix_initial=False, fix_final=False)
ascent.add_parameter('S', units='m**2', static_target=True)
ascent.add_parameter('m', units='kg', static_target=True)
# Limit the muzzle energy
ascent.add_boundary_constraint('ke', loc='initial',
upper=400000, lower=0, ref=100000)
# Second Phase (descent)
transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
traj.add_phase('descent', descent)
# All initial states and time are free, since
# they will be linked to the final states of ascent.
# Final altitude is fixed, because we will set
# it to zero so that the phase ends at ground impact)
descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
duration_ref=100, units='s')
descent.add_state('r')
descent.add_state('h', fix_initial=False, fix_final=True)
descent.add_state('gam', fix_initial=False, fix_final=False)
descent.add_state('v', fix_initial=False, fix_final=False)
descent.add_parameter('S', units='m**2', static_target=True)
descent.add_parameter('m', units='kg', static_target=True)
descent.add_objective('r', loc='final', scaler=-1.0)
# Add internally-managed design parameters to the trajectory.
traj.add_parameter('CD',
targets={'ascent': ['CD'], 'descent': ['CD']},
val=0.5, units=None, opt=False, static_target=True)
# Add externally-provided design parameters to the trajectory.
# In this case, we connect 'm' to pre-existing input parameters
# named 'mass' in each phase.
traj.add_parameter('m', units='kg', val=1.0,
targets={'ascent': 'm', 'descent': 'm'},
static_target=True)
# In this case, by omitting targets, we're connecting these
# parameters to parameters with the same name in each phase.
traj.add_parameter('S', units='m**2', val=0.005, static_target=True)
# Link Phases (link time and all state variables)
traj.link_phases(phases=['ascent', 'descent'], vars=['*'])
# Issue Connections
p.model.connect('size_comp.mass', 'traj.parameters:m')
p.model.connect('size_comp.S', 'traj.parameters:S')
# A linear solver at the top level can improve performance.
p.model.linear_solver = om.DirectSolver()
# Finish Problem Setup
p.setup()
#############################################
# Set constants and initial guesses
#############################################
p.set_val('radius', 0.05, units='m')
p.set_val('dens', 7.87, units='g/cm**3')
p.set_val('traj.parameters:CD', 0.5)
p.set_val('traj.ascent.t_initial', 0.0)
p.set_val('traj.ascent.t_duration', 10.0)
p.set_val('traj.ascent.states:r', ascent.interp('r', [0, 100]))
p.set_val('traj.ascent.states:h', ascent.interp('h', [0, 100]))
p.set_val('traj.ascent.states:v', ascent.interp('v', [200, 150]))
p.set_val('traj.ascent.states:gam', ascent.interp('gam', [25, 0]), units='deg')
p.set_val('traj.descent.t_initial', 10.0)
p.set_val('traj.descent.t_duration', 10.0)
p.set_val('traj.descent.states:r', descent.interp('r', [100, 200]))
p.set_val('traj.descent.states:h', descent.interp('h', [100, 0]))
p.set_val('traj.descent.states:v', descent.interp('v', [150, 200]))
p.set_val('traj.descent.states:gam', descent.interp('gam', [0, -45]), units='deg')
#####################################################
# Run the optimization and final explicit simulation
#####################################################
dm.run_problem(p, simulate=True)
--- Constraint Report [traj] ---
--- ascent ---
[initial] 0.0000e+00 <= ke <= 4.0000e+05 [J]
--- descent ---
None
Model viewer data has already been recorded for Driver.
Model viewer data has already been recorded for Driver.
Coloring for 'traj.phases.ascent.rhs_all' (class CannonballODE)
Jacobian shape: (100, 63) (4.44% nonzero)
FWD solves: 6 REV solves: 0
Total colors vs. total size: 6 vs 63 (90.48% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.0214 sec
Time to compute coloring: 0.0088 sec
Memory to compute coloring: 0.2500 MB
Coloring for 'traj.phases.descent.rhs_disc' (class CannonballODE)
Jacobian shape: (50, 33) (8.48% nonzero)
FWD solves: 6 REV solves: 0
Total colors vs. total size: 6 vs 33 (81.82% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.0113 sec
Time to compute coloring: 0.0061 sec
Memory to compute coloring: 0.1250 MB
Coloring for 'traj.phases.descent.rhs_col' (class CannonballODE)
Jacobian shape: (25, 18) (15.56% nonzero)
FWD solves: 6 REV solves: 0
Total colors vs. total size: 6 vs 18 (66.67% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.0065 sec
Time to compute coloring: 0.0027 sec
Memory to compute coloring: 0.0000 MB
Full total jacobian was computed 3 times, taking 0.115493 seconds.
Total jacobian shape: (87, 88)
Jacobian shape: (87, 88) (13.32% nonzero)
FWD solves: 27 REV solves: 0
Total colors vs. total size: 27 vs 88 (69.32% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.1155 sec
Time to compute coloring: 0.0610 sec
Memory to compute coloring: 0.3750 MB
Coloring created on: 2023-12-08 21:03:53
Optimization Problem -- Optimization using pyOpt_sparse
================================================================================
Objective Function: _objfunc
Solution:
--------------------------------------------------------------------------------
Total Time: 1.5119
User Objective Time : 0.2139
User Sensitivity Time : 0.9138
Interface Time : 0.2255
Opt Solver Time: 0.1587
Calls to Objective Function : 118
Calls to Sens Function : 89
Objectives
Index Name Value
0 traj.descent.states:r -3.183118E+03
Variables (c - continuous, i - integer, d - discrete)
Index Name Type Lower Bound Value Upper Bound Status
0 radius_0 c 0.000000E+00 3.539225E-01 1.000000E+00
1 traj.ascent.t_duration_0 c 1.000000E-02 1.064895E-01 1.000000E+00
2 traj.ascent.states:r_0 c -1.000000E+21 3.231043E+02 1.000000E+21
3 traj.ascent.states:r_1 c -1.000000E+21 6.616876E+02 1.000000E+21
4 traj.ascent.states:r_2 c -1.000000E+21 7.552526E+02 1.000000E+21
5 traj.ascent.states:r_3 c -1.000000E+21 9.434454E+02 1.000000E+21
6 traj.ascent.states:r_4 c -1.000000E+21 1.165633E+03 1.000000E+21
7 traj.ascent.states:r_5 c -1.000000E+21 1.229354E+03 1.000000E+21
8 traj.ascent.states:r_6 c -1.000000E+21 1.364715E+03 1.000000E+21
9 traj.ascent.states:r_7 c -1.000000E+21 1.532445E+03 1.000000E+21
10 traj.ascent.states:r_8 c -1.000000E+21 1.581810E+03 1.000000E+21
11 traj.ascent.states:r_9 c -1.000000E+21 1.688890E+03 1.000000E+21
12 traj.ascent.states:r_10 c -1.000000E+21 1.825079E+03 1.000000E+21
13 traj.ascent.states:r_11 c -1.000000E+21 1.865797E+03 1.000000E+21
14 traj.ascent.states:r_12 c -1.000000E+21 1.955144E+03 1.000000E+21
15 traj.ascent.states:r_13 c -1.000000E+21 2.070557E+03 1.000000E+21
16 traj.ascent.states:r_14 c -1.000000E+21 2.105398E+03 1.000000E+21
17 traj.ascent.states:h_0 c -1.000000E+21 1.995146E+02 1.000000E+21
18 traj.ascent.states:h_1 c -1.000000E+21 4.002917E+02 1.000000E+21
19 traj.ascent.states:h_2 c -1.000000E+21 4.537845E+02 1.000000E+21
20 traj.ascent.states:h_3 c -1.000000E+21 5.574056E+02 1.000000E+21
21 traj.ascent.states:h_4 c -1.000000E+21 6.711526E+02 1.000000E+21
22 traj.ascent.states:h_5 c -1.000000E+21 7.016712E+02 1.000000E+21
23 traj.ascent.states:h_6 c -1.000000E+21 7.625774E+02 1.000000E+21
24 traj.ascent.states:h_7 c -1.000000E+21 8.292929E+02 1.000000E+21
25 traj.ascent.states:h_8 c -1.000000E+21 8.467852E+02 1.000000E+21
26 traj.ascent.states:h_9 c -1.000000E+21 8.807923E+02 1.000000E+21
27 traj.ascent.states:h_10 c -1.000000E+21 9.151926E+02 1.000000E+21
28 traj.ascent.states:h_11 c -1.000000E+21 9.233143E+02 1.000000E+21
29 traj.ascent.states:h_12 c -1.000000E+21 9.371859E+02 1.000000E+21
30 traj.ascent.states:h_13 c -1.000000E+21 9.461925E+02 1.000000E+21
31 traj.ascent.states:h_14 c -1.000000E+21 9.467365E+02 1.000000E+21
32 traj.ascent.states:gam_0 c -1.000000E+21 5.588459E-01 1.000000E+21
33 traj.ascent.states:gam_1 c -1.000000E+21 5.461378E-01 1.000000E+21
34 traj.ascent.states:gam_2 c -1.000000E+21 5.228406E-01 1.000000E+21
35 traj.ascent.states:gam_3 c -1.000000E+21 5.140981E-01 1.000000E+21
36 traj.ascent.states:gam_4 c -1.000000E+21 4.911599E-01 1.000000E+21
37 traj.ascent.states:gam_5 c -1.000000E+21 4.530243E-01 1.000000E+21
38 traj.ascent.states:gam_6 c -1.000000E+21 4.393302E-01 1.000000E+21
39 traj.ascent.states:gam_7 c -1.000000E+21 4.048232E-01 1.000000E+21
40 traj.ascent.states:gam_8 c -1.000000E+21 3.498022E-01 1.000000E+21
41 traj.ascent.states:gam_9 c -1.000000E+21 3.305300E-01 1.000000E+21
42 traj.ascent.states:gam_10 c -1.000000E+21 2.829154E-01 1.000000E+21
43 traj.ascent.states:gam_11 c -1.000000E+21 2.091969E-01 1.000000E+21
44 traj.ascent.states:gam_12 c -1.000000E+21 1.839162E-01 1.000000E+21
45 traj.ascent.states:gam_13 c -1.000000E+21 1.226062E-01 1.000000E+21
46 traj.ascent.states:gam_14 c -1.000000E+21 3.070674E-02 1.000000E+21
47 traj.ascent.states:v_0 c -1.000000E+21 5.753380E+02 1.000000E+21
48 traj.ascent.states:v_1 c -1.000000E+21 4.377128E+02 1.000000E+21
49 traj.ascent.states:v_2 c -1.000000E+21 3.333818E+02 1.000000E+21
50 traj.ascent.states:v_3 c -1.000000E+21 3.071231E+02 1.000000E+21
51 traj.ascent.states:v_4 c -1.000000E+21 2.628446E+02 1.000000E+21
52 traj.ascent.states:v_5 c -1.000000E+21 2.189809E+02 1.000000E+21
53 traj.ascent.states:v_6 c -1.000000E+21 2.076302E+02 1.000000E+21
54 traj.ascent.states:v_7 c -1.000000E+21 1.856231E+02 1.000000E+21
55 traj.ascent.states:v_8 c -1.000000E+21 1.616748E+02 1.000000E+21
56 traj.ascent.states:v_9 c -1.000000E+21 1.552355E+02 1.000000E+21
57 traj.ascent.states:v_10 c -1.000000E+21 1.422744E+02 1.000000E+21
58 traj.ascent.states:v_11 c -1.000000E+21 1.276786E+02 1.000000E+21
59 traj.ascent.states:v_12 c -1.000000E+21 1.237104E+02 1.000000E+21
60 traj.ascent.states:v_13 c -1.000000E+21 1.156806E+02 1.000000E+21
61 traj.ascent.states:v_14 c -1.000000E+21 1.067082E+02 1.000000E+21
62 traj.ascent.states:v_15 c -1.000000E+21 1.043171E+02 1.000000E+21
63 traj.descent.t_initial_0 c 5.000000E-01 1.064895E+01 1.000000E+02
64 traj.descent.t_duration_0 c 5.000000E-03 1.617267E-01 1.000000E+00
65 traj.descent.states:r_0 c -1.000000E+21 2.105398E+03 1.000000E+21
66 traj.descent.states:r_1 c -1.000000E+21 2.410813E+03 1.000000E+21
67 traj.descent.states:r_2 c -1.000000E+21 2.663597E+03 1.000000E+21
68 traj.descent.states:r_3 c -1.000000E+21 2.873278E+03 1.000000E+21
69 traj.descent.states:r_4 c -1.000000E+21 3.045011E+03 1.000000E+21
70 traj.descent.states:r_5 c -1.000000E+21 3.183118E+03 1.000000E+21
71 traj.descent.states:h_0 c -1.000000E+21 9.467365E+02 1.000000E+21
72 traj.descent.states:h_1 c -1.000000E+21 8.986347E+02 1.000000E+21
73 traj.descent.states:h_2 c -1.000000E+21 7.651714E+02 1.000000E+21
74 traj.descent.states:h_3 c -1.000000E+21 5.609894E+02 1.000000E+21
75 traj.descent.states:h_4 c -1.000000E+21 3.009916E+02 1.000000E+21
76 traj.descent.states:gam_0 c -1.000000E+21 0.000000E+00 1.000000E+21
77 traj.descent.states:gam_1 c -1.000000E+21 -3.251038E-01 1.000000E+21
78 traj.descent.states:gam_2 c -1.000000E+21 -6.397506E-01 1.000000E+21
79 traj.descent.states:gam_3 c -1.000000E+21 -8.901606E-01 1.000000E+21
80 traj.descent.states:gam_4 c -1.000000E+21 -1.071739E+00 1.000000E+21
81 traj.descent.states:gam_5 c -1.000000E+21 -1.201117E+00 1.000000E+21
82 traj.descent.states:v_0 c -1.000000E+21 1.043171E+02 1.000000E+21
83 traj.descent.states:v_1 c -1.000000E+21 9.031921E+01 1.000000E+21
84 traj.descent.states:v_2 c -1.000000E+21 8.870119E+01 1.000000E+21
85 traj.descent.states:v_3 c -1.000000E+21 9.333723E+01 1.000000E+21
86 traj.descent.states:v_4 c -1.000000E+21 9.960765E+01 1.000000E+21
87 traj.descent.states:v_5 c -1.000000E+21 1.050638E+02 1.000000E+21
Constraints (i - inequality, e - equality)
Index Name Type Lower Value Upper Status Lagrange Multiplier (N/A)
0 traj.linkages.ascent:time_final|descent:time_initial e 0.000000E+00 -1.776357E-15 0.000000E+00 9.00000E+100
1 traj.linkages.ascent:r_final|descent:r_initial e 0.000000E+00 0.000000E+00 0.000000E+00 9.00000E+100
2 traj.linkages.ascent:h_final|descent:h_initial e 0.000000E+00 0.000000E+00 0.000000E+00 9.00000E+100
3 traj.linkages.ascent:gam_final|descent:gam_initial e 0.000000E+00 0.000000E+00 0.000000E+00 9.00000E+100
4 traj.linkages.ascent:v_final|descent:v_initial e 0.000000E+00 0.000000E+00 0.000000E+00 9.00000E+100
5 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 1.343816E-11 0.000000E+00 9.00000E+100
6 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 7.203338E-12 0.000000E+00 9.00000E+100
7 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 3.450338E-12 0.000000E+00 9.00000E+100
8 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 2.723951E-12 0.000000E+00 9.00000E+100
9 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 1.452774E-12 0.000000E+00 9.00000E+100
10 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 4.842580E-13 0.000000E+00 9.00000E+100
11 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -1.815968E-13 0.000000E+00 9.00000E+100
12 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -1.513306E-13 0.000000E+00 9.00000E+100
13 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -1.815968E-13 0.000000E+00 9.00000E+100
14 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 8.777176E-13 0.000000E+00 9.00000E+100
15 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -8.474515E-13 0.000000E+00 9.00000E+100
16 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -8.171854E-13 0.000000E+00 9.00000E+100
17 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -9.685160E-13 0.000000E+00 9.00000E+100
18 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 -8.474515E-13 0.000000E+00 9.00000E+100
19 traj.ascent.collocation_constraint.defects:r e 0.000000E+00 1.361976E-13 0.000000E+00 9.00000E+100
20 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 9.382499E-12 0.000000E+00 9.00000E+100
21 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 3.299008E-12 0.000000E+00 9.00000E+100
22 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 1.059314E-12 0.000000E+00 9.00000E+100
23 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 5.145241E-13 0.000000E+00 9.00000E+100
24 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 9.079838E-14 0.000000E+00 9.00000E+100
25 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 2.421290E-13 0.000000E+00 9.00000E+100
26 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 1.513306E-13 0.000000E+00 9.00000E+100
27 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 1.967298E-13 0.000000E+00 9.00000E+100
28 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 2.648286E-13 0.000000E+00 9.00000E+100
29 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 9.382499E-13 0.000000E+00 9.00000E+100
30 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 4.161592E-13 0.000000E+00 9.00000E+100
31 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 3.594102E-13 0.000000E+00 9.00000E+100
32 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 3.404939E-13 0.000000E+00 9.00000E+100
33 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 1.740302E-13 0.000000E+00 9.00000E+100
34 traj.ascent.collocation_constraint.defects:h e 0.000000E+00 5.755293E-13 0.000000E+00 9.00000E+100
35 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 7.725399E-15 0.000000E+00 9.00000E+100
36 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 4.067750E-15 0.000000E+00 9.00000E+100
37 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 2.357152E-15 0.000000E+00 9.00000E+100
38 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 1.378084E-15 0.000000E+00 9.00000E+100
39 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 7.019731E-16 0.000000E+00 9.00000E+100
40 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 4.433514E-17 0.000000E+00 9.00000E+100
41 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 1.995082E-16 0.000000E+00 9.00000E+100
42 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 -1.108379E-16 0.000000E+00 9.00000E+100
43 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 8.128110E-17 0.000000E+00 9.00000E+100
44 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 2.216757E-16 0.000000E+00 9.00000E+100
45 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 3.399028E-16 0.000000E+00 9.00000E+100
46 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 7.832542E-16 0.000000E+00 9.00000E+100
47 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 9.014813E-16 0.000000E+00 9.00000E+100
48 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 1.433503E-15 0.000000E+00 9.00000E+100
49 traj.ascent.collocation_constraint.defects:gam e 0.000000E+00 2.142865E-15 0.000000E+00 9.00000E+100
50 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 3.021770E-10 0.000000E+00 9.00000E+100
51 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 9.890970E-11 0.000000E+00 9.00000E+100
52 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 3.700034E-11 0.000000E+00 9.00000E+100
53 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 2.549921E-11 0.000000E+00 9.00000E+100
54 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 1.425535E-11 0.000000E+00 9.00000E+100
55 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 7.445467E-12 0.000000E+00 9.00000E+100
56 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 6.015392E-12 0.000000E+00 9.00000E+100
57 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 4.040528E-12 0.000000E+00 9.00000E+100
58 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 2.489389E-12 0.000000E+00 9.00000E+100
59 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 2.220777E-12 0.000000E+00 9.00000E+100
60 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 1.587080E-12 0.000000E+00 9.00000E+100
61 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 1.063098E-12 0.000000E+00 9.00000E+100
62 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 1.049856E-12 0.000000E+00 9.00000E+100
63 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 7.661113E-13 0.000000E+00 9.00000E+100
64 traj.ascent.collocation_constraint.defects:v e 0.000000E+00 6.034309E-13 0.000000E+00 9.00000E+100
65 traj.descent.collocation_constraint.defects:r e 0.000000E+00 8.572566E-12 0.000000E+00 9.00000E+100
66 traj.descent.collocation_constraint.defects:r e 0.000000E+00 8.894325E-12 0.000000E+00 9.00000E+100
67 traj.descent.collocation_constraint.defects:r e 0.000000E+00 8.480635E-12 0.000000E+00 9.00000E+100
68 traj.descent.collocation_constraint.defects:r e 0.000000E+00 8.124403E-12 0.000000E+00 9.00000E+100
69 traj.descent.collocation_constraint.defects:r e 0.000000E+00 6.963774E-12 0.000000E+00 9.00000E+100
70 traj.descent.collocation_constraint.defects:h e 0.000000E+00 -1.985135E-12 0.000000E+00 9.00000E+100
71 traj.descent.collocation_constraint.defects:h e 0.000000E+00 -5.377964E-12 0.000000E+00 9.00000E+100
72 traj.descent.collocation_constraint.defects:h e 0.000000E+00 -6.733946E-12 0.000000E+00 9.00000E+100
73 traj.descent.collocation_constraint.defects:h e 0.000000E+00 -7.469394E-12 0.000000E+00 9.00000E+100
74 traj.descent.collocation_constraint.defects:h e 0.000000E+00 -8.296773E-12 0.000000E+00 9.00000E+100
75 traj.descent.collocation_constraint.defects:gam e 0.000000E+00 1.714729E-14 0.000000E+00 9.00000E+100
76 traj.descent.collocation_constraint.defects:gam e 0.000000E+00 5.902797E-15 0.000000E+00 9.00000E+100
77 traj.descent.collocation_constraint.defects:gam e 0.000000E+00 -8.910305E-15 0.000000E+00 9.00000E+100
78 traj.descent.collocation_constraint.defects:gam e 0.000000E+00 -1.264725E-14 0.000000E+00 9.00000E+100
79 traj.descent.collocation_constraint.defects:gam e 0.000000E+00 -1.163726E-14 0.000000E+00 9.00000E+100
80 traj.descent.collocation_constraint.defects:v e 0.000000E+00 7.526851E-13 0.000000E+00 9.00000E+100
81 traj.descent.collocation_constraint.defects:v e 0.000000E+00 4.379292E-13 0.000000E+00 9.00000E+100
82 traj.descent.collocation_constraint.defects:v e 0.000000E+00 1.276261E-12 0.000000E+00 9.00000E+100
83 traj.descent.collocation_constraint.defects:v e 0.000000E+00 2.268469E-12 0.000000E+00 9.00000E+100
84 traj.descent.collocation_constraint.defects:v e 0.000000E+00 2.804973E-12 0.000000E+00 9.00000E+100
85 traj.phases.ascent->initial_boundary_constraint->ke i 0.000000E+00 4.000000E+00 4.000000E+00 u 9.00000E+100
Exit Status
Inform Description
0 Optimization terminated successfully.
--------------------------------------------------------------------------------
Simulating trajectory traj
Model viewer data has already been recorded for Driver.
Done simulating trajectory traj
False
Plotting the results#
sol = om.CaseReader('dymos_solution.db').get_case('final')
sim = om.CaseReader('dymos_simulation.db').get_case('final')
#############################################
# Plot the results
#############################################
rad = p.get_val('radius', units='m')[0]
print(f'optimal radius: {rad} m ')
mass = p.get_val('size_comp.mass', units='kg')[0]
print(f'cannonball mass: {mass} kg ')
area = p.get_val('size_comp.S', units='cm**2')[0]
print(f'cannonball aerodynamic reference area: {area} cm**2 ')
angle = p.get_val('traj.ascent.timeseries.gam', units='deg')[0, 0]
print(f'launch angle: {angle} deg')
max_range = p.get_val('traj.descent.timeseries.r')[-1, 0]
print(f'maximum range: {max_range} m')
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
time_imp = {'ascent': p.get_val('traj.ascent.timeseries.time'),
'descent': p.get_val('traj.descent.timeseries.time')}
time_exp = {'ascent': sim.get_val('traj.ascent.timeseries.time'),
'descent': sim.get_val('traj.descent.timeseries.time')}
r_imp = {'ascent': p.get_val('traj.ascent.timeseries.r'),
'descent': p.get_val('traj.descent.timeseries.r')}
r_exp = {'ascent': sim.get_val('traj.ascent.timeseries.r'),
'descent': sim.get_val('traj.descent.timeseries.r')}
h_imp = {'ascent': p.get_val('traj.ascent.timeseries.h'),
'descent': p.get_val('traj.descent.timeseries.h')}
h_exp = {'ascent': sim.get_val('traj.ascent.timeseries.h'),
'descent': sim.get_val('traj.descent.timeseries.h')}
axes.plot(r_imp['ascent'], h_imp['ascent'], 'bo')
axes.plot(r_imp['descent'], h_imp['descent'], 'ro')
axes.plot(r_exp['ascent'], h_exp['ascent'], 'b--')
axes.plot(r_exp['descent'], h_exp['descent'], 'r--')
axes.set_xlabel('range (m)')
axes.set_ylabel('altitude (m)')
axes.grid(True)
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
states = ['r', 'h', 'v', 'gam']
for i, state in enumerate(states):
x_imp = {'ascent': sol.get_val(f'traj.ascent.timeseries.{state}'),
'descent': sol.get_val(f'traj.descent.timeseries.{state}')}
x_exp = {'ascent': sim.get_val(f'traj.ascent.timeseries.{state}'),
'descent': sim.get_val(f'traj.descent.timeseries.{state}')}
axes[i].set_ylabel(state)
axes[i].grid(True)
axes[i].plot(time_imp['ascent'], x_imp['ascent'], 'bo')
axes[i].plot(time_imp['descent'], x_imp['descent'], 'ro')
axes[i].plot(time_exp['ascent'], x_exp['ascent'], 'b--')
axes[i].plot(time_exp['descent'], x_exp['descent'], 'r--')
plt.show()
optimal radius: 0.04185302553823411 m
cannonball mass: 2.416817833884683 kg
cannonball aerodynamic reference area: 55.03051657316943 cm**2
launch angle: 32.01951039452202 deg
maximum range: 3183.1177268008946 m
<Figure size 1000x600 with 1 Axes>
<Figure size 1000x600 with 4 Axes>