diff --git a/dymos/examples/ssto/doc/test_doc_ssto_earth.py b/dymos/examples/ssto/doc/test_doc_ssto_earth.py index 236846846..d7ad841a1 100644 --- a/dymos/examples/ssto/doc/test_doc_ssto_earth.py +++ b/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 diff --git a/dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py b/dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py index bacf34e33..1730cf89b 100644 --- a/dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py +++ b/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 diff --git a/dymos/examples/ssto/launch_vehicle_2d_eom_comp.py b/dymos/examples/ssto/launch_vehicle_2d_eom_comp.py deleted file mode 100644 index 66eda7f1e..000000000 --- a/dymos/examples/ssto/launch_vehicle_2d_eom_comp.py +++ /dev/null @@ -1,155 +0,0 @@ -import numpy as np - -import openmdao.api as om -_g = {'earth': 9.80665, - 'moon': 1.61544} - - -class LaunchVehicle2DEOM(om.ExplicitComponent): - - def initialize(self): - self.options.declare('num_nodes', types=int) - - self.options.declare('central_body', values=['earth', 'moon'], default='earth', - desc='The central graviational body for the launch vehicle.') - - 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'] - - # Inputs - 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('rho', - val=np.zeros(nn), - desc='density', - units='kg/m**3') - - 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') - - # Setup partials - ar = np.arange(self.options['num_nodes']) - - self.declare_partials(of='xdot', wrt='vx', rows=ar, cols=ar, val=1.0) - self.declare_partials(of='ydot', wrt='vy', rows=ar, cols=ar, val=1.0) - - self.declare_partials(of='vxdot', wrt='vx', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='m', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='theta', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='rho', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='thrust', rows=ar, cols=ar) - - self.declare_partials(of='vydot', wrt='m', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='theta', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='vy', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='rho', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='thrust', rows=ar, cols=ar) - - self.declare_partials(of='mdot', wrt='thrust', rows=ar, cols=ar) - self.declare_partials(of='mdot', wrt='Isp', rows=ar, cols=ar) - - 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'] - rho = inputs['rho'] - F_T = inputs['thrust'] - Isp = inputs['Isp'] - - g = _g[self.options['central_body']] - CDA = self.options['CD'] * self.options['S'] - - outputs['xdot'] = vx - outputs['ydot'] = vy - outputs['vxdot'] = (F_T * cos_theta - 0.5 * CDA * rho * vx**2) / m - outputs['vydot'] = (F_T * sin_theta - 0.5 * CDA * rho * vy**2) / m - g - outputs['mdot'] = -F_T / (g * Isp) - - def compute_partials(self, inputs, jacobian): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - m = inputs['m'] - vx = inputs['vx'] - vy = inputs['vy'] - rho = inputs['rho'] - F_T = inputs['thrust'] - Isp = inputs['Isp'] - - g = _g[self.options['central_body']] - CDA = self.options['CD'] * self.options['S'] - - jacobian['vxdot', 'vx'] = -CDA * rho * vx / m - jacobian['vxdot', 'rho'] = -0.5 * CDA * vx**2 / m - jacobian['vxdot', 'm'] = -(F_T * cos_theta - 0.5 * CDA * rho * vx**2) / m**2 - jacobian['vxdot', 'theta'] = -(F_T / m) * sin_theta - jacobian['vxdot', 'thrust'] = cos_theta / m - - jacobian['vydot', 'vy'] = -CDA * rho * vy / m - jacobian['vydot', 'rho'] = -0.5 * CDA * vy**2 / m - jacobian['vydot', 'm'] = -(F_T * sin_theta - 0.5 * CDA * rho * vy**2) / m**2 - jacobian['vydot', 'theta'] = (F_T / m) * cos_theta - jacobian['vydot', 'thrust'] = sin_theta / m - - jacobian['mdot', 'thrust'] = -1.0 / (g * Isp) - jacobian['mdot', 'Isp'] = F_T / (g * Isp**2) diff --git a/dymos/examples/ssto/launch_vehicle_ode.py b/dymos/examples/ssto/launch_vehicle_ode.py index d238d0424..20ad7b8c0 100644 --- a/dymos/examples/ssto/launch_vehicle_ode.py +++ b/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) diff --git a/dymos/examples/ssto/log_atmosphere_comp.py b/dymos/examples/ssto/log_atmosphere_comp.py deleted file mode 100644 index 3f71ddc1d..000000000 --- a/dymos/examples/ssto/log_atmosphere_comp.py +++ /dev/null @@ -1,35 +0,0 @@ -import numpy as np - -import openmdao.api as om - - -class LogAtmosphereComp(om.ExplicitComponent): - def initialize(self): - self.options.declare('num_nodes', types=int) - self.options.declare('rho_ref', types=float, default=1.225, - desc='reference density, kg/m**3') - self.options.declare('h_scale', types=float, default=8.44E3, - desc='reference altitude, m') - - def setup(self): - nn = self.options['num_nodes'] - - self.add_input('y', val=np.zeros(nn), desc='altitude', units='m') - - self.add_output('rho', val=np.zeros(nn), desc='density', units='kg/m**3') - - # Setup partials - arange = np.arange(self.options['num_nodes']) - self.declare_partials(of='rho', wrt='y', rows=arange, cols=arange, val=1.0) - - def compute(self, inputs, outputs): - rho_ref = self.options['rho_ref'] - h_scale = self.options['h_scale'] - y = inputs['y'] - outputs['rho'] = rho_ref * np.exp(-y / h_scale) - - def compute_partials(self, inputs, jacobian): - rho_ref = self.options['rho_ref'] - h_scale = self.options['h_scale'] - y = inputs['y'] - jacobian['rho', 'y'] = -(rho_ref / h_scale) * np.exp(-y / h_scale) diff --git a/dymos/test/test_upgrade_guide.py b/dymos/test/test_upgrade_guide.py index 161ab2a69..ac9912d48 100644 --- a/dymos/test/test_upgrade_guide.py +++ b/dymos/test/test_upgrade_guide.py @@ -30,7 +30,7 @@ def test_parameters(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 @@ -40,7 +40,6 @@ def test_parameters(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) @@ -51,19 +50,19 @@ def test_parameters(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) + phase.add_parameter('thrust', units='N', opt=False, val=2100000.0) # # Set the options for our constraints and objective diff --git a/mkdocs/docs/examples/ssto_earth/ssto_earth.md b/mkdocs/docs/examples/ssto_earth/ssto_earth.md index 786862610..72c85c4c0 100644 --- a/mkdocs/docs/examples/ssto_earth/ssto_earth.md +++ b/mkdocs/docs/examples/ssto_earth/ssto_earth.md @@ -37,35 +37,21 @@ and the final conditions are m_f &= \rm{free} \end{align} -## Component and Group Definitions +## Defining the ODE -The ODE system for this problem consists of two components. The -atmosphere component computes density ($\rho$). The eom component -computes the state rates. +Generally, one could define the ODE system as a composite group of multile components. +The atmosphere component computes density ($\rho$). +The eom component computes the state rates. +Decomposing the ODE into smaller calculations makes it easier to derive the analytic derivatives. -![The XDSM diagram for the ODE system in the SSTO problem.](ssto_xdsm.png) +![The notional XDSM diagram for the ODE system in the SSTO problem.](ssto_xdsm.png) -The unconnected inputs to the EOM at the top of the diagram are provided by -the Dymos phase as states, controls, or time values. The outputs, -including the state rates, are shown on the right side of the diagram. -The Dymos phases use state rate values to ensure that the integration -technique satisfies the dynamics of the system. +However, for this example we will demonstrate the use of complex-step differentiation and define the ODE as a single component. +This saves time up front in the deevlopment of the ODE at a minor cost in execution time. -=== "launch_vehicle_2d_eom_comp.py" -{{ inline_source('dymos.examples.ssto.launch_vehicle_2d_eom_comp', -include_def=True, -include_docstring=True, -indent_level=0) -}} - -=== "log_atmosphere_comp.py" -{{ inline_source('dymos.examples.ssto.log_atmosphere_comp', -include_def=True, -include_docstring=True, -indent_level=0) -}} - -### ## Defining the ODE +The unconnected inputs to the EOM at the top of the diagram are provided by the Dymos phase as states, controls, or time values. +The outputs, including the state rates, are shown on the right side of the diagram. +The Dymos phases use state rate values to ensure that the integration technique satisfies the dynamics of the system. === "launch_vehicle_ode.py" {{ inline_source('dymos.examples.ssto.launch_vehicle_ode', diff --git a/mkdocs/docs/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.md b/mkdocs/docs/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.md index b1176c156..b140f569d 100644 --- a/mkdocs/docs/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.md +++ b/mkdocs/docs/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.md @@ -14,7 +14,7 @@ just two scalar parameters, $a$ and $b$. $$\theta = \arctan{\left(a * t + b\right)}$$ -Implementing this modified constrol scheme requires only a few changes. +Implementing this modified control scheme requires only a few changes. Rather than declaring $\theta$ as a controllable parameter for the ODE system, we implement a new component, _LinearTangentGuidanceComp_ that accepts $a$ and $b$ as parameters to be optimized. It calculates $\theta$, which is then connected to the equations of motion component.