Simplify ssto

dymos/examples/ssto/doc/test_doc_ssto_earth.py

@@ -1,8 +1,6 @@
 import unittest
 
 import matplotlib.pyplot as plt
-plt.switch_backend('Agg')
-plt.style.use('ggplot')
 
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.testing_utils import use_tempdirs
@@ -23,7 +21,7 @@ def test_doc_ssto_earth(self):
         #
         p = om.Problem(model=om.Group())
         p.driver = om.pyOptSparseDriver()
-        p.driver.declare_coloring()
+        p.driver.declare_coloring(tol=1.0E-12)
 
         from dymos.examples.ssto.launch_vehicle_ode import LaunchVehicleODE
 
@@ -33,7 +31,6 @@ def test_doc_ssto_earth(self):
         traj = dm.Trajectory()
 
         phase = dm.Phase(ode_class=LaunchVehicleODE,
-                         ode_init_kwargs={'central_body': 'earth'},
                          transcription=dm.GaussLobatto(num_segments=12, order=3, compressed=False))
 
         traj.add_phase('phase0', phase)
@@ -44,19 +41,19 @@ def test_doc_ssto_earth(self):
         #
         phase.set_time_options(initial_bounds=(0, 0), duration_bounds=(10, 500))
 
-        phase.add_state('x', fix_initial=True, ref=1.0E5, defect_ref=1.0,
-                        rate_source='eom.xdot', units='m')
-        phase.add_state('y', fix_initial=True, ref=1.0E5, defect_ref=1.0,
-                        rate_source='eom.ydot', targets=['atmos.y'], units='m')
-        phase.add_state('vx', fix_initial=True, ref=1.0E3, defect_ref=1.0,
-                        rate_source='eom.vxdot', targets=['eom.vx'], units='m/s')
-        phase.add_state('vy', fix_initial=True, ref=1.0E3, defect_ref=1.0,
-                        rate_source='eom.vydot', targets=['eom.vy'], units='m/s')
-        phase.add_state('m', fix_initial=True, ref=1.0E3, defect_ref=1.0,
-                        rate_source='eom.mdot', targets=['eom.m'], units='kg')
-
-        phase.add_control('theta', units='rad', lower=-1.57, upper=1.57, targets=['eom.theta'])
-        phase.add_parameter('thrust', units='N', opt=False, val=2100000.0, targets=['eom.thrust'])
+        phase.add_state('x', fix_initial=True, ref=1.0E5, defect_ref=10000.0,
+                        rate_source='xdot')
+        phase.add_state('y', fix_initial=True, ref=1.0E5, defect_ref=10000.0,
+                        rate_source='ydot')
+        phase.add_state('vx', fix_initial=True, ref=1.0E3, defect_ref=1000.0,
+                        rate_source='vxdot')
+        phase.add_state('vy', fix_initial=True, ref=1.0E3, defect_ref=1000.0,
+                        rate_source='vydot')
+        phase.add_state('m', fix_initial=True, ref=1.0E3, defect_ref=100.0,
+                        rate_source='mdot')
+
+        phase.add_control('theta', units='rad', lower=-1.57, upper=1.57, targets=['theta'])
+        phase.add_parameter('thrust', units='N', opt=False, val=2100000.0, targets=['thrust'])
 
         #
         # Set the options for our constraints and objective

dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py

@@ -202,16 +202,6 @@ def compute_partials(self, inputs, jacobian):
                 jacobian['theta', 'b_ctrl'] = 1.0 / denom
                 jacobian['theta', 'time_phase'] = a / denom
 
