Vectorize all the things (OpenMDAO#103)

* added ODE and a notional spacecraft reorientation test problem to test out vectorized states, controls, and parameters.

* fixed indexing for RPM phases, PhaseSimulationResults remains to be fixed.  uncovered an issue with computing state value continuity when using "compressed" transcription.  ex_min_time_climb.py demonstrates the issue.

* added a brachistochrone in which the position is a vector, uncovered a bug in simulated results

* fixed recording for parameters, now design and input parameters dont have their entire time-history saved to the file.  Updated the commercial aircraft example to integrate climb_rate to get altitude, and not surprisingly its much more consistent in its result now.

* error in sparsity of test_aircraft_cruise when using gauss-lobatto, debugging

* fixed test by making climb_rate a control and not a design parameter

* removed spacecraft reorientation example for now

* restored use of the Agg backend in examples

* fix to install openmdao documentation

* retry travis fix

* fixed an error in table formatting

* fixed an error in table formatting

* Added unit tests for indexing, finite burn orbit EOM, and testing of the vector states brachistochrone example.

* utils.misc unit tests

* updated double integrator example to use other states/controls as rate sources.  parameters no longer need to have targets, for degenerate cases such as this where the parameter is purely there as a rate source.

* pep8 fixes

* disc -> state_disc in gauss lobatto get rate source path

* disc -> state_disc

* sizing issue with design and input parameters in phase simulation results

* double integrator example switched back to SLSQP

Author: robfalck

.travis.yml

@@ -113,7 +113,7 @@ install:
 # install OpenMDAO in developer mode so we have access to its sphinx extensions
 - git clone https://github.com/OpenMDAO/OpenMDAO.git;
 - cd OpenMDAO;
-- pip install -e .;
+- pip install -e .[all];
 - cd ..;
 
 # install dymos itself in developer mode.

dymos/docs/examples/commercial_aircraft.rst

@@ -21,7 +21,7 @@ In this case we also use the steady flight assumption to allow the problem
 to be parameterized with a different control set.  Rather than providing
 throttle/thrust and angle of attack as control variables, we use a *differential
 inclusion* approach.  Whereas altitude and velocity are often treated as
-state variables to be integrated, here we specify them as controls.
+state variables to be integrated, here we take a slightly different approach.
 
 Rather than using the differential equations to compute rates for the
 flight path angle and true airspeed, we use a nonlinear solver to determine
@@ -41,6 +41,7 @@ presence of nonlinear solvers can lead to inaccurate derivatives based on the pr
 of the solver.  With OpenMDAO's unified derivatives framework, the sensitivy to
 solver tolerances is greatly reduced.
 
+
 Pitfalls of Differential Inclusion
 ==================================
 
@@ -51,6 +52,17 @@ of the vehicle.  In our case, the throttle parameter :math:`$\tau$` is path cons
 within the optimizer, so that the optimizer can make the solution physically attainable
 by our model.
 
+We have two options for how to tackle this problem.  We could specify altitude as a control
+variable, and use it's approximate derivative (`alt_rate`) to provide the altitude rate for
+the flight equilibrium conditions.  In practice, however, this noisy derivative of altitude
+can cause numerical difficulties.  Other times, the optimizer will shape the control to satisfy
+our equilibrium conditions at the collocation nodes but be wildly wrong in between these nodes.
+One simple way to get around this is to still integrate altitude, but to provide `climb_rate`
+as the control variable.  This makes it easy to constraint `climb_rate` with bounds constraints,
+and the time-history of altitude will be smoother since it is the time integral of the rate of
+climb.
+
+
 Solver Options
 ##############
 
@@ -105,16 +117,17 @@ Name       Description                     Fixed Initial Value Fixed Final Value
 =========  ==============================  =================== ===================
 range      distance flown along ground     True (0 NM)         False
 mass_fuel  mass of fuel aboard aircraft    True (30000 lbm)    True (0 lbm)
+alt        aircraft altitude               True (10000 ft)     True (10000 ft)
 =========  ==============================  =================== ===================
 
 Dynamic Controls
 ================
 
-=====  ==============================  =========  =================== ===================
-Name   Description                     Optimized  Fixed Initial Value Fixed Final Value
-=====  ==============================  =========  =================== ===================
-alt    aircraft altitude               True       True (10000 ft)     True (10000 ft)
-=====  ==============================  =========  =================== ===================
+==========       ==============================  =========  =================== ===================
+Name             Description                     Optimized  Fixed Initial Value Fixed Final Value
+==========       ==============================  =========  =================== ===================
+climb_rate       aircraft rate of climb          True       False               False
+==========       ==============================  =========  =================== ===================
 
 Design Parameters
 =================
@@ -144,7 +157,6 @@ Nonlinear Path Constraints
 Name          Location (initial or final)     Lower         Upper
 ============  ==============================  ============  ==============
 tau           engine throttle parameter       0.01          1.0
-alt_rate      climb rate                      -3000 ft/min  3000 ft/min
 ============  ==============================  ============  ==============
 
 Nonlinear Boundary Constraints

dymos/docs/examples/double_integrator.rst

@@ -30,6 +30,11 @@ in the control.
 1. The ODE System: double_integrator_ode.py
 -------------------------------------------
 
+This problem is unique in that we don't actually have to calculate anything in the ODE.  For the
+sake of Dymos, we create an ExplicitComponent and provide it with the `num_nodes` option, but
+it has no inputs and no outputs.  The rates for the states are entirely provided by the states
+and controls.
+
 .. embed-code::
     dymos.examples.double_integrator.double_integrator_ode
     :layout: code

dymos/examples/aircraft_steady_flight/aircraft_ode.py

@@ -18,7 +18,9 @@
 @declare_state('range', rate_source='range_rate_comp.dXdt:range', units='m')
 @declare_state('mass_fuel', targets=['mass_comp.mass_fuel'],
                rate_source='propulsion.dXdt:mass_fuel', units='kg')
-@declare_parameter('alt', targets=['atmos.h', 'aero.alt', 'propulsion.alt'], units='m')
+@declare_state('alt', targets=['atmos.h', 'aero.alt', 'propulsion.alt'],
+               rate_source='climb_rate', units='m')
+# @declare_parameter('alt', targets=['atmos.h', 'aero.alt', 'propulsion.alt'], units='m')
 @declare_parameter('climb_rate', targets=['gam_comp.climb_rate'], units='m/s')
 @declare_parameter('mach', targets=['tas_comp.mach', 'aero.mach'], units='m/s')
 @declare_parameter('S', targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'], units='m**2')

dymos/examples/aircraft_steady_flight/ex_aircraft_steady_flight.py

@@ -27,7 +27,7 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto'):
         p.driver.opt_settings["Linesearch tolerance"] = 0.10
         p.driver.opt_settings['iSumm'] = 6
 
-    num_seg = 15
+    num_seg = 20
     seg_ends, _ = lgl(num_seg + 1)
 
     phase = Phase(transcription,
@@ -53,12 +53,11 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto'):
                             defect_scaler=1.0E-2)
     phase.set_state_options('mass_fuel', units='lbm', fix_initial=True, fix_final=True,
                             upper=1.5E5, lower=0.0, scaler=1.0E-5, defect_scaler=1.0E-1)
+    phase.set_state_options('alt', units='kft', fix_initial=True, fix_final=True, lower=0.0,
+                            upper=60)
 
-    phase.add_control('alt', units='kft', opt=True, lower=0.0, upper=50.0,
-                      rate_param='climb_rate',
-                      rate_continuity=True, rate_continuity_scaler=1.0,
-                      rate2_continuity=True, rate2_continuity_scaler=1.0, ref=1.0,
-                      fix_initial=True, fix_final=True)
+    phase.add_control('climb_rate', units='ft/min', opt=True, lower=-3000, upper=3000, ref0=-3000,
+                      ref=3000)
 
     phase.add_control('mach', units=None, opt=False)
 
@@ -67,7 +66,6 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto'):
     phase.add_input_parameter('mass_payload', units='kg')
 
     phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
-    phase.add_path_constraint('alt_rate', units='ft/min', lower=-3000, upper=3000, ref=3000)
 
     p.model.connect('assumptions.S', 'phase0.input_parameters:S')
     p.model.connect('assumptions.mass_empty', 'phase0.input_parameters:mass_empty')
@@ -83,9 +81,9 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto'):
     p['phase0.t_duration'] = 3600.0
     p['phase0.states:range'] = phase.interpolate(ys=(0, 1000.0), nodes='state_input')
     p['phase0.states:mass_fuel'] = phase.interpolate(ys=(30000, 0), nodes='state_input')
+    p['phase0.states:alt'][:] = 10.0
 
     p['phase0.controls:mach'][:] = 0.8
-    p['phase0.controls:alt'][:] = 10.0
 
     p['assumptions.S'] = 427.8
     p['assumptions.mass_empty'] = 0.15E6
@@ -101,9 +99,9 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto'):
     plt.suptitle('altitude vs time')
     plt.figure()
     plt.plot(phase.get_values('time', nodes='all'),
-             phase.get_values('alt_rate', nodes='all', units='ft/min'), 'ro')
-    plt.plot(exp_out.get_values('time'), exp_out.get_values('alt_rate', units='ft/min'), 'b-')
-    plt.suptitle('altitude rate vs time')
+             phase.get_values('climb_rate', nodes='all', units='ft/min'), 'ro')
+    plt.plot(exp_out.get_values('time'), exp_out.get_values('climb_rate', units='ft/min'), 'b-')
+    plt.suptitle('climb rate vs time')
     plt.figure()
     plt.plot(phase.get_values('time', nodes='all'), phase.get_values('mass_fuel', nodes='all'),
              'ro')
@@ -131,10 +129,10 @@ def ex_aircraft_steady_flight(optimizer='SLSQP', transcription='gauss-lobatto'):
     print(phase.get_values('alt', nodes='all').T)
 
     print('alt_rate')
-    print(phase.get_values('alt_rate', nodes='all').T)
+    print(phase.get_values('climb_rate', nodes='all').T)
 
-    print('alt_rate2')
-    print(phase.get_values('alt_rate2', nodes='all').T)
+    print('climb_rate_rate')
+    print(phase.get_values('climb_rate_rate', nodes='all').T)
 
     print('range')
     print(phase.get_values('range', nodes='all').T)

dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py

@@ -58,10 +58,10 @@ def test_cruise_results(self):
                                 defect_scaler=0.01)
         phase.set_state_options('mass_fuel', fix_final=True, upper=20000.0, lower=0.0,
                                 scaler=1.0E-4, defect_scaler=1.0E-2)
