Phase-specific ode options (#155)

* ODEOptions are now optional on an ODE system.  All ODE options can be set via the phase set_time_options, set_state_options, add_control, add_design_parameter, and add_input_parameter interfaces.
So far only brachistochrone has been tested.  More tests that include input parameters and traj parameters need to be added.

* trajectories now starting to work under the new ODE options procedures.  traj_parameters removed from phases -> trajectories now add the appropriate input parameters to each phase when setting up their input and design parameters.

* Phase option 'ode_class' can now be set in phase.setup methods if inheriting from a Phase class.
Phase time, state, control, and parameter options can now be set in phase.setup if inheriting from a Phase class.
If doing this then the call to super.setup should occur after these have been set.
Temporarily removed rate_param operability from simulate since rate_params have been changed to rate_targets.
Proper rate_target handling will occur in SimulationPhase refactor.

* fixes for polynomial controls

* add_objective is now a method on PhaseBase which caches objective information until setup.
During setup, the _setup_objective method determines the path to the objective output and adds it through the standard openMDAO methods.

* Moved time_options checking to the Phase setup stack.
Removed deprecated time options opt_initial and opt_duration.

* working out issues with steady flight example

* Non-opt controls automatically have continuity/rate continuity disabled.
Only trajectory linkages remain to be moved over to the new system in which time/state/control/parameter options are not available until setup.

* two phase cannonball working with constrained linkages

* Trajectory phase linkages working for both connected and constrained linkages.

* successful test of segment simulation component for use in SolveIVPPhase

* SolveIVPPhase successful in integrating the simple ODE forward and backward at time.
SolveIVPPhase will no longer support the times argument.  All outputs are provided
at all nodes of the given grid.

* SolveIVPPhase now supports dense output with the output_nodes_per_seg option.
If None, outputs are provided at 'all' nodes as defined in the GridData.
If given as an integer (n), then each segment provides outputs at n equally
distributed time points in the segment.

* New SolveIVPPhase simulation is working from Trajectory.simulate().

* Tweaks to readme to reflect changes to dymos.
Trajectory linkages now issue error for invalid linkage variable names.
Linkage report now displays constraint linkages with `=` symbols and
connection linkages with `-->` symbols.

* Phase linkages now raise ValueError for invalid phase names

.github/PULL_REQUEST_TEMPLATE.md

@@ -8,7 +8,11 @@ Summary of PR.
 
 ### Status
 
-- [x] Ready for merge
+- [ ] Ready for merge
+
+### Backwards incompatibilities
+
+None
 
 ### New Dependencies
 

dymos/examples/aircraft_steady_flight/ex_aircraft_steady_flight.py

@@ -7,8 +7,8 @@
 
 import numpy as np
 
-from openmdao.api import Problem, Group, pyOptSparseDriver
-from openmdao.api import IndepVarComp
+from openmdao.api import Problem, Group, pyOptSparseDriver, IndepVarComp
+from openmdao.utils.general_utils import set_pyoptsparse_opt
 
 from dymos import Phase
 
@@ -21,6 +21,7 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto',
                               use_boundary_constraints=False, compressed=False):
     p = Problem(model=Group())
     p.driver = pyOptSparseDriver()
+    _, optimizer = set_pyoptsparse_opt(optimizer, fallback=False)
     p.driver.options['optimizer'] = optimizer
     p.driver.options['dynamic_simul_derivs'] = True
     if optimizer == 'SNOPT':
@@ -71,7 +72,7 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto',
                             solve_segments=solve_segments)
 
     phase.add_control('climb_rate', units='ft/min', opt=True, lower=-3000, upper=3000,
-                      rate_continuity=True)
+                      rate_continuity=True, rate2_continuity=False)
 
     phase.add_control('mach', units=None, opt=False)
 
@@ -104,8 +105,7 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto',
     p.run_driver()
 
     if show_plots:
