{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"active-ipynb",
"remove-input",
"remove-output"
]
},
"outputs": [],
"source": [
"# This cell is mandatory in all Dymos documentation notebooks.\n",
"missing_packages = []\n",
"try:\n",
" import openmdao.api as om\n",
"except ImportError:\n",
" if 'google.colab' in str(get_ipython()):\n",
" !python -m pip install openmdao[notebooks]\n",
" else:\n",
" missing_packages.append('openmdao')\n",
"try:\n",
" import dymos as dm\n",
"except ImportError:\n",
" if 'google.colab' in str(get_ipython()):\n",
" !python -m pip install dymos\n",
" else:\n",
" missing_packages.append('dymos')\n",
"try:\n",
" import pyoptsparse\n",
"except ImportError:\n",
" if 'google.colab' in str(get_ipython()):\n",
" !pip install -q condacolab\n",
" import condacolab\n",
" condacolab.install_miniconda()\n",
" !conda install -c conda-forge pyoptsparse\n",
" else:\n",
" missing_packages.append('pyoptsparse')\n",
"if missing_packages:\n",
" raise EnvironmentError('This notebook requires the following packages '\n",
" 'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Commercial Aircraft Range Maximization by Differential Inclusion\n",
"\n",
"In this example we seek to maximize the range of a commercial aircraft\n",
"flying at a constant Mach number with a fixed initial fuel load. The\n",
"initial and final altitudes are constrained to 10000 ft (that is, in\n",
"this simplified example we're ignoring the aircraft behavior near the\n",
"ground).\n",
"\n",
"This example is simplified in several ways. The Mach number is fixed at\n",
"0.8. When only considering fuel burn, the aircraft will tend to fly\n",
"slower to conserve fuel, ignoring the running costs of being airborne.\n",
"\n",
"For a more detailed optimization of commercial aircraft performance, see\n",
"Betts {cite}`ex_commercial_aircraft-betts1995application`.\n",
"\n",
"## Differential Inclusion\n",
"\n",
"In this case we also use the steady flight assumption to allow the\n",
"problem to be parameterized with a different control set. Rather than\n",
"providing throttle/thrust and angle of attack as control variables, we\n",
"use a _differential inclusion_ approach. Whereas altitude and velocity\n",
"are often treated as state variables to be integrated, here we take a\n",
"slightly different approach.\n",
"\n",
"Rather than using the differential equations to compute rates for the\n",
"flight path angle and true airspeed, we use a nonlinear solver to\n",
"determine the alpha and thrust time-history which is needed to make the\n",
"given altitude and airspeed time-history possible.\n",
"This technique is known as differential inclusion.\n",
"It was demonstrated by Seywald {cite}`ex_commercial_aircraft-Seywald1994`.\n",
"Since the collocation techniques in Dymos are based on polynomial approximations, rates of the control variables can be easily calculated, which can then be applied to the reordered dynamics equations to solve for angle of attack and thrust.\n",
"\n",
"The use of collocation/psuedospectral techniques in solving problems via\n",
"differential inclusion has been examined by Fahroo and Ross {cite}`ex_commercial_aircraft-fahroo2001second`\n",
"and Kumar and Seywald {cite}`ex_commercial_aircraft-Kumar1996`. Since the differential inclusion\n",
"approach demonstrated here relies on nonlinear solvers running within\n",
"the ODE model, obtaining accurate derivatives is paramount. Computing\n",
"derivatives with finite differences in the presence of nonlinear solvers\n",
"can lead to inaccurate derivatives based on the properties of the\n",
"solver. With OpenMDAO's unified derivatives framework, the sensitivity to\n",
"solver tolerances is greatly reduced.\n",
"\n",
"### Pitfalls of Differential Inclusion\n",
"\n",
"Differential inclusion works well in this case but it is highly\n",
"dependent on the parameterization of the problem. For instance, it is\n",
"possible for the optimizer to suggest a flight profile that requires\n",
"thrust well beyond the capabilities of the vehicle. In our case, the\n",
"throttle parameter $\\tau$ is path constrained within the optimizer, so\n",
"that the optimizer can make the solution physically attainable by our\n",
"model.\n",
"\n",
"We have two options for how to tackle this problem. We could specify\n",
"altitude as a control variable, and use it's approximate derivative\n",
"($alt\\_rate$) to provide the altitude rate for the flight\n",
"equilibrium conditions. In practice, however, this noisy derivative of\n",
"altitude can cause numerical difficulties. Other times, the optimizer\n",
"will shape the control to satisfy our equilibrium conditions at the\n",
"collocation nodes but be wildly wrong in between these nodes. One simple\n",
"way to get around this is to still integrate altitude, but to provide\n",
"$climb\\_rate$ as the control variable. This makes it easy to\n",
"constraint $climb\\_rate$ with bounds constraints, and the\n",
"time-history of altitude will be smoother since it is the time integral\n",
"of the rate of climb.\n",
"\n",
"#### Solver Options\n",
"\n",
"In general, _openmdao.api.AnalysisError_ should be raised\n",
"when the design variables provided by the optimizer result in\n",
"nonconvergence of the solvers. This informs pyoptsparse to backtrack\n",
"along its line search for a solution. This error can be raised by the\n",
"_compute_ method of a component. The nonlinear solver can\n",
"also raise this error if it fails to converge by setting the following\n",
"options:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```python \n",
"self.nonlinear_solver.options['maxiter'] = 150\n",
"self.nonlinear_solver.options['err_on_maxiter'] = True\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the course of solving, an outer NewtonSolver won't necessarily force\n",
"subsystems to be solved, instead relying on the derivative at the at the\n",
"previous solution. In the context of OpenMDAO, even explicit components\n",
"constitute subsystems. This can lead to poorly converged solutions.\n",
"\n",
"Setting the _solve\\_subsystems_ option to True will force\n",
"subsystems to be updated, making sure that all derivative information is\n",
"up-to-date."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```python\n",
"self.nonlinear_solver.options['solve_subsystems'] = True\n",
"self.nonlinear_solver.options['max_sub_solves'] = 10\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### The Impact of Choice of Transcription\n",
"\n",
"The high-order Gauss-Lobatto method used by Dymos evaluates the ordinary\n",
"differential equations in a two-step process. First, the ODEs are\n",
"evaluated at the *state discretization* nodes. Using the state values\n",
"and derivatives, values and rates at the *collocation* nodes are then\n",
"interpolated. If a large rate is found at one of the state\n",
"discretization nodes, the interpolation at the collocation node can be\n",
"outside of the expected range, leading to nonconvergence of solvers\n",
"within the ODEs.\n",
"\n",
"There is also an overhead cost to invoking solvers. While differential\n",
"inclusion reduces the number of design variables for the optimizer, it\n",
"also slows down each evaluation of the ODE somewhat. Since evaluating\n",
"the collocation constraints in the Gauss-Lobatto method requires two\n",
"evaluations of the ODE, it typically suffers greater slowdown that the\n",
"Radau Pseudospectral method, which evaluates the defects of the ODE with\n",
"a single evaluation that encompasses all the nodes.\n",
"\n",
"For these reasons the Radau Pseudospectral transcription is typically\n",
"preferable when using a differential inclusion approach to problems.\n",
"\n",
"## Problem Formulation\n",
"\n",
"### State Variables\n",
"\n",
"|Name |Description |Fixed Initial Value |Fixed Final Value|\n",
"|------------|------------------------------|---------------------|-----------------|\n",
"|range |distance flown along ground |True (0 NM) |False |\n",
"|mass\\_fuel |mass of fuel aboard aircraft |True (30000 lbm) |True (0 lbm) |\n",
"|alt |aircraft altitude |True (10000 ft) |True (10000 ft) |\n",
"\n",
"### Dynamic Controls\n",
"\n",
"|Name |Description |Optimized |Fixed Initial Value |Fixed Final Value|\n",
"|-------------|------------------------|-----------|---------------------|-----------------|\n",
"|climb\\_rate |aircraft rate of climb |True |False |False |\n",
"\n",
"### Parameters\n",
"\n",
"|Name |Description |Optimized |\n",
"|---------------|----------------------------|----------------------|\n",
"|mach |mach number |False (0.8) |\n",
"|S |aerodynamic reference area |False (427.8 m\\*\\*2) |\n",
"|mass\\_empty |aircraft empty mass |False (330693.393 lbm)|\n",
"|mass\\_payload |aircraft payload mass |False (74100 lbm) |\n",
"\n",
"### Objective\n",
"\n",
"|Name | Description | Location | Minimized or Maximized |\n",
"|-------|-------------|----------|------------------------|\n",
"|r | range | final | Maximized |\n",
"\n",
"### Nonlinear Path Constraints\n",
"\n",
"|Name | Description |Lower |Upper|\n",
"|------|-----------------------------|-------|-----|\n",
"|tau | engine throttle parameter |0.01 |1.0 |\n",
"\n",
"### Nonlinear Boundary Constraints\n",
"\n",
"None\n",
"\n",
"## Models\n",
"\n",
"### Atmosphere\n",
"\n",
"This problem uses an analytic fit to the 1976 standard atmosphere.\n",
"\n",
"|Name |Description |Input or Output|\n",
"|------|----------------------|---------------|\n",
"|alt |altitude (m) |input |\n",
"|pres |static pressure (Pa) |output |\n",
"|temp |temperature (K) |output |\n",
"|sos |speed of sound (m/s) |output |\n",
"|rho |density (kg/m\\*\\*3) |output |\n",
"\n",
"### True Airspeed\n",
"\n",
"_TrueAirspeedComp_ uses the Mach number, provided as a\n",
"control, and the speed of sound from the atmosphere model to compute the\n",
"true airspeed of the aircraft.\n",
"\n",
"\\begin{align}\n",
" TAS &= mach \\cdot sos\n",
"\\end{align}\n",
"\n",
"|Name |Description |Input or Output|\n",
"|------|----------------------|---------------|\n",
"|mach |Mach number |input |\n",
"|sos |speed of sound (m/s) |input |\n",
"|TAS |true airspeed (m/s) |output |\n",
"\n",
"### Flight Path Angle\n",
"\n",
"_SteadyFlightPathAngleComp_ uses the true airspeed and the\n",
"climb rate, obtained by differentiating the altitude time history at the\n",
"nodes, to compute the flight path angle.\n",
"\n",
"\\begin{align}\n",
" \\gamma &= \\arcsin \\frac{\\dot h}{TAS}\n",
"\\end{align}\n",
"\n",
"|Name |Description |Input or Output|\n",
"|-----------|-------------------------|---------------|\n",
"|TAS |true airspeed (m/s) |input |\n",
"|alt\\_rate |climb rate (m/s) |input |\n",
"|gamma |flight path angle (rad) |output |\n",
"\n",
"### Range Rate\n",
"\n",
"_RangeRateComp_ uses the true airspeed and the flight path\n",
"angle to determine the velocity projected along the ground. This is the\n",
"derivative of the state variable _range_.\n",
"\n",
"\\begin{align}\n",
" \\dot{range} &= TAS \\cos \\gamma\n",
"\\end{align}\n",
"\n",
"|Name |Description |Input or Output|\n",
"|------------|-------------------------|---------------|\n",
"|TAS |true airspeed (m/s) |input |\n",
"|gamma |flight path angle (rad) |input |\n",
"|dXdt:range |range rate (m/s) |output |\n",
"\n",
"### Mass\n",
"\n",
"The component _MassComp_ defined in\n",
"`mass_comp.py` computes the aircraft total mass based on\n",
"its empty mass, payload mass, and current fuel mass. It also computes\n",
"total weight which simplifies some equations later on.\n",
"\n",
"\\begin{align}\n",
" mass_{total} &= mass_{empty} + mass_{payload} + mass_{fuel} \\\\\n",
" W_{total} &= 9.80665 \\, mass_{total}\n",
"\\end{align}\n",
"\n",
"|Name |Description |Input or Output|\n",
"|---------------|---------------------------|---------------|\n",
"|mass\\_empty |aircraft empty mass (kg) |input |\n",
"|mass\\_payload |payload mass (kg) |input |\n",
"|mass\\_fuel |fuel mass (kg) |input |\n",
"|mass\\_total |total aircraft mass (kg) |output |\n",
"|W\\_total |total aircraft weight (N) |output |\n",
"\n",
"### Dynamic Pressure\n",
"\n",
"The _DynamicPressureComp_ computes the dynamic pressure from\n",
"true airspeed and atmospheric density.\n",
"\n",
"\\begin{align}\n",
" q &= \\frac{1}{2} \\rho TAS^2\n",
"\\end{align}\n",
"\n",
"|Name |Description |Input or Output|\n",
"|------|---------------------------------|---------------|\n",
"|TAS |true airspeed (m/s) |input |\n",
"|rho |atmospheric density (kg/m\\*\\*3) |input |\n",
"|q |dynamic pressure (Pa) |output |\n",
"\n",
"### Aerodynamics\n",
"\n",
"The aerodynamics group computes the aerodynamic coefficients and forces\n",
"on the vehicle. It consists of an interpolation component which outputs\n",
"lift, drag, and moment coefficients as a function of Mach number, angle\n",
"of attack, altitude, and tail rotation angle. A second component then\n",
"uses these coefficients, along with dynamic pressure and aerodynamic\n",
"reference area, to compute the lift and drag forces on the vehicle.\n",
"\n",
"The aerodynamics group resides within the flight equilibrium group. As\n",
"that group iterates to find the combination of thrust coefficient, angle\n",
"of attack, and tail rotation angle, aerodynamics needs to update the\n",
"values of the interpolated coefficients and resulting forces.\n",
"\n",
"Organizationally speaking, we logically could have put the dynamic\n",
"pressure component within the aerodynamics group. However, since that\n",
"group doesn't need to be updated with changes in alpha and tail angle,\n",
"it's more efficient to leave it outside of flight equilibrium group.\n",
"\n",
" |Name |Description |Input or Output|\n",
" |-------|---------------------------|---------------|\n",
" |mach |mach number |input |\n",
" |alt |altitude (m) |input |\n",
" |alpha |angle of attack (deg) |input |\n",
" |eta |tail rotation angle (deg) |input |\n",
" |CL |lift coefficient |output |\n",
" |CD |drag coefficient |output |\n",
" |CM |moment coefficient |output |\n",
" |L |lift force (N) |output |\n",
" |D |drag force (N) |output |\n",
"\n",
"```{Note}\n",
"This example uses _openmdao.api.MetaModelStructuredComp_ to\n",
"interpolate aerodynamic properties of the vehicle. This component is\n",
"somewhat easier to use since it is distributed as part of OpenMDAO, but\n",
"it can be significantly slower than alternatives such as\n",
"[SMT](https://github.com/SMTorg/smt).\n",
"```\n",
"\n",
"### Flight Equilibrium\n",
"\n",
"The steady flight equilibrium group uses balances to solve for the angle\n",
"of attack and tail plane rotation angle such that the aircraft is in\n",
"steady flight (the rates of change in flight path angle and true\n",
"airspeed are zero) and the aerodynamic moment in the pitch axis is zero.\n",
"\n",
"Of course, in reality the vehicle will accelerate, but the flight\n",
"profile being modeled is so benign that assuming steady flight at\n",
"discrete points (nodes) in the trajectory is not terribly inaccurate.\n",
"\n",
"The thrust necessary for steady flight is computed by balancing the drag\n",
"equation\n",
"\n",
"\\begin{align}\n",
" C_T &= W_{total} \\frac{sin \\gamma}{q \\cdot S \\cdot \\cos \\alpha} + \\frac{C_D}{\\cos \\alpha}\n",
"\\end{align}\n",
"\n",
"The lift coefficient required for steady flight is found by balancing\n",
"lift and weight:\n",
"\n",
"\\begin{align}\n",
" \\tilde{C_L} &= W_{total} \\frac{cos \\gamma}{q \\cdot S} - C_T \\sin \\alpha\n",
"\\end{align}\n",
"\n",
"Using coefficients in the balance equations is better scaled from a numerical standpoint.\n",
"\n",
"### Propulsion\n",
"\n",
"Having determined the thrust, the propulsion group then computes the\n",
"rate of fuel burn. In addition, by normalizing thrust at any point by\n",
"the maximum possible thrust, we obtain the throttle parameter $\\tau$.\n",
"The propulsion group uses a number of components to perform these\n",
"calculations.\n",
"\n",
"\\begin{align}\n",
" T &= C_T \\cdot q \\cdot S\n",
"\\end{align}\n",
"\n",
"Maximum thrust is computed by multiplying sea-level thrust by the ratio\n",
"of pressure to sea-level atmospheric pressure.\n",
"\n",
"\\begin{align}\n",
" T_{max} &= T_{max,sl} \\frac{P}{P_{sl}}\n",
"\\end{align}\n",
"\n",
"The throttle parameter is then the ratio current thrust to maximum\n",
"possible thrust.\n",
"\n",
"\\begin{align}\n",
" \\tau &= \\frac{T}{T_{max}}\n",
"\\end{align}\n",
"\n",
"The thrust specific fuel consumption is computed as follows:\n",
"\n",
"\\begin{align}\n",
" TSFC &= TSFC_{sl} - 1.5 E - 10 \\cdot 9.80665 \\cdot alt\n",
"\\end{align}\n",
"\n",
"Finally, fuel burn rate is:\n",
"\n",
"\\begin{align}\n",
" \\dot{mass_{fuel}} &= -TSFC \\frac{T}{9.80665}\n",
"\\end{align}\n",
"\n",
"## The ODE System: aircraft_ode.py"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```python\n",
"class AircraftODE(om.Group):\n",
"\n",
" def initialize(self):\n",
" self.options.declare('num_nodes', types=int,\n",
" desc='Number of nodes to be evaluated in the RHS')\n",
"\n",
" def setup(self):\n",
" nn = self.options['num_nodes']\n",
"\n",
" self.add_subsystem(name='mass_comp',\n",
" subsys=MassComp(num_nodes=nn),\n",
" promotes_inputs=['mass_fuel'])\n",
"\n",
" self.connect('mass_comp.W_total', 'flight_equilibrium.W_total')\n",
"\n",
" self.add_subsystem(name='atmos',\n",
" subsys=USatm1976Comp(num_nodes=nn),\n",
" promotes_inputs=[('h', 'alt')])\n",
"\n",
" self.connect('atmos.pres', 'propulsion.pres')\n",
" self.connect('atmos.sos', 'tas_comp.sos')\n",
" self.connect('atmos.rho', 'q_comp.rho')\n",
"\n",
" self.add_subsystem(name='tas_comp',\n",
" subsys=TrueAirspeedComp(num_nodes=nn))\n",
"\n",
" self.connect('tas_comp.TAS',\n",
" ('gam_comp.TAS', 'q_comp.TAS', 'range_rate_comp.TAS'))\n",
"\n",
" self.add_subsystem(name='gam_comp',\n",
" subsys=SteadyFlightPathAngleComp(num_nodes=nn))\n",
"\n",
" self.connect('gam_comp.gam', ('flight_equilibrium.gam', 'range_rate_comp.gam'))\n",
"\n",
" self.add_subsystem(name='q_comp',\n",
" subsys=DynamicPressureComp(num_nodes=nn))\n",
"\n",
" self.connect('q_comp.q', ('aero.q', 'flight_equilibrium.q', 'propulsion.q'))\n",
"\n",
" self.add_subsystem(name='flight_equilibrium',\n",
" subsys=SteadyFlightEquilibriumGroup(num_nodes=nn),\n",
" promotes_inputs=['aero.*', 'alt'],\n",
" promotes_outputs=['aero.*'])\n",
"\n",
" self.connect('flight_equilibrium.CT', 'propulsion.CT')\n",
"\n",
" self.add_subsystem(name='propulsion', subsys=PropulsionGroup(num_nodes=nn),\n",
" promotes_inputs=['alt'])\n",
"\n",
" self.add_subsystem(name='range_rate_comp', subsys=RangeRateComp(num_nodes=nn))\n",
"\n",
" # We promoted multiple inputs to the group named 'alt'.\n",
" # In order to automatically create a source variable for 'alt', we must specify their units\n",
" # since they have different units in different components.\n",
" self.set_input_defaults('alt', val=np.ones(nn,), units='m')\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case the system has only two integrated states:\n",
"_range_ and _mass\\_fuel_. There are six parameters.\n",
"Two of them (_alt_ and _climb\\_rate_) will be varied dynamically in the phase\n",
"The other four (_mach_, _S_, _mass\\_empty_, and _mass\\_payload_), will be set to fixed values as non-optimized parameters.\n",
"More details on the various models involved can be found in the examples code.\n",
"\n",
"## Building the problem\n",
"\n",
"In the following code we define and solve the optimal control problem.\n",
"Note that we demonstrate the use of externally-connected design parameters in this case.\n",
"The four parameters are connected to a source provided by the _assumptions_ IndepVarComp."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"remove-input",
"hide-output"
]
},
"outputs": [],
"source": [
"om.display_source(\"dymos.examples.aircraft_steady_flight.aircraft_ode\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"active-ipynb",
"remove-input",
"remove-output"
]
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"import openmdao.api as om\n",
"import dymos as dm\n",
"\n",
"from dymos.examples.aircraft_steady_flight.aircraft_ode import AircraftODE\n",
"from dymos.examples.plotting import plot_results\n",
"from dymos.utils.lgl import lgl\n",
"\n",
"from bokeh.io import output_notebook\n",
"\n",
"p = om.Problem(model=om.Group())\n",
"p.driver = om.pyOptSparseDriver()\n",
"p.driver.options['optimizer'] = 'IPOPT'\n",
"p.driver.opt_settings['print_level'] = 0\n",
"p.driver.declare_coloring()\n",
"\n",
"num_seg = 15\n",
"seg_ends, _ = lgl(num_seg + 1)\n",
"\n",
"traj = p.model.add_subsystem('traj', dm.Trajectory())\n",
"\n",
"phase = traj.add_phase('phase0',\n",
" dm.Phase(ode_class=AircraftODE,\n",
" transcription=dm.Radau(num_segments=num_seg,\n",
" segment_ends=seg_ends,\n",
" order=3, compressed=False)))\n",
"\n",
"# Pass Reference Area from an external source\n",
"assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())\n",
"assumptions.add_output('S', val=427.8, units='m**2')\n",
"assumptions.add_output('mass_empty', val=1.0, units='kg')\n",
"assumptions.add_output('mass_payload', val=1.0, units='kg')\n",
"\n",
"phase.set_time_options(fix_initial=True,\n",
" duration_bounds=(300, 10000),\n",
" duration_ref=1000)\n",
"\n",
"phase.add_state('range', units='NM',\n",
" rate_source='range_rate_comp.dXdt:range',\n",
" fix_initial=True, fix_final=False, ref=1e-3,\n",
" defect_ref=1e-3, lower=0, upper=2000)\n",
"\n",
"phase.add_state('mass_fuel', units='lbm',\n",
" rate_source='propulsion.dXdt:mass_fuel',\n",
" fix_initial=True, fix_final=True,\n",
" upper=1.5E5, lower=0.0, ref=1e2, defect_ref=1e2)\n",
"\n",
"phase.add_state('alt', units='kft',\n",
" rate_source='climb_rate',\n",
" fix_initial=True, fix_final=True,\n",
" lower=0.0, upper=60, ref=1e-3, defect_ref=1e-3)\n",
"\n",
"phase.add_control('climb_rate', units='ft/min', opt=True, lower=-3000, upper=3000,\n",
" targets=['gam_comp.climb_rate'],\n",
" rate_continuity=True, rate2_continuity=False)\n",
"\n",
"phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], units=None, opt=False)\n",
"\n",
"phase.add_parameter('S',\n",
" targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],\n",
" units='m**2')\n",
"\n",
"phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')\n",
"phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')\n",
"\n",
"phase.add_path_constraint('propulsion.tau', lower=0.01, upper=2.0)\n",
"\n",
"p.model.connect('assumptions.S', 'traj.phase0.parameters:S')\n",
"p.model.connect('assumptions.mass_empty', 'traj.phase0.parameters:mass_empty')\n",
"p.model.connect('assumptions.mass_payload', 'traj.phase0.parameters:mass_payload')\n",
"\n",
"phase.add_objective('range', loc='final', ref=-1.0)\n",
"\n",
"phase.add_timeseries_output('aero.CL')\n",
"phase.add_timeseries_output('aero.CD')\n",
"\n",
"p.setup()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set the initial guesses and assumptions"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p['traj.phase0.t_initial'] = 0.0\n",
"p['traj.phase0.t_duration'] = 3600.0\n",
"p['traj.phase0.states:range'] = phase.interp('range', ys=(0, 724.0))\n",
"p['traj.phase0.states:mass_fuel'] = phase.interp('mass_fuel', ys=(30000, 1e-3))\n",
"p['traj.phase0.states:alt'][:] = 10.0\n",
"\n",
"p['traj.phase0.controls:mach'][:] = 0.8\n",
"\n",
"p['assumptions.S'] = 427.8\n",
"p['assumptions.mass_empty'] = 0.15E6\n",
"p['assumptions.mass_payload'] = 84.02869 * 400"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Run the problem\n",
"\n",
"The dymos `run_problem` method is used to automatically run the driver and follow the driver execution with a simulation run to help verify the results.\n",
"Plots of the resulting solution and simulation will be generated."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dm.run_problem(p, simulate=True, make_plots=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"from IPython.display import IFrame\n",
"IFrame(src=str(p.get_reports_dir() / 'traj_results_report.html'), width=800, height=1200)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"remove-input",
"remove-output"
]
},
"outputs": [],
"source": [
"from openmdao.utils.assert_utils import assert_near_equal\n",
"\n",
"assert_near_equal(p.get_val('traj.phase0.timeseries.range', units='NM')[-1],\n",
" 726.85, tolerance=1.0E-2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"```{bibliography}\n",
":filter: docname in docnames\n",
":labelprefix: C\n",
":keyprefix: ex_commercial_aircraft-\n",
"```"
]
}
],
"metadata": {
"celltoolbar": "Edit Metadata",
"jupytext": {
"cell_metadata_filter": "-all",
"notebook_metadata_filter": "-all",
"text_representation": {
"extension": ".md",
"format_name": "markdown"
}
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}