-
-        phase.add_control('alt', units='km', opt=False, rate_param='climb_rate')
+        phase.set_state_options('alt', units='km', fix_initial=True)
 
         phase.add_control('mach', units=None, opt=False)
+        phase.add_control('climb_rate', units='m/s', opt=False)
 
         phase.add_input_parameter('S', units='m**2')
         phase.add_input_parameter('mass_empty', units='kg')
@@ -84,8 +84,9 @@ def test_cruise_results(self):
         p['phase0.t_duration'] = 1.515132 * 3600.0
         p['phase0.states:range'] = phase.interpolate(ys=(0, 1296.4), nodes='state_input')
         p['phase0.states:mass_fuel'] = phase.interpolate(ys=(12236.594555, 0), nodes='state_input')
+        p['phase0.states:alt'] = 5.0
         p['phase0.controls:mach'] = 0.8
-        p['phase0.controls:alt'] = 5.0
+        p['phase0.controls:climb_rate'] = 0.0
 
         p['assumptions.S'] = 427.8
         p['assumptions.mass_empty'] = 0.15E6

dymos/examples/aircraft_steady_flight/test/test_doc_aircraft_steady_flight.py

@@ -62,19 +62,18 @@ def test_steady_aircraft_for_docs(self):
         phase.set_state_options('mass_fuel', units='lbm', fix_initial=True, fix_final=True,
                                 upper=1.5E5, lower=0.0, scaler=1.0E-5, defect_scaler=1.0E-1)
 