-        exp_out = phase.simulate(times=np.linspace(0, p['phase0.t_duration'], 500), record=True,
-                                 record_file='test_ex_aircraft_steady_flight_rec.db')
+        exp_out = phase.simulate()
 
         t_imp = p.get_val('phase0.timeseries.time')
         t_exp = exp_out.get_val('phase0.timeseries.time')
@@ -119,7 +119,8 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto',
         mass_fuel_imp = p.get_val('phase0.timeseries.states:mass_fuel', units='kg')
         mass_fuel_exp = exp_out.get_val('phase0.timeseries.states:mass_fuel', units='kg')
 
-        plt.plot(t_imp, alt_imp, 'ro')
+        plt.show()
+        plt.plot(t_imp, alt_imp, 'b-')
         plt.plot(t_exp, alt_exp, 'b-')
         plt.suptitle('altitude vs time')
 

dymos/examples/aircraft_steady_flight/test/test_doc_aircraft_steady_flight.py

@@ -96,8 +96,7 @@ def test_steady_aircraft_for_docs(self):
 
         p.run_driver()
 
-        exp_out = phase.simulate(times=np.linspace(0, p['phase0.t_duration'], 500), record=True,
-                                 record_file='test_doc_aircraft_steady_flight_rec.db')
+        exp_out = phase.simulate()
 
         time_imp = p.get_val('phase0.timeseries.time', units='s')
         time_exp = exp_out.get_val('phase0.timeseries.time', units='s')

dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py

@@ -29,5 +29,6 @@ def test_ex_aircraft_steady_flight_solve(self):
         assert_rel_error(self, p.get_val('phase0.timeseries.states:range', units='NM')[-1],
                          726.85, tolerance=1.0E-2)
 
+
 if __name__ == '__main__':
     unittest.main()

dymos/examples/brachistochrone/ex_brachistochrone.py

@@ -68,7 +68,7 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 
     # Plot results
     if SHOW_PLOTS:
-        exp_out = phase.simulate(times=50, record=False)
+        exp_out = phase.simulate()
 
         fig, ax = plt.subplots()
         fig.suptitle('Brachistochrone Solution')

dymos/examples/brachistochrone/ex_brachistochrone_vector_states.py

@@ -79,7 +79,7 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     # Plot results
     if SHOW_PLOTS:
         p.run_driver()
-        exp_out = phase.simulate(times=50, record_file=sim_record)
+        exp_out = phase.simulate(record_file=sim_record)
 
         fig, ax = plt.subplots()
         fig.suptitle('Brachistochrone Solution')

dymos/examples/brachistochrone/test/test_brachistochrone_integrated_control.py

@@ -4,6 +4,7 @@
 import unittest
 
 import numpy as np
+from scipy.interpolate import interp1d
 
 from openmdao.api import ExplicitComponent
 from dymos import declare_time, declare_state, declare_parameter
