How to use the dymos.run_problem function in dymos

p.setup()

        p['traj.phase0.t_initial'] = 0.0
        p['traj.phase0.t_duration'] = 3600.0
        p['traj.phase0.states:range'][:] = phase.interpolate(ys=(0, 724.0), nodes='state_input')
        p['traj.phase0.states:mass_fuel'][:] = phase.interpolate(ys=(30000, 1e-3), nodes='state_input')
        p['traj.phase0.states:alt'][:] = 10.0

        p['traj.phase0.controls:mach'][:] = 0.8

        p['assumptions.S'] = 427.8
        p['assumptions.mass_empty'] = 0.15E6
        p['assumptions.mass_payload'] = 84.02869 * 400

        dm.run_problem(p)

        assert_near_equal(p.get_val('traj.phase0.timeseries.states:range', units='NM')[-1],
                          726.85, tolerance=1.0E-2)

        exp_out = traj.simulate()

        plot_results([('traj.phase0.timeseries.states:range', 'traj.phase0.timeseries.states:alt',
                       'range (NM)', 'altitude (kft)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.states:mass_fuel',
                       'time (s)', 'fuel mass (lbm)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.CL',
                       'time (s)', 'lift coefficient'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.CD',
                       'time (s)', 'drag coefficient')],
                     title='Commercial Aircraft Optimization',
                     p_sol=p, p_sim=exp_out)

#
        # Set the initial values
        #
        p['traj.phase0.t_initial'] = 0.0
        p['traj.phase0.t_duration'] = 2.0

        p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input')
        p['traj.phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input')
        p['traj.phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input')
        p['traj.phase0.controls:theta'] = phase.interpolate(ys=[5, 100.5], nodes='control_input')

        #
        # Solve for the optimal trajectory
        #
        dm.run_problem(p)

        # Test the results
        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016,
                          tolerance=1.0E-3)

        # Generate the explicitly simulated trajectory
        exp_out = traj.simulate()

        plot_results([('traj.phase0.timeseries.states:x', 'traj.phase0.timeseries.states:y',
                       'x (m)', 'y (m)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:theta',
                       'time (s)', 'theta (deg)')],
                     title='Brachistochrone Solution\nRadau Pseudospectral Method',
                     p_sol=p, p_sim=exp_out)

        plt.show()