-        phase.add_control('alt', units='kft', opt=True, lower=0.0, upper=50.0,
-                          rate_param='climb_rate',
-                          rate_continuity=True, rate_continuity_scaler=1.0,
-                          rate2_continuity=True, rate2_continuity_scaler=1.0, ref=1.0,
-                          fix_initial=True, fix_final=True)
+        phase.set_state_options('alt', units='kft', fix_initial=True, fix_final=True,
+                                upper=60, lower=0)
+
+        phase.add_control('climb_rate', units='ft/min', opt=True, lower=-3000, upper=3000,
+                          ref0=-3000, ref=3000)
 
         phase.add_input_parameter('mach', units=None)
         phase.add_input_parameter('S', units='m**2')
         phase.add_input_parameter('mass_empty', units='kg')
         phase.add_input_parameter('mass_payload', units='kg')
 
         phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
-        phase.add_path_constraint('alt_rate', units='ft/min', lower=-3000, upper=3000, ref=3000)
 
         p.model.connect('assumptions.mach', 'phase0.input_parameters:mach')
         p.model.connect('assumptions.S', 'phase0.input_parameters:S')
@@ -91,8 +90,7 @@ def test_steady_aircraft_for_docs(self):
         p['phase0.t_duration'] = 3600.0
         p['phase0.states:range'] = phase.interpolate(ys=(0, 1000.0), nodes='state_input')
         p['phase0.states:mass_fuel'] = phase.interpolate(ys=(30000, 0), nodes='state_input')
