How to use dymos - 10 common examples

compressed=True):

    if transcription == 'gauss-lobatto':
        t = dm.GaussLobatto(num_segments=num_segments,
                            order=transcription_order,
                            compressed=compressed)
    elif transcription == 'radau-ps':
        t = dm.Radau(num_segments=num_segments,
                     order=transcription_order,
                     compressed=compressed)
    elif transcription == 'runge-kutta':
        t = dm.RungeKutta(num_segments=num_segments,
                          order=transcription_order,
                          compressed=compressed)

    phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)

    return phase

import dymos as dm
        import matplotlib.pyplot as plt
        plt.switch_backend('Agg')  # disable plotting to the screen

        from oscillator_ode import OscillatorODE

        # Instantiate an OpenMDAO Problem instance.
        prob = om.Problem()

        # Instantiate a Dymos Trajectory and add it to the Problem model.
        traj = dm.Trajectory()
        prob.model.add_subsystem('traj', traj)

        # Instantiate a Phase and add it to the Trajectory.
        # Here the transcription is necessary but not particularly relevant.
        phase = dm.Phase(ode_class=OscillatorODE, transcription=dm.Radau(num_segments=4))
        traj.add_phase('phase0', phase)

        # Tell Dymos the states to be propagated using the given ODE.
        phase.add_state('x', rate_source='v', targets=['x'], units='m')
        phase.add_state('v', rate_source='v_dot', targets=['v'], units='m/s')

        # The spring constant, damping coefficient, and mass are inputs to the system that are constant throughout the phase.
        phase.add_input_parameter('k', units='N/m', targets=['k'])
        phase.add_input_parameter('c', units='N*s/m', targets=['c'])
        phase.add_input_parameter('m', units='kg', targets=['m'])

        # Setup the OpenMDAO problem
        prob.setup()

        # Assign values to the times and states
        prob.set_val('traj.phase0.t_initial', 0.0)

plt.switch_backend('Agg')  # disable plotting to the screen

        from oscillator_ode import OscillatorODE

        # Instantiate an OpenMDAO Problem instance.
        prob = om.Problem()

        # We need an optimization driver.  To solve this simple problem ScipyOptimizerDriver will work.
        prob.driver = om.ScipyOptimizeDriver()

        # Instantiate a Dymos Trajectory and add it to the Problem model.
        traj = dm.Trajectory()
        prob.model.add_subsystem('traj', traj)

        # Instantiate a Phase and add it to the Trajectory.
        phase = dm.Phase(ode_class=OscillatorODE, transcription=dm.Radau(num_segments=10))
        traj.add_phase('phase0', phase)

        # Tell Dymos that the duration of the phase is bounded.
        phase.set_time_options(fix_initial=True, fix_duration=True)

        # Tell Dymos the states to be propagated using the given ODE.
        phase.add_state('x', fix_initial=True, rate_source='v', targets=['x'], units='m')
        phase.add_state('v', fix_initial=True, rate_source='v_dot', targets=['v'], units='m/s')

        # The spring constant, damping coefficient, and mass are inputs to the system that are constant throughout the phase.
        phase.add_input_parameter('k', units='N/m', targets=['k'])
        phase.add_input_parameter('c', units='N*s/m', targets=['c'])
        phase.add_input_parameter('m', units='kg', targets=['m'])

        # Since we're using an optimization driver, an objective is required.  We'll minimize the final time in this case.
        phase.add_objective('time', loc='final')

import openmdao.api as om
        import dymos as dm
        import matplotlib.pyplot as plt
        plt.switch_backend('Agg')  # disable plotting to the screen

        from oscillator_ode import OscillatorODE

        # Instantiate an OpenMDAO Problem instance.
        prob = om.Problem()

        # Instantiate a Dymos Trajectory and add it to the Problem model.
        traj = dm.Trajectory()
        prob.model.add_subsystem('traj', traj)

        # Instantiate a Phase and add it to the Trajectory.
        phase = dm.Phase(ode_class=OscillatorODE, transcription=dm.Radau(num_segments=4, solve_segments=True))
        traj.add_phase('phase0', phase)

        # Tell Dymos the states to be propagated using the given ODE.
        phase.add_state('x', fix_initial=True, rate_source='v', targets=['x'], units='m')
        phase.add_state('v', fix_initial=True, rate_source='v_dot', targets=['v'], units='m/s')

        # The spring constant, damping coefficient, and mass are inputs to the system that are constant throughout the phase.
        phase.add_input_parameter('k', units='N/m', targets=['k'])
        phase.add_input_parameter('c', units='N*s/m', targets=['c'])
        phase.add_input_parameter('m', units='kg', targets=['m'])

        # Setup the OpenMDAO problem
        prob.setup()

        # Assign values to the times and states
        prob.set_val('traj.phase0.t_initial', 0.0)