@@ -144,10 +145,33 @@ def test_brachistochrone_integrated_control_gauss_lobatto(self):
         # Test the results
         assert_rel_error(self, p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
 
-        # Generate the explicitly simulated trajectory
-        t0 = p['phase0.t_initial']
-        tf = t0 + p['phase0.t_duration']
-        phase.simulate(times=np.linspace(t0, tf, 50))
+        sim_out = phase.simulate(times_per_seg=20)
+
+        x_sol = p.get_val('phase0.timeseries.states:x')
+        y_sol = p.get_val('phase0.timeseries.states:y')
+        v_sol = p.get_val('phase0.timeseries.states:v')
+        theta_sol = p.get_val('phase0.timeseries.states:theta')
+        theta_dot_sol = p.get_val('phase0.timeseries.controls:theta_dot')
+        time_sol = p.get_val('phase0.timeseries.time')
+
+        x_sim = sim_out.get_val('phase0.timeseries.states:x')
+        y_sim = sim_out.get_val('phase0.timeseries.states:y')
+        v_sim = sim_out.get_val('phase0.timeseries.states:v')
+        theta_sim = sim_out.get_val('phase0.timeseries.states:theta')
+        theta_dot_sim = sim_out.get_val('phase0.timeseries.controls:theta_dot')
+        time_sim = sim_out.get_val('phase0.timeseries.time')
+
+        x_interp = interp1d(time_sim[:, 0], x_sim[:, 0])
+        y_interp = interp1d(time_sim[:, 0], y_sim[:, 0])
+        v_interp = interp1d(time_sim[:, 0], v_sim[:, 0])
+        theta_interp = interp1d(time_sim[:, 0], theta_sim[:, 0])
+        theta_dot_interp = interp1d(time_sim[:, 0], theta_dot_sim[:, 0])
+
+        assert_rel_error(self, x_interp(time_sol), x_sol, tolerance=1.0E-5)
+        assert_rel_error(self, y_interp(time_sol), y_sol, tolerance=1.0E-5)
+        assert_rel_error(self, v_interp(time_sol), v_sol, tolerance=1.0E-5)
+        assert_rel_error(self, theta_interp(time_sol), theta_sol, tolerance=1.0E-5)
+        assert_rel_error(self, theta_dot_interp(time_sol), theta_dot_sol, tolerance=1.0E-5)
 
     def test_brachistochrone_integrated_control_radau_ps(self):
         import numpy as np
@@ -197,10 +221,32 @@ def test_brachistochrone_integrated_control_radau_ps(self):
         p.run_driver()
 
         # Test the results
-        tf = p.get_val('phase0.time')[-1]
-        assert_rel_error(self, tf, 1.8016, tolerance=1.0E-3)
+        assert_rel_error(self, p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
 
-        # Generate the explicitly simulated trajectory
-        t0 = p['phase0.t_initial']
-        tf = t0 + p['phase0.t_duration']
-        phase.simulate(times=np.linspace(t0, tf, 50))
+        sim_out = phase.simulate(times_per_seg=20)
+
+        x_sol = p.get_val('phase0.timeseries.states:x')
+        y_sol = p.get_val('phase0.timeseries.states:y')
+        v_sol = p.get_val('phase0.timeseries.states:v')
+        theta_sol = p.get_val('phase0.timeseries.states:theta')
+        theta_dot_sol = p.get_val('phase0.timeseries.controls:theta_dot')
+        time_sol = p.get_val('phase0.timeseries.time')
+
+        x_sim = sim_out.get_val('phase0.timeseries.states:x')
+        y_sim = sim_out.get_val('phase0.timeseries.states:y')
+        v_sim = sim_out.get_val('phase0.timeseries.states:v')
+        theta_sim = sim_out.get_val('phase0.timeseries.states:theta')
+        theta_dot_sim = sim_out.get_val('phase0.timeseries.controls:theta_dot')
+        time_sim = sim_out.get_val('phase0.timeseries.time')
+
+        x_interp = interp1d(time_sim[:, 0], x_sim[:, 0])
+        y_interp = interp1d(time_sim[:, 0], y_sim[:, 0])
+        v_interp = interp1d(time_sim[:, 0], v_sim[:, 0])
+        theta_interp = interp1d(time_sim[:, 0], theta_sim[:, 0])
+        theta_dot_interp = interp1d(time_sim[:, 0], theta_dot_sim[:, 0])
+
+        assert_rel_error(self, x_interp(time_sol), x_sol, tolerance=1.0E-5)
+        assert_rel_error(self, y_interp(time_sol), y_sol, tolerance=1.0E-5)
+        assert_rel_error(self, v_interp(time_sol), v_sol, tolerance=1.0E-5)
+        assert_rel_error(self, theta_interp(time_sol), theta_sol, tolerance=1.0E-5)
+        assert_rel_error(self, theta_dot_interp(time_sol), theta_dot_sol, tolerance=1.0E-5)