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.227715 sec.
Time to compute coloring: 0.073810 sec.
Memory to compute coloring: 0.375000 MB.
Full total jacobian was computed 3 times, taking 0.495898 seconds.
Total jacobian shape: (98, 100) 


Jacobian shape: (98, 100)  (18.77% 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.495898 sec.
Time to compute coloring: 0.144609 sec.
Memory to compute coloring: 2.929688 MB.
Optimization Problem -- Optimization using pyOpt_sparse
================================================================================
    Objective Function: _objfunc

    Solution: 
--------------------------------------------------------------------------------
    Total Time:                   30.5223
       User Objective Time :      26.3832
       User Sensitivity Time :     3.5949
       Interface Time :            0.3116
       Opt Solver Time:            0.2326
    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.433787E-12    0.000000E+00              9.00000E+100
          1  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00   -3.751666E-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   -4.320100E-12    0.000000E+00              9.00000E+100
          4  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00   -3.410605E-12    0.000000E+00              9.00000E+100
          5  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00   -2.273737E-12    0.000000E+00              9.00000E+100
          6  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00   -2.273737E-12    0.000000E+00              9.00000E+100
          7  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00   -6.821210E-13    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    0.000000E+00    0.000000E+00              9.00000E+100
         10  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00   -1.364242E-12    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    0.000000E+00    0.000000E+00              9.00000E+100
         13  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00    1.364242E-12    0.000000E+00              9.00000E+100
         14  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00    9.094947E-13    0.000000E+00              9.00000E+100
         15  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00    1.364242E-12    0.000000E+00              9.00000E+100
         16  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00    2.273737E-12    0.000000E+00              9.00000E+100
         17  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00    1.818989E-12    0.000000E+00              9.00000E+100
         18  traj.phases.phase.continuity_comp.defect_states:r      e   0.000000E+00    1.364242E-12    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    9.094947E-13    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.364242E-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.194156E-11    0.000000E+00              9.00000E+100
         25  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00   -1.307399E-11    0.000000E+00              9.00000E+100
         26  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00   -7.162271E-12    0.000000E+00              9.00000E+100
         27  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00   -4.206413E-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   -2.273737E-13    0.000000E+00              9.00000E+100
         30  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00   -2.273737E-13    0.000000E+00              9.00000E+100
         31  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    9.094947E-13    0.000000E+00              9.00000E+100
         32  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    3.410605E-13    0.000000E+00              9.00000E+100
         33  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    4.547474E-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    0.000000E+00    0.000000E+00              9.00000E+100
         36  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    4.547474E-13    0.000000E+00              9.00000E+100
         37  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    2.273737E-13    0.000000E+00              9.00000E+100
         38  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    5.684342E-13    0.000000E+00              9.00000E+100
         39  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    1.023182E-12    0.000000E+00              9.00000E+100
         40  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    1.818989E-12    0.000000E+00              9.00000E+100
         41  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    1.818989E-12    0.000000E+00              9.00000E+100
         42  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    2.387424E-12    0.000000E+00              9.00000E+100
         43  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    2.899014E-12    0.000000E+00              9.00000E+100
         44  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    3.126388E-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.467449E-12    0.000000E+00              9.00000E+100
         47  traj.phases.phase.continuity_comp.defect_states:h      e   0.000000E+00    3.524292E-12    0.000000E+00              9.00000E+100
         48  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.887379E-15    0.000000E+00              9.00000E+100
         49  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    7.771561E-16    0.000000E+00              9.00000E+100
         50  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.054712E-15    0.000000E+00              9.00000E+100
         51  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.887379E-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    5.051515E-15    0.000000E+00              9.00000E+100
         54  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    6.883383E-15    0.000000E+00              9.00000E+100
         55  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    9.714451E-15    0.000000E+00              9.00000E+100
         56  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.149081E-14    0.000000E+00              9.00000E+100
         57  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.311885E-14    0.000000E+00              9.00000E+100
         58  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.301736E-14    0.000000E+00              9.00000E+100
         59  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    1.132427E-14    0.000000E+00              9.00000E+100
         60  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    8.493206E-15    0.000000E+00              9.00000E+100
         61  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00    4.829470E-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   -8.881784E-16    0.000000E+00              9.00000E+100
         64  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00   -2.220446E-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   -3.219647E-15    0.000000E+00              9.00000E+100
         67  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00   -2.997602E-15    0.000000E+00              9.00000E+100
         68  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00   -2.664535E-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   -2.220446E-15    0.000000E+00              9.00000E+100
         71  traj.phases.phase.continuity_comp.defect_states:gam    e   0.000000E+00   -2.220446E-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.242296E-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.705303E-12    0.000000E+00              9.00000E+100
         76  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00    1.023182E-12    0.000000E+00              9.00000E+100
         77  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00    6.252776E-13    0.000000E+00              9.00000E+100
         78  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00    3.694822E-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    2.842171E-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   -1.989520E-13    0.000000E+00              9.00000E+100
         83  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00   -3.268497E-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.126388E-13    0.000000E+00              9.00000E+100
         86  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00   -2.415845E-13    0.000000E+00              9.00000E+100
         87  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00   -2.131628E-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    7.105427E-14    0.000000E+00              9.00000E+100
         90  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00    5.684342E-14    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.700062E-13    0.000000E+00              9.00000E+100
         93  traj.phases.phase.continuity_comp.defect_states:v      e   0.000000E+00    3.410605E-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    3.979039E-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.04186845473893394 m 
cannonball mass: 2.4194917134068668 kg 
launch angle: 32.01931307430318 deg
maximum range: 3182.962437214954 m
../../_images/55edd217d3bb77531a38783ae8d8bce1e5e41b652fe744be7ae7db496fec0b24.png ../../_images/bedce7f7fc5075383f285f2e23c18c6888eca4ae8e87dd1f0b63996821437ea6.png