phase.add_boundary_constraint('x', loc='final', equals=1)

        phase.add_objective('xL', loc='final')

        p.setup(check=True)

        p.set_val('traj.phase0.states:x', phase.interpolate(ys=[1.5, 1], nodes='state_input'))
        p.set_val('traj.phase0.states:xL', phase.interpolate(ys=[0, 1], nodes='state_input'))
        p.set_val('traj.phase0.t_initial', 0)
        p.set_val('traj.phase0.t_duration', tf)
        p.set_val('traj.phase0.controls:u', phase.interpolate(ys=[-0.6, 2.4], nodes='control_input'))

        #
        # Solve the problem.
        #
        dm.run_problem(p)

        #
        # Verify that the results are correct.
        #
        ui, uf, J = solution()

        assert_near_equal(p.get_val('traj.phase0.timeseries.controls:u')[0],
                          ui,
                          tolerance=1.5e-2)

        assert_near_equal(p.get_val('traj.phase0.timeseries.controls:u')[-1],
                          uf,
                          tolerance=1.5e-2)

        assert_near_equal(p.get_val('traj.phase0.timeseries.states:xL')[-1],
                          J,

def test_vanderpol_for_docs_optimize(self):
        import dymos as dm
        from dymos.examples.plotting import plot_results
        from dymos.examples.vanderpol.vanderpol_dymos import vanderpol

        # Create the Dymos problem instance
        p = vanderpol(transcription='gauss-lobatto', num_segments=75,
                      transcription_order=3, compressed=True, optimizer='SLSQP')

        # Find optimal control solution to stop oscillation
        dm.run_problem(p)

        # check validity by using scipy.integrate.solve_ivp to integrate the solution
        exp_out = p.model.traj.simulate()

        # Display the results
        plot_results([('traj.phase0.timeseries.time',
                       'traj.phase0.timeseries.states:x1',
                       'time (s)',
                       'x1 (V)'),
                     ('traj.phase0.timeseries.time',
                      'traj.phase0.timeseries.states:x0',
                      'time (s)',
                      'x0 (V/s)'),
                      ('traj.phase0.timeseries.states:x0',
                       'traj.phase0.timeseries.states:x1',
                       'x0 vs x1',

p.set_val('traj.burn2.states:theta', value=burn2.interpolate(ys=[0, 4.0],
                                                                     nodes='state_input'))
        p.set_val('traj.burn2.states:vr', value=burn2.interpolate(ys=[0, 0],
                                                                  nodes='state_input'))
        p.set_val('traj.burn2.states:vt',
                  value=burn2.interpolate(ys=[1, np.sqrt(1 / 3.)],
                                          nodes='state_input'))
        p.set_val('traj.burn2.states:deltav',
                  value=burn2.interpolate(ys=[0.1, 0.2], nodes='state_input'))
        p.set_val('traj.burn2.states:accel', value=burn2.interpolate(ys=[0.1, 0],
                                                                     nodes='state_input'))

        p.set_val('traj.burn2.controls:u1', value=burn2.interpolate(ys=[0, 0],
                                                                    nodes='control_input'))

        dm.run_problem(p)

        assert_near_equal(p.get_val('traj.burn2.states:deltav')[-1], 0.3995,
                          tolerance=2.0E-3)

        #
        # Plot results
        #
        traj = p.model.traj
        exp_out = traj.simulate()

        fig = plt.figure(figsize=(8, 4))
        fig.suptitle('Two Burn Orbit Raise Solution')
        ax_u1 = plt.subplot2grid((2, 2), (0, 0))
        ax_deltav = plt.subplot2grid((2, 2), (1, 0))
        ax_xy = plt.subplot2grid((2, 2), (0, 1), rowspan=2)

def test_vanderpol_for_docs_optimize_refine(self):
        import dymos as dm
        from dymos.examples.plotting import plot_results
        from dymos.examples.vanderpol.vanderpol_dymos import vanderpol

        # Create the Dymos problem instance
        p = vanderpol(transcription='gauss-lobatto', num_segments=15,
                      transcription_order=3, compressed=True, optimizer='SLSQP')

        # Enable grid refinement and find optimal control solution to stop oscillation
        p.model.traj.phases.phase0.set_refine_options(refine=True)
        dm.run_problem(p, refine=True, refine_iteration_limit=10)

        # check validity by using scipy.integrate.solve_ivp to integrate the solution
        exp_out = p.model.traj.simulate()

        # Display the results
        plot_results([('traj.phase0.timeseries.time',
                       'traj.phase0.timeseries.states:x1',
                       'time (s)',
                       'x1 (V)'),
                     ('traj.phase0.timeseries.time',
                      'traj.phase0.timeseries.states:x0',
                      'time (s)',
                      'x0 (V/s)'),
                      ('traj.phase0.timeseries.states:x0',
                       'traj.phase0.timeseries.states:x1',
                       'x0 vs x1',

#
        p.setup(check=True)

        p['traj.phase0.t_initial'] = 0.0
        p['traj.phase0.t_duration'] = 150.0
        p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 1.15E5], nodes='state_input')
        p['traj.phase0.states:y'] = phase.interpolate(ys=[0, 1.85E5], nodes='state_input')
        p['traj.phase0.states:vx'] = phase.interpolate(ys=[0, 7796.6961], nodes='state_input')
        p['traj.phase0.states:vy'] = phase.interpolate(ys=[1.0E-6, 0], nodes='state_input')
        p['traj.phase0.states:m'] = phase.interpolate(ys=[117000, 1163], nodes='state_input')
        p['traj.phase0.controls:theta'] = phase.interpolate(ys=[1.5, -0.76], nodes='control_input')

        #
        # Solve the Problem
        #
        dm.run_problem(p)

        #
        # Check the results.
        #
        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 143, tolerance=0.05)
        assert_near_equal(p.get_val('traj.phase0.timeseries.states:y')[-1], 1.85E5, 1e-4)
        assert_near_equal(p.get_val('traj.phase0.timeseries.states:vx')[-1], 7796.6961, 1e-4)
        assert_near_equal(p.get_val('traj.phase0.timeseries.states:vy')[-1], 0, 1e-4)
        #
        # Get the explitly simulated results
        #
        exp_out = traj.simulate()

        #
        # Plot the results
        #

phase.add_timeseries_output('guidance.theta', units='deg')

        p.setup(check=True)

        p['traj.phase0.t_initial'] = 0.0
        p['traj.phase0.t_duration'] = 500.0
        p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 350000.0], nodes='state_input')
        p['traj.phase0.states:y'] = phase.interpolate(ys=[0, 185000.0], nodes='state_input')
        p['traj.phase0.states:vx'] = phase.interpolate(ys=[0, 1627.0], nodes='state_input')
        p['traj.phase0.states:vy'] = phase.interpolate(ys=[1.0E-6, 0], nodes='state_input')
        p['traj.phase0.states:m'] = phase.interpolate(ys=[50000, 50000], nodes='state_input')
        p['traj.phase0.design_parameters:a_ctrl'] = -0.01
        p['traj.phase0.design_parameters:b_ctrl'] = 3.0

        dm.run_problem(p)

        #
        # Check the results.
        #
        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 481, tolerance=0.01)

        #
        # Get the explitly simulated results
        #
        exp_out = traj.simulate()

        #
        # Plot the results
        #
        plot_results([('traj.phase0.timeseries.states:x', 'traj.phase0.timeseries.states:y',
                       'range (m)', 'altitude (m)'),

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

        dm.run_problem(p)

        assert_near_equal(p.get_val('traj.descent.states:r')[-1],
                          3183.25, tolerance=1.0E-2)

        exp_out = traj.simulate()

        print('optimal radius: {0:6.4f} m '.format(p.get_val('external_params.radius',
                                                             units='m')[0]))
        print('cannonball mass: {0:6.4f} kg '.format(p.get_val('size_comp.mass',
                                                               units='kg')[0]))
        print('launch angle: {0:6.4f} '
              'deg '.format(p.get_val('traj.ascent.timeseries.states:gam',  units='deg')[0, 0]))
        print('maximum range: {0:6.4f} '
              'm '.format(p.get_val('traj.descent.timeseries.states:r')[-1, 0]))

        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))

p.model.linear_solver = om.DirectSolver()
        p.setup(check=True)

        p['phase0.t_initial'] = 0.0
        p['phase0.t_duration'] = 2.0

        p['phase0.states:x'] = phase0.interpolate(ys=[0, 10], nodes='state_input')
        p['phase0.states:y'] = phase0.interpolate(ys=[10, 5], nodes='state_input')
        p['phase0.states:v'] = phase0.interpolate(ys=[0, 9.9], nodes='state_input')
        p['phase0.controls:theta'] = phase0.interpolate(ys=[5, 100], nodes='control_input')
        p['phase0.input_parameters:g'] = 9.80665

        p['phase1.states:S'] = 0.0

        dm.run_problem(p)

        expected = np.sqrt((10-0)**2 + (10 - 5)**2)
        assert_near_equal(p['phase1.timeseries.states:S'][-1], expected, tolerance=1.0E-3)

        fig, (ax0, ax1) = plt.subplots(2, 1)
        fig.tight_layout()
        ax0.plot(p.get_val('phase0.timeseries.states:x'), p.get_val('phase0.timeseries.states:y'), '.')
        ax0.set_xlabel('x (m)')
        ax0.set_ylabel('y (m)')
        ax1.plot(p.get_val('phase1.timeseries.time'), p.get_val('phase1.timeseries.states:S'), '+')
        ax1.set_xlabel('t (s)')
        ax1.set_ylabel('S (m)')
        plt.show()