Cannonball Implicit Duration#
This example demonstrates determining the duration of a phase using a nonlinear solver to satisfy an endpoint condition. As in the multi-phase cannonball problem, we will simultaneously find the optimal design for the cannonball (its radius) and the optimal flight profile (its launch angle). However, here we will use a single phase that terminates at the condition
(91)#\[\begin{align}
h(t_f) = 0
\end{align}\]
The cannonball sizing component#
This component 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
Generating an initial guess#
Since we will be using a nonlinear solver to converge the dynamics and the duration, a simple linear initial guess will not be sufficient. The initial guess we will use is generated from the kinematics equations that assume no atmospheric drag. We compute the state values for some duration and initial flight path angle
def initial_guess(t_dur, gam0=np.pi/3):
t = np.linspace(0, t_dur, num=100)
g = -9.80665
v0 = -0.5*g*t_dur/np.sin(gam0)
r = v0*np.cos(gam0)*t
h = v0*np.sin(gam0)*t + 0.5*g*t**2
v = np.sqrt(v0*np.cos(gam0)**2 + (v0*np.sin(gam0) + g*t)**2)
gam = np.arctan2(v0*np.sin(gam0) + g*t, v0*np.cos(gam0)**2)
guess = {'t': t, 'r': r, 'h': h, 'v': v, 'gam': gam}
return guess
Building and running the problem#
The following code defines the components for the physical cannonball calculations and ODE problem, sets up trajectory using a single phase, and specifies the stopping condition for the phase. 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'] = 'IPOPT'
p.driver.declare_coloring()
p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'monotone'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.set_solver_print(level=0, depth=99)
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=25, order=3, compressed=False)
phase = dm.Phase(ode_class=CannonballODE, transcription=transcription)
phase = traj.add_phase('phase', phase)
# Duration is bounded to be greater than 1 to ensure the nonlinear solve
# does not converge to the initial point.
phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s')
# All initial states except flight path angle are fixed
# 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.
phase.add_state('r', fix_initial=True, solve_segments='forward')
phase.add_state('h', fix_initial=True, solve_segments='forward')
phase.add_state('gam', fix_initial=False, solve_segments='forward')
phase.add_state('v', fix_initial=False, solve_segments='forward')
phase.add_parameter('S', units='m**2', static_target=True)
phase.add_parameter('m', units='kg', static_target=True)
phase.add_parameter('CD', val=0.5, opt=False, static_target=True)
phase.add_boundary_constraint('ke', loc='initial',
upper=400000, lower=0, ref=100000)
# A duration balance is added setting altitude to zero.
# A nonlinear solver is used to find the duration of required to satisfy the condition.
phase.set_duration_balance('h', val=0.0)
phase.add_objective('r', loc='final', scaler=-1.0)
p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
p.model.connect('size_comp.S', 'traj.phase.parameters:S')
# 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.phase.parameters:CD', 0.5)
guess = initial_guess(t_dur=30.0)
p.set_val('traj.phase.t_initial', 0.0)
p.set_val('traj.phase.t_duration', 30.0)
p.set_val('traj.phase.states:r', phase.interp('r', ys=guess['r'], xs=guess['t']))
p.set_val('traj.phase.states:h', phase.interp('h', ys=guess['h'], xs=guess['t']))
p.set_val('traj.phase.states:v', phase.interp('v', ys=guess['v'], xs=guess['t']))
p.set_val('traj.phase.states:gam', phase.interp('gam', ys=guess['gam'], xs=guess['t']), units='rad')
#####################################################
# Run the optimization and final explicit simulation
#####################################################
dm.run_problem(p)
--- Constraint Report [traj] ---
--- phase ---
[initial] 0.0000e+00 <= ke <= 4.0000e+05 [J]
Model viewer data has already been recorded for Driver.
Model viewer data has already been recorded for Driver.
Coloring for 'traj.phases.phase.rhs_all' (class CannonballODE)
Jacobian shape: (500, 303) ( 0.92% nonzero)
FWD solves: 6 REV solves: 0
Total colors vs. total size: 6 vs 303 (98.0% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.144650 sec.
Time to compute coloring: 0.049749 sec.
Memory to compute coloring: 0.000000 MB.
Full total jacobian was computed 3 times, taking 0.447014 seconds.
Total jacobian shape: (98, 100)
Jacobian shape: (98, 100) (18.76% nonzero)
FWD solves: 21 REV solves: 0
Total colors vs. total size: 21 vs 100 (79.0% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.447014 sec.
Time to compute coloring: 0.104608 sec.
Memory to compute coloring: 2.621094 MB.
Optimization Problem -- Optimization using pyOpt_sparse
================================================================================
Objective Function: _objfunc
Solution:
--------------------------------------------------------------------------------
Total Time: 27.7786
User Objective Time : 24.1241
User Sensitivity Time : 3.1751
Interface Time : 0.2535
Opt Solver Time: 0.2259
Calls to Objective Function : 233
Calls to Sens Function : 54
Objectives
Index Name Value
0 traj.phases.phase.indep_states.states:r -3.182962E+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.540939E-01 1.000000E+00
1 traj.phase.t_duration_0 c 1.000000E+00 1.000000E+02 1.000000E+02 u
2 traj.phases.phase.indep_states.states:r_0 c -1.000000E+21 4.358669E+02 1.000000E+21
3 traj.phases.phase.indep_states.states:r_1 c -1.000000E+21 7.590741E+02 1.000000E+21
4 traj.phases.phase.indep_states.states:r_2 c -1.000000E+21 1.017857E+03 1.000000E+21
5 traj.phases.phase.indep_states.states:r_3 c -1.000000E+21 1.234831E+03 1.000000E+21
6 traj.phases.phase.indep_states.states:r_4 c -1.000000E+21 1.422453E+03 1.000000E+21
7 traj.phases.phase.indep_states.states:r_5 c -1.000000E+21 1.588316E+03 1.000000E+21
8 traj.phases.phase.indep_states.states:r_6 c -1.000000E+21 1.737381E+03 1.000000E+21
9 traj.phases.phase.indep_states.states:r_7 c -1.000000E+21 1.873054E+03 1.000000E+21
10 traj.phases.phase.indep_states.states:r_8 c -1.000000E+21 1.997754E+03 1.000000E+21
11 traj.phases.phase.indep_states.states:r_9 c -1.000000E+21 2.113240E+03 1.000000E+21
12 traj.phases.phase.indep_states.states:r_10 c -1.000000E+21 2.220813E+03 1.000000E+21
13 traj.phases.phase.indep_states.states:r_11 c -1.000000E+21 2.321446E+03 1.000000E+21
14 traj.phases.phase.indep_states.states:r_12 c -1.000000E+21 2.415862E+03 1.000000E+21
15 traj.phases.phase.indep_states.states:r_13 c -1.000000E+21 2.504609E+03 1.000000E+21
16 traj.phases.phase.indep_states.states:r_14 c -1.000000E+21 2.588097E+03 1.000000E+21
17 traj.phases.phase.indep_states.states:r_15 c -1.000000E+21 2.666639E+03 1.000000E+21
18 traj.phases.phase.indep_states.states:r_16 c -1.000000E+21 2.740483E+03 1.000000E+21
19 traj.phases.phase.indep_states.states:r_17 c -1.000000E+21 2.809828E+03 1.000000E+21
20 traj.phases.phase.indep_states.states:r_18 c -1.000000E+21 2.874850E+03 1.000000E+21
21 traj.phases.phase.indep_states.states:r_19 c -1.000000E+21 2.935708E+03 1.000000E+21
22 traj.phases.phase.indep_states.states:r_20 c -1.000000E+21 2.992555E+03 1.000000E+21
23 traj.phases.phase.indep_states.states:r_21 c -1.000000E+21 3.045544E+03 1.000000E+21
24 traj.phases.phase.indep_states.states:r_22 c -1.000000E+21 3.094833E+03 1.000000E+21
25 traj.phases.phase.indep_states.states:r_23 c -1.000000E+21 3.140584E+03 1.000000E+21
26 traj.phases.phase.indep_states.states:h_0 c -1.000000E+21 2.675246E+02 1.000000E+21
27 traj.phases.phase.indep_states.states:h_1 c -1.000000E+21 4.559862E+02 1.000000E+21
28 traj.phases.phase.indep_states.states:h_2 c -1.000000E+21 5.966683E+02 1.000000E+21
29 traj.phases.phase.indep_states.states:h_3 c -1.000000E+21 7.042126E+02 1.000000E+21
30 traj.phases.phase.indep_states.states:h_4 c -1.000000E+21 7.866604E+02 1.000000E+21
31 traj.phases.phase.indep_states.states:h_5 c -1.000000E+21 8.488963E+02 1.000000E+21
32 traj.phases.phase.indep_states.states:h_6 c -1.000000E+21 8.941007E+02 1.000000E+21
33 traj.phases.phase.indep_states.states:h_7 c -1.000000E+21 9.244542E+02 1.000000E+21
34 traj.phases.phase.indep_states.states:h_8 c -1.000000E+21 9.415152E+02 1.000000E+21
35 traj.phases.phase.indep_states.states:h_9 c -1.000000E+21 9.464412E+02 1.000000E+21
36 traj.phases.phase.indep_states.states:h_10 c -1.000000E+21 9.401270E+02 1.000000E+21
37 traj.phases.phase.indep_states.states:h_11 c -1.000000E+21 9.232973E+02 1.000000E+21
38 traj.phases.phase.indep_states.states:h_12 c -1.000000E+21 8.965701E+02 1.000000E+21
39 traj.phases.phase.indep_states.states:h_13 c -1.000000E+21 8.605033E+02 1.000000E+21
40 traj.phases.phase.indep_states.states:h_14 c -1.000000E+21 8.156262E+02 1.000000E+21
41 traj.phases.phase.indep_states.states:h_15 c -1.000000E+21 7.624607E+02 1.000000E+21
42 traj.phases.phase.indep_states.states:h_16 c -1.000000E+21 7.015341E+02 1.000000E+21
43 traj.phases.phase.indep_states.states:h_17 c -1.000000E+21 6.333833E+02 1.000000E+21
44 traj.phases.phase.indep_states.states:h_18 c -1.000000E+21 5.585557E+02 1.000000E+21
45 traj.phases.phase.indep_states.states:h_19 c -1.000000E+21 4.776046E+02 1.000000E+21
46 traj.phases.phase.indep_states.states:h_20 c -1.000000E+21 3.910839E+02 1.000000E+21
47 traj.phases.phase.indep_states.states:h_21 c -1.000000E+21 2.995404E+02 1.000000E+21
48 traj.phases.phase.indep_states.states:h_22 c -1.000000E+21 2.035072E+02 1.000000E+21
49 traj.phases.phase.indep_states.states:h_23 c -1.000000E+21 1.034977E+02 1.000000E+21
50 traj.phases.phase.indep_states.states:gam_0 c -1.000000E+21 5.588424E-01 1.000000E+21
51 traj.phases.phase.indep_states.states:gam_1 c -1.000000E+21 5.398186E-01 1.000000E+21
52 traj.phases.phase.indep_states.states:gam_2 c -1.000000E+21 5.135922E-01 1.000000E+21
53 traj.phases.phase.indep_states.states:gam_3 c -1.000000E+21 4.797863E-01 1.000000E+21
54 traj.phases.phase.indep_states.states:gam_4 c -1.000000E+21 4.379029E-01 1.000000E+21
55 traj.phases.phase.indep_states.states:gam_5 c -1.000000E+21 3.873723E-01 1.000000E+21
56 traj.phases.phase.indep_states.states:gam_6 c -1.000000E+21 3.276275E-01 1.000000E+21
57 traj.phases.phase.indep_states.states:gam_7 c -1.000000E+21 2.582142E-01 1.000000E+21
58 traj.phases.phase.indep_states.states:gam_8 c -1.000000E+21 1.789451E-01 1.000000E+21
59 traj.phases.phase.indep_states.states:gam_9 c -1.000000E+21 9.009145E-02 1.000000E+21
60 traj.phases.phase.indep_states.states:gam_10 c -1.000000E+21 -7.423019E-03 1.000000E+21
61 traj.phases.phase.indep_states.states:gam_11 c -1.000000E+21 -1.118970E-01 1.000000E+21
62 traj.phases.phase.indep_states.states:gam_12 c -1.000000E+21 -2.208805E-01 1.000000E+21
63 traj.phases.phase.indep_states.states:gam_13 c -1.000000E+21 -3.314120E-01 1.000000E+21
64 traj.phases.phase.indep_states.states:gam_14 c -1.000000E+21 -4.404248E-01 1.000000E+21
65 traj.phases.phase.indep_states.states:gam_15 c -1.000000E+21 -5.451909E-01 1.000000E+21
66 traj.phases.phase.indep_states.states:gam_16 c -1.000000E+21 -6.436480E-01 1.000000E+21
67 traj.phases.phase.indep_states.states:gam_17 c -1.000000E+21 -7.345244E-01 1.000000E+21
68 traj.phases.phase.indep_states.states:gam_18 c -1.000000E+21 -8.172791E-01 1.000000E+21
69 traj.phases.phase.indep_states.states:gam_19 c -1.000000E+21 -8.919360E-01 1.000000E+21
70 traj.phases.phase.indep_states.states:gam_20 c -1.000000E+21 -9.588926E-01 1.000000E+21
71 traj.phases.phase.indep_states.states:gam_21 c -1.000000E+21 -1.018753E+00 1.000000E+21
72 traj.phases.phase.indep_states.states:gam_22 c -1.000000E+21 -1.072207E+00 1.000000E+21
73 traj.phases.phase.indep_states.states:gam_23 c -1.000000E+21 -1.119950E+00 1.000000E+21
74 traj.phases.phase.indep_states.states:gam_24 c -1.000000E+21 -1.162642E+00 1.000000E+21
75 traj.phases.phase.indep_states.states:v_0 c -1.000000E+21 5.750200E+02 1.000000E+21
76 traj.phases.phase.indep_states.states:v_1 c -1.000000E+21 3.997227E+02 1.000000E+21
77 traj.phases.phase.indep_states.states:v_2 c -1.000000E+21 3.060109E+02 1.000000E+21
78 traj.phases.phase.indep_states.states:v_3 c -1.000000E+21 2.471802E+02 1.000000E+21
79 traj.phases.phase.indep_states.states:v_4 c -1.000000E+21 2.066205E+02 1.000000E+21
80 traj.phases.phase.indep_states.states:v_5 c -1.000000E+21 1.769413E+02 1.000000E+21
81 traj.phases.phase.indep_states.states:v_6 c -1.000000E+21 1.543746E+02 1.000000E+21
82 traj.phases.phase.indep_states.states:v_7 c -1.000000E+21 1.368146E+02 1.000000E+21
83 traj.phases.phase.indep_states.states:v_8 c -1.000000E+21 1.230084E+02 1.000000E+21
84 traj.phases.phase.indep_states.states:v_9 c -1.000000E+21 1.121738E+02 1.000000E+21
85 traj.phases.phase.indep_states.states:v_10 c -1.000000E+21 1.037991E+02 1.000000E+21
86 traj.phases.phase.indep_states.states:v_11 c -1.000000E+21 9.752646E+01 1.000000E+21
87 traj.phases.phase.indep_states.states:v_12 c -1.000000E+21 9.307779E+01 1.000000E+21
88 traj.phases.phase.indep_states.states:v_13 c -1.000000E+21 9.020964E+01 1.000000E+21
89 traj.phases.phase.indep_states.states:v_14 c -1.000000E+21 8.868763E+01 1.000000E+21
90 traj.phases.phase.indep_states.states:v_15 c -1.000000E+21 8.827834E+01 1.000000E+21
91 traj.phases.phase.indep_states.states:v_16 c -1.000000E+21 8.875205E+01 1.000000E+21
92 traj.phases.phase.indep_states.states:v_17 c -1.000000E+21 8.989109E+01 1.000000E+21
93 traj.phases.phase.indep_states.states:v_18 c -1.000000E+21 9.149866E+01 1.000000E+21
94 traj.phases.phase.indep_states.states:v_19 c -1.000000E+21 9.340455E+01 1.000000E+21
95 traj.phases.phase.indep_states.states:v_20 c -1.000000E+21 9.546736E+01 1.000000E+21
96 traj.phases.phase.indep_states.states:v_21 c -1.000000E+21 9.757524E+01 1.000000E+21
97 traj.phases.phase.indep_states.states:v_22 c -1.000000E+21 9.964293E+01 1.000000E+21
98 traj.phases.phase.indep_states.states:v_23 c -1.000000E+21 1.016075E+02 1.000000E+21
99 traj.phases.phase.indep_states.states:v_24 c -1.000000E+21 1.034255E+02 1.000000E+21
Constraints (i - inequality, e - equality)
Index Name Type Lower Value Upper Status Lagrange Multiplier (N/A)
0 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 4.263256E-12 0.000000E+00 9.00000E+100
1 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -3.637979E-12 0.000000E+00 9.00000E+100
2 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -4.320100E-12 0.000000E+00 9.00000E+100
3 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -3.637979E-12 0.000000E+00 9.00000E+100
4 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -2.728484E-12 0.000000E+00 9.00000E+100
5 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -2.501110E-12 0.000000E+00 9.00000E+100
6 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -1.591616E-12 0.000000E+00 9.00000E+100
7 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -1.364242E-12 0.000000E+00 9.00000E+100
8 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -6.821210E-13 0.000000E+00 9.00000E+100
9 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -1.364242E-12 0.000000E+00 9.00000E+100
10 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -4.547474E-13 0.000000E+00 9.00000E+100
11 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 4.547474E-13 0.000000E+00 9.00000E+100
12 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 4.547474E-13 0.000000E+00 9.00000E+100
13 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 4.547474E-13 0.000000E+00 9.00000E+100
14 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 1.818989E-12 0.000000E+00 9.00000E+100
15 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 -4.547474E-13 0.000000E+00 9.00000E+100
16 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 1.364242E-12 0.000000E+00 9.00000E+100
17 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 2.273737E-12 0.000000E+00 9.00000E+100
18 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 9.094947E-13 0.000000E+00 9.00000E+100
19 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 2.273737E-12 0.000000E+00 9.00000E+100
20 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 2.728484E-12 0.000000E+00 9.00000E+100
21 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 1.364242E-12 0.000000E+00 9.00000E+100
22 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 1.818989E-12 0.000000E+00 9.00000E+100
23 traj.phases.phase.continuity_comp.defect_states:r e 0.000000E+00 2.273737E-12 0.000000E+00 9.00000E+100
24 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 -2.211209E-11 0.000000E+00 9.00000E+100
25 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 -1.273293E-11 0.000000E+00 9.00000E+100
26 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 -7.503331E-12 0.000000E+00 9.00000E+100
27 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 -4.092726E-12 0.000000E+00 9.00000E+100
28 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 -1.705303E-12 0.000000E+00 9.00000E+100
29 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 -4.547474E-13 0.000000E+00 9.00000E+100
30 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 5.684342E-13 0.000000E+00 9.00000E+100
31 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.410605E-13 0.000000E+00 9.00000E+100
32 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 6.821210E-13 0.000000E+00 9.00000E+100
33 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.410605E-13 0.000000E+00 9.00000E+100
34 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 4.547474E-13 0.000000E+00 9.00000E+100
35 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 2.273737E-13 0.000000E+00 9.00000E+100
36 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 0.000000E+00 0.000000E+00 9.00000E+100
37 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.410605E-13 0.000000E+00 9.00000E+100
38 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.410605E-13 0.000000E+00 9.00000E+100
39 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 1.250555E-12 0.000000E+00 9.00000E+100
40 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 1.477929E-12 0.000000E+00 9.00000E+100
41 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 1.932676E-12 0.000000E+00 9.00000E+100
42 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 2.614797E-12 0.000000E+00 9.00000E+100
43 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.069545E-12 0.000000E+00 9.00000E+100
44 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.183231E-12 0.000000E+00 9.00000E+100
45 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.353762E-12 0.000000E+00 9.00000E+100
46 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.609557E-12 0.000000E+00 9.00000E+100
47 traj.phases.phase.continuity_comp.defect_states:h e 0.000000E+00 3.552714E-12 0.000000E+00 9.00000E+100
48 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.665335E-15 0.000000E+00 9.00000E+100
49 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 8.881784E-16 0.000000E+00 9.00000E+100
50 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.387779E-15 0.000000E+00 9.00000E+100
51 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.998401E-15 0.000000E+00 9.00000E+100
52 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 3.275158E-15 0.000000E+00 9.00000E+100
53 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 4.718448E-15 0.000000E+00 9.00000E+100
54 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 7.382983E-15 0.000000E+00 9.00000E+100
55 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 9.686696E-15 0.000000E+00 9.00000E+100
56 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.179612E-14 0.000000E+00 9.00000E+100
57 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.316048E-14 0.000000E+00 9.00000E+100
58 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.294798E-14 0.000000E+00 9.00000E+100
59 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.143530E-14 0.000000E+00 9.00000E+100
60 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 8.437695E-15 0.000000E+00 9.00000E+100
61 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 4.940492E-15 0.000000E+00 9.00000E+100
62 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 1.554312E-15 0.000000E+00 9.00000E+100
63 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -6.661338E-16 0.000000E+00 9.00000E+100
64 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -2.109424E-15 0.000000E+00 9.00000E+100
65 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -2.553513E-15 0.000000E+00 9.00000E+100
66 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -2.886580E-15 0.000000E+00 9.00000E+100
67 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -2.664535E-15 0.000000E+00 9.00000E+100
68 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -2.886580E-15 0.000000E+00 9.00000E+100
69 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -2.220446E-15 0.000000E+00 9.00000E+100
70 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -1.998401E-15 0.000000E+00 9.00000E+100
71 traj.phases.phase.continuity_comp.defect_states:gam e 0.000000E+00 -1.776357E-15 0.000000E+00 9.00000E+100
72 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 2.512479E-11 0.000000E+00 9.00000E+100
73 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 8.355983E-12 0.000000E+00 9.00000E+100
74 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 3.552714E-12 0.000000E+00 9.00000E+100
75 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 1.762146E-12 0.000000E+00 9.00000E+100
76 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 9.663381E-13 0.000000E+00 9.00000E+100
77 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 5.968559E-13 0.000000E+00 9.00000E+100
78 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 4.547474E-13 0.000000E+00 9.00000E+100
79 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 1.989520E-13 0.000000E+00 9.00000E+100
80 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 7.105427E-14 0.000000E+00 9.00000E+100
81 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -5.684342E-14 0.000000E+00 9.00000E+100
82 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -2.131628E-13 0.000000E+00 9.00000E+100
83 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -2.984279E-13 0.000000E+00 9.00000E+100
84 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -3.410605E-13 0.000000E+00 9.00000E+100
85 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -3.552714E-13 0.000000E+00 9.00000E+100
86 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -2.557954E-13 0.000000E+00 9.00000E+100
87 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -1.705303E-13 0.000000E+00 9.00000E+100
88 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 -1.136868E-13 0.000000E+00 9.00000E+100
89 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 0.000000E+00 0.000000E+00 9.00000E+100
90 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 1.278977E-13 0.000000E+00 9.00000E+100
91 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 1.847411E-13 0.000000E+00 9.00000E+100
92 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 2.415845E-13 0.000000E+00 9.00000E+100
93 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 3.268497E-13 0.000000E+00 9.00000E+100
94 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 3.552714E-13 0.000000E+00 9.00000E+100
95 traj.phases.phase.continuity_comp.defect_states:v e 0.000000E+00 4.405365E-13 0.000000E+00 9.00000E+100
96 traj.phases.phase->initial_boundary_constraint->ke i 0.000000E+00 4.000000E+00 4.000000E+00 u 9.00000E+100
Exit Status
Inform Description
0 Solve Succeeded
--------------------------------------------------------------------------------
False
Plotting the results#
exp_out = traj.simulate()
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 ')
angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]
print(f'launch angle: {angle} deg')
max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]
print(f'maximum range: {max_range} m')
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
time_imp = p.get_val('traj.phase.timeseries.time')
time_exp = exp_out.get_val('traj.phase.timeseries.time')
r_imp = p.get_val('traj.phase.timeseries.r')
r_exp = exp_out.get_val('traj.phase.timeseries.r')
h_imp = p.get_val('traj.phase.timeseries.h')
h_exp = exp_out.get_val('traj.phase.timeseries.h')
axes.plot(r_imp, h_imp, 'bo')
axes.plot(r_exp, h_exp, 'r--')
axes.set_xlabel('range (m)')
axes.set_ylabel('altitude (m)')
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
states = ['r', 'h', 'v', 'gam']
for i, state in enumerate(states):
x_imp = p.get_val(f'traj.phase.timeseries.{state}')
x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')
axes[i].set_ylabel(state)
axes[i].plot(time_imp, x_imp, 'bo')
axes[i].plot(time_exp, x_exp, 'r--')
plt.show()
Simulating trajectory traj
Done simulating trajectory traj
optimal radius: 0.041868454738899524 m
cannonball mass: 2.4194917134009 kg
launch angle: 32.01931307429335 deg
maximum range: 3182.9624372149538 m
