How to use the dymos.Radau function in dymos

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

def double_integrator_direct_collocation(transcription='gauss-lobatto', compressed=True):
    p = Problem(model=Group())
    p.driver = pyOptSparseDriver()
    p.driver.options['dynamic_simul_derivs'] = True

    if transcription == 'gauss-lobatto':
        transcription = GaussLobatto(num_segments=30, order=3, compressed=compressed)
    elif transcription == "radau-ps":
        transcription = Radau(num_segments=30, order=3, compressed=compressed)

    phase = Phase(ode_class=DoubleIntegratorODE, transcription=transcription)
    p.model.add_subsystem('phase0', phase)

    phase.set_time_options(fix_initial=True, fix_duration=True)

    phase.set_state_options('x', fix_initial=True)
    phase.set_state_options('v', fix_initial=True, fix_final=True)

    phase.add_control('u', units='m/s**2', scaler=0.01, continuity=False, rate_continuity=False,
                      rate2_continuity=False, lower=-1.0, upper=1.0)

    # Maximize distance travelled in one second.
    phase.add_objective('x', loc='final', scaler=-1)

    p.model.linear_solver = DirectSolver()

p.driver.options['optimizer'] = optimizer
    p.driver.options['dynamic_simul_derivs'] = True
    if optimizer == 'SNOPT':
        p.driver.opt_settings['Major iterations limit'] = 20
        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
        p.driver.opt_settings["Linesearch tolerance"] = 0.10
        p.driver.opt_settings['iSumm'] = 6
    if optimizer == 'SLSQP':
        p.driver.opt_settings['MAXIT'] = 50

    num_seg = 15
    seg_ends, _ = lgl(num_seg + 1)

    phase = Phase(ode_class=AircraftODE,
                  transcription=Radau(num_segments=num_seg, segment_ends=seg_ends,
                                      order=3, compressed=compressed, solve_segments=solve_segments))

    # Pass Reference Area from an external source
    assumptions = p.model.add_subsystem('assumptions', IndepVarComp())
    assumptions.add_output('S', val=427.8, units='m**2')
    assumptions.add_output('mass_empty', val=1.0, units='kg')
    assumptions.add_output('mass_payload', val=1.0, units='kg')

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

    phase.set_time_options(initial_bounds=(0, 0),
                           duration_bounds=(300, 10000),
                           duration_ref=5600)

    fix_final = True
    if use_boundary_constraints:

p.driver = pyOptSparseDriver()
    p.driver.options['optimizer'] = optimizer
    p.driver.options['dynamic_simul_derivs'] = True

    if optimizer == 'SNOPT':
        p.driver.opt_settings['Major iterations limit'] = 1000
        p.driver.opt_settings['iSumm'] = 6
        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
        p.driver.opt_settings['Function precision'] = 1.0E-12
        p.driver.opt_settings['Linesearch tolerance'] = 0.1
        p.driver.opt_settings['Major step limit'] = 0.5
        # p.driver.opt_settings['Verify level'] = 3

    t = {'gauss-lobatto': GaussLobatto(num_segments=num_seg, order=transcription_order),
         'radau-ps': Radau(num_segments=num_seg, order=transcription_order),
         'runge-kutta': RungeKutta(num_segments=num_seg)}

    phase = Phase(ode_class=MinTimeClimbODE, transcription=t[transcription])

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

    phase.set_time_options(fix_initial=True, duration_bounds=(50, 400),
                           duration_ref=100.0)

    phase.set_state_options('r', fix_initial=True, lower=0, upper=1.0E6,
                            ref=1.0E3, defect_ref=1.0E3, units='m')

    phase.set_state_options('h', fix_initial=True, lower=0, upper=20000.0,
                            ref=1.0E2, defect_ref=1.0E2, units='m')

    phase.set_state_options('v', fix_initial=True, lower=10.0,

import numpy as np

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

import openmdao.api as om
import dymos as dm
from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE

p = om.Problem(model=om.Group())
phase = dm.Phase(ode_class=BrachistochroneODE,
                 transcription=dm.Radau(num_segments=4, order=[3, 5, 3, 5]))
p.model.add_subsystem('phase0', phase)

p.setup()
p['phase0.t_initial'] = 1.0
p['phase0.t_duration'] = 9.0
p.run_model()

grid_data = phase.options['transcription'].grid_data

t_all = p.get_val('phase0.timeseries.time')
t_disc = t_all[grid_data.subset_node_indices['state_disc'], 0]
t_col = t_all[grid_data.subset_node_indices['col'], 0]


def f(x):  # pragma: no cover
    return np.sin(x) / x + 1