-        # @dm.declare_time(units='s', targets=['guidance.time_phase'])
-        # @dm.declare_state('x', rate_source='eom.xdot', units='m')
-        # @dm.declare_state('y', rate_source='eom.ydot', units='m')
-        # @dm.declare_state('vx', rate_source='eom.vxdot', targets=['eom.vx'], units='m/s')
-        # @dm.declare_state('vy', rate_source='eom.vydot', targets=['eom.vy'], units='m/s')
-        # @dm.declare_state('m', rate_source='eom.mdot', targets=['eom.m'], units='kg')
-        # @dm.declare_parameter('thrust', targets=['eom.thrust'], units='N')
-        # @dm.declare_parameter('a_ctrl', targets=['guidance.a_ctrl'], units='1/s')
-        # @dm.declare_parameter('b_ctrl', targets=['guidance.b_ctrl'], units=None)
-        # @dm.declare_parameter('Isp', targets=['eom.Isp'], units='s')
         class LaunchVehicleLinearTangentODE(om.Group):
             """
             The LaunchVehicleLinearTangentODE for this case consists of a guidance component and

dymos/examples/ssto/launch_vehicle_ode.py

@@ -1,34 +1,121 @@
 import openmdao.api as om
+import numpy as np
 
-from .log_atmosphere_comp import LogAtmosphereComp
-from .launch_vehicle_2d_eom_comp import LaunchVehicle2DEOM
 
-
-class LaunchVehicleODE(om.Group):
+class LaunchVehicleODE(om.ExplicitComponent):
 
     def initialize(self):
         self.options.declare('num_nodes', types=int,
                              desc='Number of nodes to be evaluated in the RHS')
 
-        self.options.declare('central_body', values=['earth', 'moon'], default='earth',
-                             desc='The central gravitational body for the launch vehicle.')
+        self.options.declare('g', types=float, default=9.80665,
+                             desc='Gravitational acceleration, m/s**2')
+
+        self.options.declare('rho_ref', types=float, default=1.225,
+                             desc='Reference atmospheric density, kg/m**3')
+
+        self.options.declare('h_scale', types=float, default=8.44E3,
+                             desc='Reference altitude, m')
+
+        self.options.declare('CD', types=float, default=0.5,
+                             desc='coefficient of drag')
+
+        self.options.declare('S', types=float, default=7.069,
+                             desc='aerodynamic reference area (m**2)')
 
     def setup(self):
         nn = self.options['num_nodes']
-        cb = self.options['central_body']
 
-        if cb == 'earth':
-            rho_ref = 1.225
-            h_scale = 8.44E3
-        elif cb == 'moon':
-            rho_ref = 0.0
-            h_scale = 1.0
-        else:
-            raise RuntimeError('Unrecognized value for central_body: {0}'.format(cb))
+        self.add_input('y',
+                       val=np.zeros(nn),
+                       desc='altitude',
+                       units='m')
+
+        self.add_input('vx',
+                       val=np.zeros(nn),
+                       desc='x velocity',
+                       units='m/s')
+
+        self.add_input('vy',
+                       val=np.zeros(nn),
+                       desc='y velocity',
+                       units='m/s')
+
+        self.add_input('m',
+                       val=np.zeros(nn),
+                       desc='mass',
+                       units='kg')
+
+        self.add_input('theta',
+                       val=np.zeros(nn),
+                       desc='pitch angle',
+                       units='rad')
+
+        self.add_input('thrust',
+                       val=2100000 * np.ones(nn),
+                       desc='thrust',
+                       units='N')
+
+        self.add_input('Isp',
+                       val=265.2 * np.ones(nn),
+                       desc='specific impulse',
+                       units='s')
+        # Outputs
+        self.add_output('xdot',
+                        val=np.zeros(nn),
+                        desc='velocity component in x',
+                        units='m/s')
+
+        self.add_output('ydot',
+                        val=np.zeros(nn),
+                        desc='velocity component in y',
+                        units='m/s')
+
+        self.add_output('vxdot',
+                        val=np.zeros(nn),
+                        desc='x acceleration magnitude',
+                        units='m/s**2')
+
+        self.add_output('vydot',
+                        val=np.zeros(nn),
+                        desc='y acceleration magnitude',
+                        units='m/s**2')
+
+        self.add_output('mdot',
+                        val=np.zeros(nn),
+                        desc='mass rate of change',
+                        units='kg/s')
+
+        self.add_output('rho',
+                        val=np.zeros(nn),
+                        desc='density',
+                        units='kg/m**3')
+
+        # Setup partials
+        # Complex-step derivatives
+        self.declare_coloring(wrt='*', method='cs', show_sparsity=True)
+
+    def compute(self, inputs, outputs):
+
+        theta = inputs['theta']
+        cos_theta = np.cos(theta)
+        sin_theta = np.sin(theta)
+        vx = inputs['vx']
+        vy = inputs['vy']
+        m = inputs['m']
+        F_T = inputs['thrust']
+        Isp = inputs['Isp']
+        y = inputs['y']
 
-        self.add_subsystem('atmos',
-                           LogAtmosphereComp(num_nodes=nn, rho_ref=rho_ref, h_scale=h_scale))
+        g = self.options['g']
+        rho_ref = self.options['rho_ref']
+        h_scale = self.options['h_scale']
 
-        self.add_subsystem('eom', LaunchVehicle2DEOM(num_nodes=nn, central_body=cb))
+        CDA = self.options['CD'] * self.options['S']
 
-        self.connect('atmos.rho', 'eom.rho')
+        outputs['rho'] = rho_ref * np.exp(-y / h_scale)
+        outputs['xdot'] = vx
+        outputs['ydot'] = vy
+        outputs['vxdot'] = (F_T * cos_theta - 0.5 * CDA * outputs['rho'] * vx**2) / m
+        outputs['vydot'] = (F_T * sin_theta - 0.5 * CDA * outputs['rho'] * vy**2) / m - g
+        outputs['mdot'] = -F_T / (g * Isp)