import dymos as dm
        import matplotlib.pyplot as plt
        plt.switch_backend('Agg')  # disable plotting to the screen

        from projectile_ode import ProjectileODE

        # Instnatiate an OpenMDAO Problem instance.
        prob = om.Problem()

        # Instantiate a Dymos Trajectory and add it to the Problem model.
        traj = dm.Trajectory()
        prob.model.add_subsystem('traj', traj)

        # Instantiate a Phase and add it to the Trajectory.
        # Here the transcription is necessary but not particularly relevant.
        phase = dm.Phase(ode_class=ProjectileODE, transcription=dm.Radau(num_segments=10))
        traj.add_phase('phase0', phase)

        # Tell Dymos the states to be propagated using the given ODE.
        phase.add_state('x', rate_source='x_dot', targets=None, units='m')
        phase.add_state('y', rate_source='y_dot', targets=None, units='m')
        phase.add_state('vx', rate_source='vx_dot', targets=['vx'], units='m/s')
        phase.add_state('vy', rate_source='vy_dot', targets=['vy'], units='m/s')

        # Setup the OpenMDAO problem
        prob.setup()

        # Assign values to the times and states
        prob.set_val('traj.phase0.t_initial', 0.0)
        prob.set_val('traj.phase0.t_duration', 15.0)

        prob.set_val('traj.phase0.states:x', 0.0)

from openmdao.utils.assert_utils import assert_near_equal
        import dymos as dm
        from dymos.examples.shuttle_reentry.shuttle_ode import ShuttleODE
        from dymos.examples.plotting import plot_results

        # Instantiate the problem, add the driver, and allow it to use coloring
        p = om.Problem(model=om.Group())
        p.driver = om.pyOptSparseDriver()
        p.driver.declare_coloring()
        p.driver.options['optimizer'] = 'SLSQP'

        # Instantiate the trajectory and add a phase to it
        traj = p.model.add_subsystem('traj', dm.Trajectory())
        phase0 = traj.add_phase('phase0',
                                dm.Phase(ode_class=ShuttleODE,
                                         transcription=dm.Radau(num_segments=15, order=3)))

        phase0.set_time_options(fix_initial=True, units='s', duration_ref=200)
        phase0.add_state('h', fix_initial=True, fix_final=True, units='ft', rate_source='hdot',
                         targets=['h'], lower=0, ref0=75000, ref=300000, defect_ref=1000)
        phase0.add_state('gamma', fix_initial=True, fix_final=True, units='rad',
                         rate_source='gammadot', targets=['gamma'],
                         lower=-89. * np.pi / 180, upper=89. * np.pi / 180)
        phase0.add_state('phi', fix_initial=True, fix_final=False, units='rad',
                         rate_source='phidot', lower=0, upper=89. * np.pi / 180)
        phase0.add_state('psi', fix_initial=True, fix_final=False, units='rad',
                         rate_source='psidot', targets=['psi'], lower=0, upper=90. * np.pi / 180)
        phase0.add_state('theta', fix_initial=True, fix_final=False, units='rad',
                         rate_source='thetadot', targets=['theta'],
                         lower=-89. * np.pi / 180, upper=89. * np.pi / 180)
        phase0.add_state('v', fix_initial=True, fix_final=True, units='ft/s',
                         rate_source='vdot', targets=['v'], lower=0, ref0=2500, ref=25000)

def make_brachistochrone_phase(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
                               compressed=True):

    if transcription == 'gauss-lobatto':
        t = dm.GaussLobatto(num_segments=num_segments,
                            order=transcription_order,
                            compressed=compressed)
    elif transcription == 'radau-ps':
        t = dm.Radau(num_segments=num_segments,
                     order=transcription_order,
                     compressed=compressed)
    elif transcription == 'runge-kutta':
        t = dm.RungeKutta(num_segments=num_segments,
                          order=transcription_order,
                          compressed=compressed)

    phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)

    return phase