-
-        p['phase0.controls:alt'][:] = 10.0
+        p['phase0.states:alt'][:] = 10.0
 
         p['assumptions.mach'][:] = 0.8
         p['assumptions.S'] = 427.8
@@ -122,10 +120,10 @@ def test_steady_aircraft_for_docs(self):
 
         plt.figure()
         plt.plot(phase.get_values('time', nodes='all'),
-                 phase.get_values('alt_rate', nodes='all', units='ft/min'),
+                 phase.get_values('climb_rate', nodes='all', units='ft/min'),
                  'ro')
         plt.plot(exp_out.get_values('time'),
-                 exp_out.get_values('alt_rate', units='ft/min'),
+                 exp_out.get_values('climb_rate', units='ft/min'),
                  'b-')
 
         plt.suptitle('Climb Rate vs Time')

dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py

@@ -3,9 +3,6 @@
 import os
 import unittest
 
-import matplotlib
-matplotlib.use('Agg')
-
 from openmdao.utils.assert_utils import assert_rel_error
 
 from dymos.examples.aircraft_steady_flight.ex_aircraft_steady_flight import \

dymos/examples/brachistochrone/brachistochrone_vector_states_ode.py

@@ -0,0 +1,78 @@
+from __future__ import print_function, division, absolute_import
+
+import numpy as np
+from openmdao.api import ExplicitComponent
+from dymos import declare_time, declare_state, declare_parameter
+
+
+@declare_time(units='s')
+@declare_state('pos', rate_source='pos_dot', units='m', shape=(2,))
+@declare_state('v', rate_source='vdot', targets=['v'], units='m/s')
+@declare_parameter('theta', targets=['theta'], units='rad')
+@declare_parameter('g', units='m/s**2', targets=['g'])
+class BrachistochroneVectorStatesODE(ExplicitComponent):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        # Inputs
+        self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')
+
+        self.add_input('g', val=9.80665 * np.ones(nn), desc='grav. acceleration', units='m/s/s')
+
+        self.add_input('theta', val=np.zeros(nn), desc='angle of wire', units='rad')
+
+        self.add_output('pos_dot', val=np.zeros((nn, 2)), desc='velocity components', units='m/s')
+
+        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2')
+
+        self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant',
+                        units='m/s')
+
+        # Setup partials
+        arange = np.arange(self.options['num_nodes'])
+
+        self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)
+        self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange)
+
+        rs = np.arange(2*nn, dtype=int)
+        cs = np.repeat(np.arange(nn, dtype=int), 2)
+        self.declare_partials(of='pos_dot', wrt='v', rows=rs, cols=cs, val=5)
+        self.declare_partials(of='pos_dot', wrt='theta', rows=rs, cols=cs, val=7)
+
+        self.declare_partials(of='check', wrt='v', rows=arange, cols=arange)
+        self.declare_partials(of='check', wrt='theta', rows=arange, cols=arange)
+
+    def compute(self, inputs, outputs):
+        theta = inputs['theta']
+        cos_theta = np.cos(theta)
+        sin_theta = np.sin(theta)
+        g = inputs['g']
+        v = inputs['v']
+
+        outputs['vdot'] = g * cos_theta
+        outputs['pos_dot'][:, 0] = v * sin_theta
+        outputs['pos_dot'][:, 1] = -v * cos_theta
+        outputs['check'] = v / sin_theta
+
+    def compute_partials(self, inputs, jacobian):
+        theta = inputs['theta']
+        cos_theta = np.cos(theta)
+        sin_theta = np.sin(theta)
+        g = inputs['g']
+        v = inputs['v']
+
+        jacobian['vdot', 'g'] = cos_theta
+        jacobian['vdot', 'theta'] = -g * sin_theta
+
+        jacobian['pos_dot', 'v'][0::2] = sin_theta
+        jacobian['pos_dot', 'v'][1::2] = -cos_theta
+
+        jacobian['pos_dot', 'theta'][0::2] = v * cos_theta
+        jacobian['pos_dot', 'theta'][1::2] = v * sin_theta
+
+        jacobian['check', 'v'] = 1 / sin_theta
+        jacobian['check', 'theta'] = -v * cos_theta / sin_theta**2