Skip to content

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

\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

\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

\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

\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.

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.

import numpy as np
from scipy.interpolate import interp1d

import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal

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

#############################################
# Build the ODE class
#############################################
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=['state_rate_source:v'])
        self.add_output('gam_dot', shape=nn, units='rad/s', tags=['state_rate_source:gam'])
        self.add_output('h_dot', shape=nn, units='m/s', tags=['state_rate_source:h'])
        self.add_output('r_dot', shape=nn, units='m/s', tags=['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

#############################################
# Setup the Dymos problem
#############################################

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', dynamic=False)
ascent.add_parameter('m', units='kg', dynamic=False)

# 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', dynamic=False)
descent.add_parameter('m', units='kg', dynamic=False)

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, dynamic=False)

# 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': 'mass', 'descent': 'mass'}, dynamic=False)

# 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, dynamic=False)

# 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.interpolate(ys=[0, 100],
          nodes='state_input'))
p.set_val('traj.ascent.states:h', ascent.interpolate(ys=[0, 100],
          nodes='state_input'))
p.set_val('traj.ascent.states:v', ascent.interpolate(ys=[200, 150],
          nodes='state_input'))
p.set_val('traj.ascent.states:gam', ascent.interpolate(ys=[25, 0],
          nodes='state_input'), 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.interpolate(ys=[100, 200],
          nodes='state_input'))
p.set_val('traj.descent.states:h', descent.interpolate(ys=[100, 0],
          nodes='state_input'))
p.set_val('traj.descent.states:v', descent.interpolate(ys=[150, 200],
          nodes='state_input'))
p.set_val('traj.descent.states:gam', descent.interpolate(ys=[0, -45],
          nodes='state_input'), units='deg')

#####################################################
# Run the optimization and final explicit simulation
#####################################################
dm.run_problem(p)



# use the explicit simulation to check the final collocation solution accuracy
exp_out = traj.simulate()

#############################################
# 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 ')
angle = p.get_val('traj.ascent.timeseries.states:gam', units='deg')[0, 0]
print(f'launch angle: {angle} deg')
max_range = p.get_val('traj.descent.timeseries.states: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': exp_out.get_val('traj.ascent.timeseries.time'),
            'descent': exp_out.get_val('traj.descent.timeseries.time')}

r_imp = {'ascent': p.get_val('traj.ascent.timeseries.states:r'),
         'descent': p.get_val('traj.descent.timeseries.states:r')}

r_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.states:r'),
         'descent': exp_out.get_val('traj.descent.timeseries.states:r')}

h_imp = {'ascent': p.get_val('traj.ascent.timeseries.states:h'),
         'descent': p.get_val('traj.descent.timeseries.states:h')}

h_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.states:h'),
         'descent': exp_out.get_val('traj.descent.timeseries.states: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)')

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': p.get_val(f'traj.ascent.timeseries.states:{state}'),
             'descent': p.get_val(f'traj.descent.timeseries.states:{state}')}

    x_exp = {'ascent': exp_out.get_val(f'traj.ascent.timeseries.states:{state}'),
             'descent': exp_out.get_val(f'traj.descent.timeseries.states:{state}')}

    axes[i].set_ylabel(state)

    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--')

params = ['m', 'S']
fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(12, 6))
for i, param in enumerate(params):
    p_imp = {
        'ascent': p.get_val(f'traj.ascent.timeseries.parameters:{param}'),
        'descent': p.get_val(f'traj.descent.timeseries.parameters:{param}')}

    p_exp = {'ascent': exp_out.get_val(f'traj.ascent.timeseries.parameters:{param}'),
             'descent': exp_out.get_val(f'traj.descent.timeseries.parameters:{param}')}

    axes[i].set_ylabel(param)

    axes[i].plot(time_imp['ascent'], p_imp['ascent'], 'bo')
    axes[i].plot(time_imp['descent'], p_imp['descent'], 'ro')
    axes[i].plot(time_exp['ascent'], p_exp['ascent'], 'b--')
    axes[i].plot(time_exp['descent'], p_exp['descent'], 'r--')

plt.show()
--- Linkage Report [traj] ---
    --- ascent - descent ---
        time [final]  ==  time [initial]
        r    [final]  ==  r    [initial]
        h    [final]  ==  h    [initial]
        gam  [final]  ==  gam  [initial]
        v    [final]  ==  v    [initial]