p.driver = om.pyOptSparseDriver()
        p.driver.options['optimizer'] = 'SLSQP'
        p.driver.declare_coloring()

        #
        # First Phase: Standard Brachistochrone
        #
        num_segments = 10
        transcription_order = 3
        compressed = False

        tx0 = dm.GaussLobatto(num_segments=num_segments,
                              order=transcription_order,
                              compressed=compressed)

        tx1 = dm.Radau(num_segments=num_segments*2,
                       order=transcription_order*3,
                       compressed=compressed)

        phase0 = dm.Phase(ode_class=BrachistochroneODE, transcription=tx0)

        p.model.add_subsystem('phase0', phase0)

        phase0.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

        phase0.add_state('x', rate_source=BrachistochroneODE.states['x']['rate_source'],
                         units=BrachistochroneODE.states['x']['units'],
                         fix_initial=True, fix_final=False, solve_segments=False)

        phase0.add_state('y', rate_source=BrachistochroneODE.states['y']['rate_source'],
                         units=BrachistochroneODE.states['y']['units'],
                         fix_initial=True, fix_final=False, solve_segments=False)

'units': 'rad'},
                                                           tan_theta={'value': np.ones(nn)}))

                self.add_subsystem('eom', LaunchVehicle2DEOM(num_nodes=nn))

                self.connect('guidance.theta', 'eom.theta')

        #
        # Setup and solve the optimal control problem
        #
        p = om.Problem(model=om.Group())

        traj = p.model.add_subsystem('traj', dm.Trajectory())

        phase = dm.Phase(ode_class=LaunchVehicleLinearTangentODE,
                         transcription=dm.Radau(num_segments=20, order=3, compressed=False))
        traj.add_phase('phase0', phase)

        phase.set_time_options(initial_bounds=(0, 0), duration_bounds=(10, 1000), units='s')

        #
        # Set the state options.  We include rate_source, units, and targets here since the ODE
        # is not decorated with their default values.
        #
        phase.add_state('x', fix_initial=True, lower=0, rate_source='eom.xdot', units='m')
        phase.add_state('y', fix_initial=True, lower=0, rate_source='eom.ydot', units='m')
        phase.add_state('vx', fix_initial=True, lower=0, rate_source='eom.vxdot',
                        units='m/s', targets=['eom.vx'])
        phase.add_state('vy', fix_initial=True, rate_source='eom.vydot',
                        units='m/s', targets=['eom.vy'])
        phase.add_state('m', fix_initial=True, rate_source='eom.mdot',
                        units='kg', targets=['eom.m'])

self.add_subsystem('guidance', om.ExecComp('theta=arctan(tan_theta)',
                                                           theta={'value': np.ones(nn),
                                                                  'units': 'rad'},
                                                           tan_theta={'value': np.ones(nn)}))

                self.add_subsystem('eom', LaunchVehicle2DEOM(num_nodes=nn))

                self.connect('guidance.theta', 'eom.theta')

        #
        # Setup and solve the optimal control problem
        #
        p = om.Problem(model=om.Group())

        traj = p.model.add_subsystem('traj', dm.Trajectory())

        phase = dm.Phase(ode_class=LaunchVehicleLinearTangentODE,
                         transcription=dm.Radau(num_segments=20, order=3, compressed=False))
        traj.add_phase('phase0', phase)

        phase.set_time_options(initial_bounds=(0, 0), duration_bounds=(10, 1000), units='s')

        #
        # Set the state options.  We include rate_source, units, and targets here since the ODE
        # is not decorated with their default values.
        #
        phase.add_state('x', fix_initial=True, lower=0, rate_source='eom.xdot', units='m')
        phase.add_state('y', fix_initial=True, lower=0, rate_source='eom.ydot', units='m')
        phase.add_state('vx', fix_initial=True, lower=0, rate_source='eom.vxdot',
                        units='m/s', targets=['eom.vx'])
        phase.add_state('vy', fix_initial=True, rate_source='eom.vydot',