How to use the dymos.examples.plotting.plot_results function in dymos

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)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta',
                       'range (m)', 'altitude (m)')],
                     title='Single Stage to Orbit Solution Using Linear Tangent Guidance',
                     p_sol=p, p_sim=exp_out)

        plt.show()

dm.run_problem(p)

        # Check the validity of the solution
        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 2008.59,
                          tolerance=1e-3)
        assert_near_equal(p.get_val('traj.phase0.timeseries.states:theta', units='deg')[-1],
                          34.1412, tolerance=1e-3)
        # assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 2181.90371131, tolerance=1e-3)
        # assert_near_equal(p.get_val('traj.phase0.timeseries.states:theta')[-1], .53440626, tolerance=1e-3)

        # Run the simulation to check if the model is physically valid
        sim_out = traj.simulate()

        # Plot the results

        plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:alpha',
                       'time (s)', 'alpha (rad)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:beta',
                       'time (s)', 'beta (rad)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.states:theta',
                       'time (s)', 'theta (rad)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.q',
                       'time (s)', 'q (btu/ft/ft/s')], title='Reentry Solution', p_sol=p,
                     p_sim=sim_out)

        plt.show()

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',
                       'x0 vs x1'),
                     ('traj.phase0.timeseries.time',
                      'traj.phase0.timeseries.controls:u',
                      'time (s)',
                      'control u'),
                      ],

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)

        plt.show()

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,
                          tolerance=1e-2)

        #
        # Simulate the explicit solution and plot the results.
        #
        exp_out = traj.simulate()

        plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.states:x',
                       'time (s)', 'x $(m)$'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:u',
                       'time (s)', 'u $(m/s^2)$')],
                     title='Hyper Sensitive Problem Solution\nRadau Pseudospectral Method',
                     p_sol=p, p_sim=exp_out)

        plt.show()

#
        x = p.get_val('traj.phase0.timeseries.states:x')
        v = p.get_val('traj.phase0.timeseries.states:v')

        assert_near_equal(x[0], 0.0, tolerance=1.0E-4)
        assert_near_equal(x[-1], 0.25, tolerance=1.0E-4)

        assert_near_equal(v[0], 0.0, tolerance=1.0E-4)
        assert_near_equal(v[-1], 0.0, tolerance=1.0E-4)

        #
        # Simulate the explicit solution and plot the results.
        #
        exp_out = traj.simulate()

        plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.states:x',
                       'time (s)', 'x $(m)$'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.states:v',
                       'time (s)', 'v $(m/s)$'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:u',
                       'time (s)', 'u $(m/s^2)$')],
                     title='Double Integrator Solution\nRadau Pseudospectral Method',
                     p_sol=p, p_sim=exp_out)

        plt.show()

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',
                       'x0 vs x1'),
                     ('traj.phase0.timeseries.time',
                      'traj.phase0.timeseries.controls:u',
                      'time (s)',
                      'control u'),
                      ],

def test_vanderpol_for_docs_simulation(self):
        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)

        # Run the problem (simulate only)
        p.run_model()

        # 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',
                       'x0 vs x1'),
                     ('traj.phase0.timeseries.time',
                      'traj.phase0.timeseries.controls:u',
                      'time (s)',
                      'control u'),
                      ],

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\nRK4 Shooting Method',
                     p_sol=p, p_sim=exp_out)

        plt.show()

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

        #
        # Test the results
        #
        assert_near_equal(p.get_val('traj.phase0.t_duration'), 321.0, tolerance=1.0E-1)

        #
        # Get the explicitly simulated solution and plot the results
        #
        exp_out = traj.simulate()

        plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.states:h',
                       'time (s)', 'altitude (m)'),
                      ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:alpha',
                       'time (s)', 'alpha (deg)')],
                     title='Supersonic Minimum Time-to-Climb Solution',
                     p_sol=p, p_sim=exp_out)

        plt.show()