----------------------------

Approx 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.5% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.033778 sec.
Time to compute coloring: 0.001997 sec.

Approx 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.8% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.019570 sec.
Time to compute coloring: 0.001112 sec.

Approx 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.7% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.012182 sec.
Time to compute coloring: 0.000723 sec.
Full total jacobian was computed 3 times, taking 0.600696 seconds.
Total jacobian shape: (87, 88) 


Jacobian shape: (87, 88)  (13.98% nonzero)
FWD solves: 27   REV solves: 0
Total colors vs. total size: 27 vs 88  (69.3% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.600696 sec.
Time to compute coloring: 0.004195 sec.


Optimization Problem -- Optimization using pyOpt_sparse
================================================================================
    Objective Function: _objfunc

    Solution: 
--------------------------------------------------------------------------------
    Total Time:                   10.7333
       User Objective Time :       0.5899
       User Sensitivity Time :     9.0918
       Interface Time :            0.8561
       Opt Solver Time:            0.1955
    Calls to Objective Function :     118
    Calls to Sens Function :           89


   Objectives
      Index  Name                                                 Value          Optimum
          0  traj.phases.descent.indep_states.states:r    -3.183118E+03     0.000000E+00

   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.phases.ascent.time_extents.t_duration_0       c     1.000000E-02     1.064895E-01     1.000000E+00           
          2  traj.phases.ascent.indep_states.states:r_0         c    -1.000000E+21     3.231043E+02     1.000000E+21           
          3  traj.phases.ascent.indep_states.states:r_1         c    -1.000000E+21     6.616876E+02     1.000000E+21           
          4  traj.phases.ascent.indep_states.states:r_2         c    -1.000000E+21     7.552526E+02     1.000000E+21           
          5  traj.phases.ascent.indep_states.states:r_3         c    -1.000000E+21     9.434454E+02     1.000000E+21           
          6  traj.phases.ascent.indep_states.states:r_4         c    -1.000000E+21     1.165633E+03     1.000000E+21           
          7  traj.phases.ascent.indep_states.states:r_5         c    -1.000000E+21     1.229354E+03     1.000000E+21           
          8  traj.phases.ascent.indep_states.states:r_6         c    -1.000000E+21     1.364715E+03     1.000000E+21           
          9  traj.phases.ascent.indep_states.states:r_7         c    -1.000000E+21     1.532445E+03     1.000000E+21           
         10  traj.phases.ascent.indep_states.states:r_8         c    -1.000000E+21     1.581810E+03     1.000000E+21           
         11  traj.phases.ascent.indep_states.states:r_9         c    -1.000000E+21     1.688890E+03     1.000000E+21           
         12  traj.phases.ascent.indep_states.states:r_10        c    -1.000000E+21     1.825079E+03     1.000000E+21           
         13  traj.phases.ascent.indep_states.states:r_11        c    -1.000000E+21     1.865797E+03     1.000000E+21           
         14  traj.phases.ascent.indep_states.states:r_12        c    -1.000000E+21     1.955144E+03     1.000000E+21           
         15  traj.phases.ascent.indep_states.states:r_13        c    -1.000000E+21     2.070557E+03     1.000000E+21           
         16  traj.phases.ascent.indep_states.states:r_14        c    -1.000000E+21     2.105398E+03     1.000000E+21           
         17  traj.phases.ascent.indep_states.states:h_0         c    -1.000000E+21     1.995146E+02     1.000000E+21           
         18  traj.phases.ascent.indep_states.states:h_1         c    -1.000000E+21     4.002917E+02     1.000000E+21           
         19  traj.phases.ascent.indep_states.states:h_2         c    -1.000000E+21     4.537845E+02     1.000000E+21           
         20  traj.phases.ascent.indep_states.states:h_3         c    -1.000000E+21     5.574056E+02     1.000000E+21           
         21  traj.phases.ascent.indep_states.states:h_4         c    -1.000000E+21     6.711526E+02     1.000000E+21           
         22  traj.phases.ascent.indep_states.states:h_5         c    -1.000000E+21     7.016712E+02     1.000000E+21           
         23  traj.phases.ascent.indep_states.states:h_6         c    -1.000000E+21     7.625774E+02     1.000000E+21           
         24  traj.phases.ascent.indep_states.states:h_7         c    -1.000000E+21     8.292929E+02     1.000000E+21           
         25  traj.phases.ascent.indep_states.states:h_8         c    -1.000000E+21     8.467852E+02     1.000000E+21           
         26  traj.phases.ascent.indep_states.states:h_9         c    -1.000000E+21     8.807923E+02     1.000000E+21           
         27  traj.phases.ascent.indep_states.states:h_10        c    -1.000000E+21     9.151926E+02     1.000000E+21           
         28  traj.phases.ascent.indep_states.states:h_11        c    -1.000000E+21     9.233143E+02     1.000000E+21           
         29  traj.phases.ascent.indep_states.states:h_12        c    -1.000000E+21     9.371859E+02     1.000000E+21           
         30  traj.phases.ascent.indep_states.states:h_13        c    -1.000000E+21     9.461925E+02     1.000000E+21           
         31  traj.phases.ascent.indep_states.states:h_14        c    -1.000000E+21     9.467365E+02     1.000000E+21           
         32  traj.phases.ascent.indep_states.states:gam_0       c    -1.000000E+21     5.588459E-01     1.000000E+21           
         33  traj.phases.ascent.indep_states.states:gam_1       c    -1.000000E+21     5.461378E-01     1.000000E+21           
         34  traj.phases.ascent.indep_states.states:gam_2       c    -1.000000E+21     5.228406E-01     1.000000E+21           
         35  traj.phases.ascent.indep_states.states:gam_3       c    -1.000000E+21     5.140981E-01     1.000000E+21           
         36  traj.phases.ascent.indep_states.states:gam_4       c    -1.000000E+21     4.911599E-01     1.000000E+21           
         37  traj.phases.ascent.indep_states.states:gam_5       c    -1.000000E+21     4.530243E-01     1.000000E+21           
         38  traj.phases.ascent.indep_states.states:gam_6       c    -1.000000E+21     4.393302E-01     1.000000E+21           
         39  traj.phases.ascent.indep_states.states:gam_7       c    -1.000000E+21     4.048232E-01     1.000000E+21           
         40  traj.phases.ascent.indep_states.states:gam_8       c    -1.000000E+21     3.498022E-01     1.000000E+21           
         41  traj.phases.ascent.indep_states.states:gam_9       c    -1.000000E+21     3.305300E-01     1.000000E+21           
         42  traj.phases.ascent.indep_states.states:gam_10      c    -1.000000E+21     2.829154E-01     1.000000E+21           
         43  traj.phases.ascent.indep_states.states:gam_11      c    -1.000000E+21     2.091969E-01     1.000000E+21           
         44  traj.phases.ascent.indep_states.states:gam_12      c    -1.000000E+21     1.839162E-01     1.000000E+21           
         45  traj.phases.ascent.indep_states.states:gam_13      c    -1.000000E+21     1.226062E-01     1.000000E+21           
         46  traj.phases.ascent.indep_states.states:gam_14      c    -1.000000E+21     3.070674E-02     1.000000E+21           
         47  traj.phases.ascent.indep_states.states:v_0         c    -1.000000E+21     5.753380E+02     1.000000E+21           
         48  traj.phases.ascent.indep_states.states:v_1         c    -1.000000E+21     4.377128E+02     1.000000E+21           
         49  traj.phases.ascent.indep_states.states:v_2         c    -1.000000E+21     3.333818E+02     1.000000E+21           
         50  traj.phases.ascent.indep_states.states:v_3         c    -1.000000E+21     3.071231E+02     1.000000E+21           
         51  traj.phases.ascent.indep_states.states:v_4         c    -1.000000E+21     2.628446E+02     1.000000E+21           
         52  traj.phases.ascent.indep_states.states:v_5         c    -1.000000E+21     2.189809E+02     1.000000E+21           
         53  traj.phases.ascent.indep_states.states:v_6         c    -1.000000E+21     2.076302E+02     1.000000E+21           
         54  traj.phases.ascent.indep_states.states:v_7         c    -1.000000E+21     1.856231E+02     1.000000E+21           
         55  traj.phases.ascent.indep_states.states:v_8         c    -1.000000E+21     1.616748E+02     1.000000E+21           
         56  traj.phases.ascent.indep_states.states:v_9         c    -1.000000E+21     1.552355E+02     1.000000E+21           
         57  traj.phases.ascent.indep_states.states:v_10        c    -1.000000E+21     1.422744E+02     1.000000E+21           
         58  traj.phases.ascent.indep_states.states:v_11        c    -1.000000E+21     1.276786E+02     1.000000E+21           
         59  traj.phases.ascent.indep_states.states:v_12        c    -1.000000E+21     1.237104E+02     1.000000E+21           
         60  traj.phases.ascent.indep_states.states:v_13        c    -1.000000E+21     1.156806E+02     1.000000E+21           
         61  traj.phases.ascent.indep_states.states:v_14        c    -1.000000E+21     1.067082E+02     1.000000E+21           
         62  traj.phases.ascent.indep_states.states:v_15        c    -1.000000E+21     1.043171E+02     1.000000E+21           
         63  traj.phases.descent.time_extents.t_initial_0       c     5.000000E-01     1.064895E+01     1.000000E+02           
         64  traj.phases.descent.time_extents.t_duration_0      c     5.000000E-03     1.617267E-01     1.000000E+00           
         65  traj.phases.descent.indep_states.states:r_0        c    -1.000000E+21     2.105398E+03     1.000000E+21           
         66  traj.phases.descent.indep_states.states:r_1        c    -1.000000E+21     2.410813E+03     1.000000E+21           
         67  traj.phases.descent.indep_states.states:r_2        c    -1.000000E+21     2.663597E+03     1.000000E+21           
         68  traj.phases.descent.indep_states.states:r_3        c    -1.000000E+21     2.873278E+03     1.000000E+21           
         69  traj.phases.descent.indep_states.states:r_4        c    -1.000000E+21     3.045011E+03     1.000000E+21           
         70  traj.phases.descent.indep_states.states:r_5        c    -1.000000E+21     3.183118E+03     1.000000E+21           
         71  traj.phases.descent.indep_states.states:h_0        c    -1.000000E+21     9.467365E+02     1.000000E+21           
         72  traj.phases.descent.indep_states.states:h_1        c    -1.000000E+21     8.986347E+02     1.000000E+21           
         73  traj.phases.descent.indep_states.states:h_2        c    -1.000000E+21     7.651714E+02     1.000000E+21           
         74  traj.phases.descent.indep_states.states:h_3        c    -1.000000E+21     5.609894E+02     1.000000E+21           
         75  traj.phases.descent.indep_states.states:h_4        c    -1.000000E+21     3.009916E+02     1.000000E+21           
         76  traj.phases.descent.indep_states.states:gam_0      c    -1.000000E+21    -1.567954E-20     1.000000E+21           
         77  traj.phases.descent.indep_states.states:gam_1      c    -1.000000E+21    -3.251038E-01     1.000000E+21           
         78  traj.phases.descent.indep_states.states:gam_2      c    -1.000000E+21    -6.397506E-01     1.000000E+21           
         79  traj.phases.descent.indep_states.states:gam_3      c    -1.000000E+21    -8.901606E-01     1.000000E+21           
         80  traj.phases.descent.indep_states.states:gam_4      c    -1.000000E+21    -1.071739E+00     1.000000E+21           
         81  traj.phases.descent.indep_states.states:gam_5      c    -1.000000E+21    -1.201117E+00     1.000000E+21           
         82  traj.phases.descent.indep_states.states:v_0        c    -1.000000E+21     1.043171E+02     1.000000E+21           
         83  traj.phases.descent.indep_states.states:v_1        c    -1.000000E+21     9.031921E+01     1.000000E+21           
         84  traj.phases.descent.indep_states.states:v_2        c    -1.000000E+21     8.870119E+01     1.000000E+21           
         85  traj.phases.descent.indep_states.states:v_3        c    -1.000000E+21     9.333723E+01     1.000000E+21           
         86  traj.phases.descent.indep_states.states:v_4        c    -1.000000E+21     9.960765E+01     1.000000E+21           
         87  traj.phases.descent.indep_states.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.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    1.355922E-11    0.000000E+00              9.00000E+100
          1  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    7.203338E-12    0.000000E+00              9.00000E+100
          2  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    3.631935E-12    0.000000E+00              9.00000E+100
          3  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    2.663419E-12    0.000000E+00              9.00000E+100
          4  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    1.271177E-12    0.000000E+00              9.00000E+100
          5  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -9.079838E-14    0.000000E+00              9.00000E+100
          6  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -8.171854E-13    0.000000E+00              9.00000E+100
          7  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    9.079838E-14    0.000000E+00              9.00000E+100
          8  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00    1.815968E-13    0.000000E+00              9.00000E+100
          9  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -6.961209E-13    0.000000E+00              9.00000E+100
         10  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -1.059314E-12    0.000000E+00              9.00000E+100
         11  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -9.079838E-13    0.000000E+00              9.00000E+100
         12  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -1.649504E-12    0.000000E+00              9.00000E+100
         13  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -6.658548E-13    0.000000E+00              9.00000E+100
         14  traj.phases.ascent.collocation_constraint.defects:r                 e   0.000000E+00   -8.777176E-13    0.000000E+00              9.00000E+100
         15  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    1.035101E-11    0.000000E+00              9.00000E+100
         16  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    3.450338E-12    0.000000E+00              9.00000E+100
         17  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    1.271177E-12    0.000000E+00              9.00000E+100
         18  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    4.237258E-13    0.000000E+00              9.00000E+100
         19  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    2.572621E-13    0.000000E+00              9.00000E+100
         20  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    1.815968E-13    0.000000E+00              9.00000E+100
         21  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    1.361976E-13    0.000000E+00              9.00000E+100
         22  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    4.691249E-13    0.000000E+00              9.00000E+100
         23  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    5.750564E-13    0.000000E+00              9.00000E+100
         24  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    7.339535E-13    0.000000E+00              9.00000E+100
         25  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    1.740302E-13    0.000000E+00              9.00000E+100
         26  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00   -2.269959E-14    0.000000E+00              9.00000E+100
         27  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    6.053225E-13    0.000000E+00              9.00000E+100
         28  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    3.404939E-14    0.000000E+00              9.00000E+100
         29  traj.phases.ascent.collocation_constraint.defects:h                 e   0.000000E+00    1.191729E-13    0.000000E+00              9.00000E+100
         30  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    7.725399E-15    0.000000E+00              9.00000E+100
         31  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    3.816517E-15    0.000000E+00              9.00000E+100
         32  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    2.475379E-15    0.000000E+00              9.00000E+100
         33  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    1.174881E-15    0.000000E+00              9.00000E+100
         34  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    6.502488E-16    0.000000E+00              9.00000E+100
         35  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    9.605948E-17    0.000000E+00              9.00000E+100
         36  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00   -2.216757E-16    0.000000E+00              9.00000E+100
         37  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    3.694595E-17    0.000000E+00              9.00000E+100
         38  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00   -1.256162E-16    0.000000E+00              9.00000E+100
         39  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    3.768487E-16    0.000000E+00              9.00000E+100
         40  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    2.364541E-16    0.000000E+00              9.00000E+100
         41  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    5.911353E-16    0.000000E+00              9.00000E+100
         42  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    9.310380E-16    0.000000E+00              9.00000E+100
         43  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    1.374389E-15    0.000000E+00              9.00000E+100
         44  traj.phases.ascent.collocation_constraint.defects:gam               e   0.000000E+00    1.995082E-15    0.000000E+00              9.00000E+100
         45  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    2.960632E-10    0.000000E+00              9.00000E+100
         46  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    9.724506E-11    0.000000E+00              9.00000E+100
         47  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    3.636475E-11    0.000000E+00              9.00000E+100
         48  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    2.489389E-11    0.000000E+00              9.00000E+100
         49  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    1.401322E-11    0.000000E+00              9.00000E+100
         50  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    7.150372E-12    0.000000E+00              9.00000E+100
         51  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    5.947294E-12    0.000000E+00              9.00000E+100
         52  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    3.957296E-12    0.000000E+00              9.00000E+100
         53  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    2.432640E-12    0.000000E+00              9.00000E+100
         54  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    2.137545E-12    0.000000E+00              9.00000E+100
         55  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    1.581405E-12    0.000000E+00              9.00000E+100
         56  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    1.104714E-12    0.000000E+00              9.00000E+100
         57  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    9.136587E-13    0.000000E+00              9.00000E+100
         58  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    7.944858E-13    0.000000E+00              9.00000E+100
         59  traj.phases.ascent.collocation_constraint.defects:v                 e   0.000000E+00    5.816771E-13    0.000000E+00              9.00000E+100
         60  traj.phases.descent.collocation_constraint.defects:r                e   0.000000E+00    8.319756E-12    0.000000E+00              9.00000E+100
         61  traj.phases.descent.collocation_constraint.defects:r                e   0.000000E+00    8.296773E-12    0.000000E+00              9.00000E+100
         62  traj.phases.descent.collocation_constraint.defects:r                e   0.000000E+00    8.618532E-12    0.000000E+00              9.00000E+100
         63  traj.phases.descent.collocation_constraint.defects:r                e   0.000000E+00    7.342989E-12    0.000000E+00              9.00000E+100
         64  traj.phases.descent.collocation_constraint.defects:r                e   0.000000E+00    6.619032E-12    0.000000E+00              9.00000E+100
         65  traj.phases.descent.collocation_constraint.defects:h                e   0.000000E+00   -1.953534E-12    0.000000E+00              9.00000E+100
         66  traj.phases.descent.collocation_constraint.defects:h                e   0.000000E+00   -5.079188E-12    0.000000E+00              9.00000E+100
         67  traj.phases.descent.collocation_constraint.defects:h                e   0.000000E+00   -6.630524E-12    0.000000E+00              9.00000E+100
         68  traj.phases.descent.collocation_constraint.defects:h                e   0.000000E+00   -7.193601E-12    0.000000E+00              9.00000E+100
         69  traj.phases.descent.collocation_constraint.defects:h                e   0.000000E+00   -7.997998E-12    0.000000E+00              9.00000E+100
         70  traj.phases.descent.collocation_constraint.defects:gam              e   0.000000E+00    1.658619E-14    0.000000E+00              9.00000E+100
         71  traj.phases.descent.collocation_constraint.defects:gam              e   0.000000E+00    5.902797E-15    0.000000E+00              9.00000E+100
         72  traj.phases.descent.collocation_constraint.defects:gam              e   0.000000E+00   -8.349203E-15    0.000000E+00              9.00000E+100
         73  traj.phases.descent.collocation_constraint.defects:gam              e   0.000000E+00   -1.219836E-14    0.000000E+00              9.00000E+100
         74  traj.phases.descent.collocation_constraint.defects:gam              e   0.000000E+00   -1.125571E-14    0.000000E+00              9.00000E+100
         75  traj.phases.descent.collocation_constraint.defects:v                e   0.000000E+00    7.656129E-13    0.000000E+00              9.00000E+100
         76  traj.phases.descent.collocation_constraint.defects:v                e   0.000000E+00    4.427771E-13    0.000000E+00              9.00000E+100
         77  traj.phases.descent.collocation_constraint.defects:v                e   0.000000E+00    1.274825E-12    0.000000E+00              9.00000E+100
         78  traj.phases.descent.collocation_constraint.defects:v                e   0.000000E+00    2.209935E-12    0.000000E+00              9.00000E+100
         79  traj.phases.descent.collocation_constraint.defects:v                e   0.000000E+00    2.736743E-12    0.000000E+00              9.00000E+100
         80  traj.linkages.ascent:time_final|descent:time_initial                e   0.000000E+00    0.000000E+00    0.000000E+00              9.00000E+100
         81  traj.linkages.ascent:r_final|descent:r_initial                      e   0.000000E+00    0.000000E+00    0.000000E+00              9.00000E+100
         82  traj.linkages.ascent:h_final|descent:h_initial                      e   0.000000E+00    0.000000E+00    0.000000E+00              9.00000E+100
         83  traj.linkages.ascent:gam_final|descent:gam_initial                  e   0.000000E+00    1.567954E-20    0.000000E+00              9.00000E+100
         84  traj.linkages.ascent:v_final|descent:v_initial                      e   0.000000E+00    0.000000E+00    0.000000E+00              9.00000E+100
         85  traj.phases.ascent.initial_boundary_constraints.initial_value:ke    i   0.000000E+00    4.000000E+00    4.000000E+00         u    9.00000E+100

--------------------------------------------------------------------------------


Simulating trajectory traj
Done simulating trajectory traj
optimal radius: 0.04185302552964484 m 
cannonball mass: 2.416817832396711 kg 
launch angle: 32.01951038834683 deg
maximum range: 3183.117726801719 m