From 0bbed232b53aa466fd2a245e060bbc7ee2c3eb1c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Thu, 6 Jan 2022 14:05:39 -0500 Subject: [PATCH 01/10] added balanced field example that uses an ExplicitFuncComp as the ODE. --- .../balanced_field_funccomp.ipynb | 590 ++++++++++++++++++ 1 file changed, 590 insertions(+) create mode 100644 docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb diff --git a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb new file mode 100644 index 000000000..f69c3da10 --- /dev/null +++ b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb @@ -0,0 +1,590 @@ +{ + "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": [ + "# Balanced Field Length using OpenMDAO's Function Wrapping Components.\n", + "\n", + "```{admonition} Things you'll learn through this example\n", + "- How to build an ODE that uses existing functions rather than rewriting the dynamics as OpenMDAO components.\n", + "```\n", + "\n", + "This example performs the same balanced field length optimization as the example found [here](examples:balanced_field).\n", + "The purpose of this example is to demonstrate the ability of OpenMDAO to wrap existing user-defined functions that was introduced with OpenMDAO 3.14.0.\n", + "\n", + "## Why use the function wrapping capability?\n", + "\n", + "In many cases, users are likely to have existing code that performs their calculations.\n", + "Rather than forcing them to rewrite this code in the form of OpenMDAO components, the function wrapping components make it easier to get models up and running in OpenMDAO.\n", + "\n", + "Also, doing so gives users another option for differenting their models by providing a format that can be used by Google's [jax package](https://github.com/google/jax), which is a python-based automatic differentiation (AD) tool that provides another option, in addition to complex-step, for providing derivatives across the user's model.\n", + "\n", + "In our case, we have two functions to wrap: the runway ODE and the climb ODE.\n", + "\n", + "### Runway ODE function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def runway_ode(rho, S, CD0, CL0, CL_max, alpha_max, h_w, AR, e, span, T, mu_r, m, v, h, alpha):\n", + " g = 9.80665\n", + " W = m * g\n", + " v_stall = np.sqrt(2 * W / rho / S / CL_max)\n", + " v_over_v_stall = v / v_stall\n", + "\n", + " CL = CL0 + (alpha / alpha_max) * (CL_max - CL0)\n", + " K_nom = 1.0 / (np.pi * AR * e)\n", + " b = span / 2.0\n", + " fact = ((h + h_w) / b) ** 1.5\n", + " K = K_nom * 33 * fact / (1.0 + 33 * fact)\n", + "\n", + " q = 0.5 * rho * v ** 2\n", + " L = q * S * CL\n", + " D = q * S * (CD0 + K * CL ** 2)\n", + "\n", + " # Compute the downward force on the landing gear\n", + " calpha = np.cos(alpha)\n", + " salpha = np.sin(alpha)\n", + "\n", + " # Runway normal force\n", + " F_r = m * g - L * calpha - T * salpha\n", + "\n", + " # Compute the dynamics\n", + " v_dot = (T * calpha - D - F_r * mu_r) / m\n", + " r_dot = v\n", + "\n", + " return CL, q, L, D, K, F_r, v_dot, r_dot, W, v_stall, v_over_v_stall" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Climb ODE function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def climb_ode(rho, S, CD0, CL0, CL_max, alpha_max, h_w, AR, e, span, T, m, v, h, alpha, gam):\n", + " g = 9.80665\n", + " W = m * g\n", + " v_stall = np.sqrt(2 * W / rho / S / CL_max)\n", + " v_over_v_stall = v / v_stall\n", + "\n", + " CL = CL0 + (alpha / alpha_max) * (CL_max - CL0)\n", + " K_nom = 1.0 / (np.pi * AR * e)\n", + " b = span / 2.0\n", + " fact = ((h + h_w) / b) ** 1.5\n", + " K = K_nom * 33 * fact / (1.0 + 33 * fact)\n", + "\n", + " q = 0.5 * rho * v ** 2\n", + " L = q * S * CL\n", + " D = q * S * (CD0 + K * CL ** 2)\n", + "\n", + " # Compute the downward force on the landing gear\n", + " calpha = np.cos(alpha)\n", + " salpha = np.sin(alpha)\n", + "\n", + " # Runway normal force\n", + " F_r = m * g - L * calpha - T * salpha\n", + "\n", + " # Compute the dynamics\n", + " cgam = np.cos(gam)\n", + " sgam = np.sin(gam)\n", + " v_dot = (T * calpha - D) / m - g * sgam\n", + " gam_dot = (T * salpha + L) / (m * v) - (g / v) * cgam\n", + " h_dot = v * sgam\n", + " r_dot = v * cgam\n", + "\n", + " return CL, q, L, D, K, F_r, v_dot, r_dot, W, v_stall, v_over_v_stall, gam_dot, h_dot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The Function Wrapping Component\n", + "\n", + "The following component wraps one of the above ODE functions depending on its 'mode'.\n", + "This allows us to declare a single ODE component for use in the problem regardless of whether the given phase models the ground roll or the climb." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import openmdao.api as om\n", + "import openmdao.func_api as omf\n", + "\n", + "\n", + "class BalancedFieldODEComp(om.ExplicitFuncComp):\n", + " \n", + " def __init__(self, **kwargs):\n", + " nn = kwargs['num_nodes'] \n", + " mode = kwargs.get('mode', 'runway')\n", + " ode_func = runway_ode if mode == 'runway' else climb_ode\n", + " \n", + " meta = (omf.wrap(ode_func)\n", + " .add_input('rho', val=1.225, desc='atmospheric density at runway', units='kg/m**3')\n", + " .add_input('S', val=124.7, desc='aerodynamic reference area', units='m**2')\n", + " .add_input('CD0', val=0.03, desc='zero-lift drag coefficient', units=None)\n", + " .add_input('CL0', val=0.5, desc='zero-alpha lift coefficient', units=None)\n", + " .add_input('CL_max', val=2.0, desc='maximum lift coefficient for linear fit', units=None)\n", + " .add_input('alpha_max', val=np.radians(10), desc='angle of attack at CL_max', units='rad')\n", + " .add_input('h_w', val=1.0, desc='height of the wing above the CG', units='m')\n", + " .add_input('AR', val=9.45, desc='wing aspect ratio', units=None)\n", + " .add_input('e', val=0.801, desc='Oswald span efficiency factor', units=None)\n", + " .add_input('span', val=35.7, desc='Wingspan', units='m')\n", + " .add_input('T', val=1.0, desc='thrust', units='N')\n", + "\n", + " # Dynamic inputs (can assume a different value at every node)\n", + " .add_input('m', shape=(nn,), desc='aircraft mass', units='kg')\n", + " .add_input('v', shape=(nn,), desc='aircraft true airspeed', units='m/s')\n", + " .add_input('h', shape=(nn,), desc='altitude', units='m')\n", + " .add_input('alpha', shape=(nn,), desc='angle of attack', units='rad')\n", + " \n", + " # Outputs\n", + " .add_output('CL', shape=(nn,), desc='lift coefficient', units=None)\n", + " .add_output('q', shape=(nn,), desc='dynamic pressure', units='Pa')\n", + " .add_output('L', shape=(nn,), desc='lift force', units='N')\n", + " .add_output('D', shape=(nn,), desc='drag force', units='N')\n", + " .add_output('K', val=np.ones(nn), desc='drag-due-to-lift factor', units=None)\n", + " .add_output('F_r', shape=(nn,), desc='runway normal force', units='N')\n", + " .add_output('v_dot', shape=(nn,), desc='rate of change of speed', units='m/s**2',\n", + " tags=['dymos.state_rate_source:v'])\n", + " .add_output('r_dot', shape=(nn,), desc='rate of change of range', units='m/s',\n", + " tags=['dymos.state_rate_source:r'])\n", + " .add_output('W', shape=(nn,), desc='aircraft weight', units='N')\n", + " .add_output('v_stall', shape=(nn,), desc='stall speed', units='m/s')\n", + " .add_output('v_over_v_stall', shape=(nn,), desc='stall speed ratio', units=None)\n", + " )\n", + " \n", + " if mode == 'runway':\n", + " meta.add_input('mu_r', val=0.05, desc='runway friction coefficient', units=None)\n", + " else:\n", + " meta.add_input('gam', shape=(nn,), desc='flight path angle', units='rad')\n", + " meta.add_output('gam_dot', shape=(nn,), desc='rate of change of flight path angle',\n", + " units='rad/s', tags=['dymos.state_rate_source:gam'])\n", + " meta.add_output('h_dot', shape=(nn,), desc='rate of change of altitude', units='m/s',\n", + " tags=['dymos.state_rate_source:h'])\n", + "\n", + " grad_method = 'jax'\n", + " meta.declare_coloring('*', method=grad_method)\n", + " meta.declare_partials(of='*', wrt='*', method=grad_method)\n", + " \n", + " # Remove kwargs that aren't options on the component\n", + " kwargs.pop('mode')\n", + " \n", + " super().__init__(meta, **kwargs)\n", + " \n", + " def initialize(self):\n", + " self.options.declare('num_nodes', types=int)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Building and running the problem\n", + "\n", + "This code is identical to that in the previous [balanced field example](examples:balanced_field)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "output_scroll" + ] + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import openmdao.api as om\n", + "from openmdao.utils.general_utils import set_pyoptsparse_opt\n", + "import dymos as dm\n", + "\n", + "p = om.Problem()\n", + "\n", + "_, optimizer = set_pyoptsparse_opt('IPOPT', fallback=True)\n", + "\n", + "p.driver = om.pyOptSparseDriver()\n", + "p.driver.declare_coloring()\n", + "\n", + "# Use IPOPT if available, with fallback to SLSQP\n", + "p.driver.options['optimizer'] = optimizer\n", + "p.driver.options['print_results'] = False\n", + "if optimizer == 'IPOPT':\n", + " p.driver.opt_settings['print_level'] = 5\n", + " p.driver.opt_settings['mu_strategy'] = 'adaptive'\n", + " p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n", + " p.driver.opt_settings['mu_init'] = 0.01\n", + " p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'\n", + "\n", + "# First Phase: Brake release to V1 - both engines operable\n", + "br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'mode': 'runway'})\n", + "br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0)\n", + "br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0)\n", + "br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0)\n", + "br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg')\n", + "br_to_v1.add_timeseries_output('*')\n", + "\n", + "# Second Phase: Rejected takeoff at V1 - no engines operable\n", + "rto = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'mode': 'runway'})\n", + "rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0)\n", + "rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", + "rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", + "rto.add_parameter('alpha', val=0.0, opt=False, units='deg')\n", + "rto.add_timeseries_output('*')\n", + "\n", + "# Third Phase: V1 to Vr - single engine operable\n", + "v1_to_vr = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'mode': 'runway'})\n", + "v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0)\n", + "v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", + "v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", + "v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg')\n", + "v1_to_vr.add_timeseries_output('*')\n", + "\n", + "# Fourth Phase: Rotate - single engine operable\n", + "rotate = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'mode': 'runway'})\n", + "rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0)\n", + "rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", + "rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", + "rotate.add_polynomial_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, val=[0, 10])\n", + "rotate.add_timeseries_output('*')\n", + "\n", + "# Fifth Phase: Climb to target speed and altitude at end of runway.\n", + "climb = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=5),\n", + " ode_init_kwargs={'mode': 'climb'})\n", + "climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0)\n", + "climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", + "climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0)\n", + "climb.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", + "climb.add_state('gam', fix_initial=True, lower=0, ref=0.05, defect_ref=0.05)\n", + "climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10)\n", + "climb.add_timeseries_output('*')\n", + "\n", + "# Instantiate the trajectory and add phases\n", + "traj = dm.Trajectory()\n", + "p.model.add_subsystem('traj', traj)\n", + "traj.add_phase('br_to_v1', br_to_v1)\n", + "traj.add_phase('rto', rto)\n", + "traj.add_phase('v1_to_vr', v1_to_vr)\n", + "traj.add_phase('rotate', rotate)\n", + "traj.add_phase('climb', climb)\n", + "\n", + "# Add parameters common to multiple phases to the trajectory\n", + "traj.add_parameter('m', val=174200., opt=False, units='lbm',\n", + " desc='aircraft mass',\n", + " targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'],\n", + " 'rotate': ['m'], 'climb': ['m']})\n", + "\n", + "traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True,\n", + " desc='nominal aircraft thrust',\n", + " targets={'br_to_v1': ['T']})\n", + "\n", + "traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True,\n", + " desc='thrust under a single engine',\n", + " targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']})\n", + "\n", + "traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True,\n", + " desc='thrust when engines are shut down for rejected takeoff',\n", + " targets={'rto': ['T']})\n", + "\n", + "traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True,\n", + " desc='nominal runway friction coefficient',\n", + " targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']})\n", + "\n", + "traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True,\n", + " desc='runway friction coefficient under braking',\n", + " targets={'rto': ['mu_r']})\n", + "\n", + "traj.add_parameter('h_runway', val=0., opt=False, units='ft',\n", + " desc='runway altitude',\n", + " targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'],\n", + " 'rotate': ['h']})\n", + "\n", + "traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True,\n", + " desc='atmospheric density',\n", + " targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'],\n", + " 'rotate': ['rho']})\n", + "\n", + "traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True,\n", + " desc='aerodynamic reference area',\n", + " targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'],\n", + " 'rotate': ['S'], 'climb': ['S']})\n", + "\n", + "traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True,\n", + " desc='zero-lift drag coefficient',\n", + " targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True,\n", + " desc='wing aspect ratio',\n", + " targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('e', val=801, opt=False, units=None, static_target=True,\n", + " desc='Oswald span efficiency factor',\n", + " targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True,\n", + " desc='wingspan',\n", + " targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True,\n", + " desc='height of wing above CG',\n", + " targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True,\n", + " desc='zero-alpha lift coefficient',\n", + " targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True,\n", + " desc='maximum lift coefficient for linear fit',\n", + " targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "traj.add_parameter('alpha_max', val=10.0, opt=False, units='deg', static_target=True,\n", + " desc='angle of attack at maximum lift',\n", + " targets={f'{phase}': ['alpha_max'] for phase in ['br_to_v1', 'v1_to_vr',\n", + " 'rto', 'rotate' 'climb']})\n", + "\n", + "# Standard \"end of first phase to beginning of second phase\" linkages\n", + "# Alpha changes from being a parameter in v1_to_vr to a polynomial control\n", + "# in rotate, to a dynamic control in `climb`.\n", + "traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v'])\n", + "traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha'])\n", + "traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha'])\n", + "traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v'])\n", + "\n", + "# Less common \"final value of r must match at ends of two phases\".\n", + "traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final',\n", + " phase_b='climb', var_b='r', loc_b='final',\n", + " ref=1000)\n", + "\n", + "# Define the constraints and objective for the optimal control problem\n", + "v1_to_vr.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=100)\n", + "\n", + "rto.add_boundary_constraint('v', loc='final', equals=0., ref=100, linear=True)\n", + "\n", + "rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000)\n", + "\n", + "climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True)\n", + "climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True)\n", + "climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg')\n", + "climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.25, ref=1.25)\n", + "\n", + "rto.add_objective('r', loc='final', ref=1000.0)\n", + "\n", + "#\n", + "# Setup the problem and set the initial guess\n", + "#\n", + "p.setup(check=True)\n", + "\n", + "p.set_val('traj.br_to_v1.t_initial', 0)\n", + "p.set_val('traj.br_to_v1.t_duration', 35)\n", + "p.set_val('traj.br_to_v1.states:r', br_to_v1.interp('r', [0, 2500.0]))\n", + "p.set_val('traj.br_to_v1.states:v', br_to_v1.interp('v', [0, 100.0]))\n", + "p.set_val('traj.br_to_v1.parameters:alpha', 0, units='deg')\n", + "\n", + "p.set_val('traj.v1_to_vr.t_initial', 35)\n", + "p.set_val('traj.v1_to_vr.t_duration', 35)\n", + "p.set_val('traj.v1_to_vr.states:r', v1_to_vr.interp('r', [2500, 300.0]))\n", + "p.set_val('traj.v1_to_vr.states:v', v1_to_vr.interp('v', [100, 110.0]))\n", + "p.set_val('traj.v1_to_vr.parameters:alpha', 0.0, units='deg')\n", + "\n", + "p.set_val('traj.rto.t_initial', 35)\n", + "p.set_val('traj.rto.t_duration', 35)\n", + "p.set_val('traj.rto.states:r', rto.interp('r', [2500, 5000.0]))\n", + "p.set_val('traj.rto.states:v', rto.interp('v', [110, 0]))\n", + "p.set_val('traj.rto.parameters:alpha', 0.0, units='deg')\n", + "\n", + "p.set_val('traj.rotate.t_initial', 70)\n", + "p.set_val('traj.rotate.t_duration', 5)\n", + "p.set_val('traj.rotate.states:r', rotate.interp('r', [1750, 1800.0]))\n", + "p.set_val('traj.rotate.states:v', rotate.interp('v', [80, 85.0]))\n", + "p.set_val('traj.rotate.polynomial_controls:alpha', 0.0, units='deg')\n", + "\n", + "p.set_val('traj.climb.t_initial', 75)\n", + "p.set_val('traj.climb.t_duration', 15)\n", + "p.set_val('traj.climb.states:r', climb.interp('r', [5000, 5500.0]), units='ft')\n", + "p.set_val('traj.climb.states:v', climb.interp('v', [160, 170.0]), units='kn')\n", + "p.set_val('traj.climb.states:h', climb.interp('h', [0, 35.0]), units='ft')\n", + "p.set_val('traj.climb.states:gam', climb.interp('gam', [0, 5.0]), units='deg')\n", + "p.set_val('traj.climb.controls:alpha', 5.0, units='deg')\n", + "\n", + "dm.run_problem(p, run_driver=True, simulate=True)\n", + "\n", + "print(p.get_val('traj.rto.states:r')[-1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim_case = om.CaseReader('dymos_solution.db').get_case('final')\n", + "\n", + "fig, axes = plt.subplots(2, 1, sharex=True, gridspec_kw={'top': 0.92}, figsize=(12, 6))\n", + "for phase in ['br_to_v1', 'rto', 'v1_to_vr', 'rotate', 'climb']:\n", + " r = sim_case.get_val(f'traj.{phase}.timeseries.states:r', units='ft')\n", + " v = sim_case.get_val(f'traj.{phase}.timeseries.states:v', units='kn')\n", + " t = sim_case.get_val(f'traj.{phase}.timeseries.time', units='s')\n", + " axes[0].plot(t, r, '-', label=phase)\n", + " axes[1].plot(t, v, '-', label=phase)\n", + "fig.suptitle('Balanced Field Length')\n", + "axes[1].set_xlabel('time (s)')\n", + "axes[0].set_ylabel('range (ft)')\n", + "axes[1].set_ylabel('airspeed (kts)')\n", + "axes[0].grid(True)\n", + "axes[1].grid(True)\n", + "\n", + "tv1 = sim_case.get_val('traj.br_to_v1.timeseries.time', units='s')[-1, 0]\n", + "v1 = sim_case.get_val('traj.br_to_v1.timeseries.states:v', units='kn')[-1, 0]\n", + "\n", + "tf_rto = sim_case.get_val('traj.rto.timeseries.time', units='s')[-1, 0]\n", + "rf_rto = sim_case.get_val('traj.rto.timeseries.states:r', units='ft')[-1, 0]\n", + "\n", + "axes[0].annotate(f'field length = {r[-1, 0]:5.1f} ft', xy=(t[-1, 0], r[-1, 0]),\n", + " xycoords='data', xytext=(0.7, 0.5),\n", + " textcoords='axes fraction', arrowprops=dict(arrowstyle='->'),\n", + " horizontalalignment='center', verticalalignment='top')\n", + "\n", + "axes[0].annotate(f'', xy=(tf_rto, rf_rto),\n", + " xycoords='data', xytext=(0.7, 0.5),\n", + " textcoords='axes fraction', arrowprops=dict(arrowstyle='->'),\n", + " horizontalalignment='center', verticalalignment='top')\n", + "\n", + "axes[1].annotate(f'$v1$ = {v1:5.1f} kts', xy=(tv1, v1), xycoords='data', xytext=(0.5, 0.5),\n", + " textcoords='axes fraction', arrowprops=dict(arrowstyle='->'),\n", + " horizontalalignment='center', verticalalignment='top')\n", + "\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "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.rto.states:r')[-1], 2188.2, tolerance=0.01)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "```{bibliography}\n", + ":filter: docname in docnames\n", + "```" + ] + } + ], + "metadata": { + "celltoolbar": "Tags", + "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.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 5f18f507918ce4d3dbfa79e566ca225ae63d93f1 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Thu, 6 Jan 2022 14:06:21 -0500 Subject: [PATCH 02/10] added a cross reference tag to the balanced field example, needs to be merged with upstream --- docs/dymos_book/_toc.yml | 1 + .../examples/balanced_field/balanced_field.ipynb | 7 +++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml index 48868ec2b..6ca4d873b 100644 --- a/docs/dymos_book/_toc.yml +++ b/docs/dymos_book/_toc.yml @@ -20,6 +20,7 @@ parts: - file: examples/brachistochrone/brachistochrone_upstream_init_duration_states - file: examples/vanderpol/vanderpol - file: examples/balanced_field/balanced_field + - file: examples/balanced_field/balanced_field_funccomp - file: examples/bryson_denham/bryson_denham - file: examples/commercial_aircraft/commercial_aircraft - file: examples/double_integrator/double_integrator diff --git a/docs/dymos_book/examples/balanced_field/balanced_field.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field.ipynb index 72b45ac20..a4d073beb 100644 --- a/docs/dymos_book/examples/balanced_field/balanced_field.ipynb +++ b/docs/dymos_book/examples/balanced_field/balanced_field.ipynb @@ -47,6 +47,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "\n", + "(examples:balanced_field)=\n", + "\n", "# Aircraft Balanced Field Length Calculation\n", "\n", "```{admonition} Things you'll learn through this example\n", @@ -733,7 +736,7 @@ } }, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -747,7 +750,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.1" + "version": "3.9.7" } }, "nbformat": 4, From cc6e1db45a1484439b8d94f9bc6562368ad5dd0a Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Sat, 8 Jan 2022 11:19:23 -0500 Subject: [PATCH 03/10] added jax and jaxlib as dependencies for [docs] --- setup.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 6ad5af8b7..2393d6dd7 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,9 @@ 'ipython', 'numpydoc>=0.9.1', 'redbaron', - 'tabulate' + 'tabulate', + 'jaxlib', + 'jax' ], 'notebooks': [ 'notebook', From 2164a7b3c07c78ffa2bdf599b8715f4a823ece7e Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Fri, 28 Jan 2022 12:53:07 -0500 Subject: [PATCH 04/10] Addressed Justin Gray's review. Added a factory function to return a copy of the ODE class. This function behaves like a class constructor, so from the dymos perspective it can be used in place of ode_class. Found an issue where Oswald's efficiency factor was 801 instead of 0.801. Cleaned up trajectory parameter settings, leaving them default where possible for brevity. --- .../balanced_field/balanced_field.ipynb | 72 ++---- .../balanced_field_funccomp.ipynb | 218 +++++++----------- 2 files changed, 102 insertions(+), 188 deletions(-) diff --git a/docs/dymos_book/examples/balanced_field/balanced_field.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field.ipynb index 714415044..a23de40c6 100644 --- a/docs/dymos_book/examples/balanced_field/balanced_field.ipynb +++ b/docs/dymos_book/examples/balanced_field/balanced_field.ipynb @@ -228,7 +228,7 @@ "|S |aerodynamic reference area |$m^2$ | 124.7 |\n", "|CD0 |zero-lift drag coefficient |- | 0.03 |\n", "|AR |wing aspect ratio |- | 9.45 |\n", - "|e |Oswald's wing efficiency |- | 801 |\n", + "|e |Oswald's wing efficiency |- | 0.801 |\n", "|span |wingspan |$m$ | 35.7 |\n", "|h_w |height of wing above CoG |$m$ | 1.0 |\n", "|CL0 |aerodynamic reference area |- | 0.5 |\n", @@ -446,7 +446,10 @@ "p.driver.options['print_results'] = False\n", "if optimizer == 'IPOPT':\n", " p.driver.opt_settings['print_level'] = 5\n", - " p.driver.opt_settings['derivative_test'] = 'first-order'" + " p.driver.opt_settings['mu_strategy'] = 'adaptive'\n", + " p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n", + " p.driver.opt_settings['mu_init'] = 0.01\n", + " p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'" ] }, { @@ -523,7 +526,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the next part of our script we define the trajectory-level parameters. We assign these to their corresponding targets in the various phases of the trajectory." + "In the next part of our script we define the trajectory-level parameters. We assign these to their corresponding targets in the various phases of the trajectory.\n", + "\n", + "We're omitting a lot of the scalar inputs to the ODE and letting them assume their default values for the sake of brevity." ] }, { @@ -532,13 +537,15 @@ "metadata": {}, "outputs": [], "source": [ + "all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']\n", + "groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate']\n", "\n", "# Add parameters common to multiple phases to the trajectory\n", "traj.add_parameter('m', val=174200., opt=False, units='lbm',\n", " desc='aircraft mass',\n", - " targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'],\n", - " 'rotate': ['m'], 'climb': ['m']})\n", + " targets={phase: ['m'] for phase in all_phases})\n", "\n", + "# Handle parameters which change from phase to phase.\n", "traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True,\n", " desc='nominal aircraft thrust',\n", " targets={'br_to_v1': ['T']})\n", @@ -561,58 +568,7 @@ "\n", "traj.add_parameter('h_runway', val=0., opt=False, units='ft',\n", " desc='runway altitude',\n", - " targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'],\n", - " 'rotate': ['h']})\n", - "\n", - "traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True,\n", - " desc='atmospheric density',\n", - " targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'],\n", - " 'rotate': ['rho']})\n", - "\n", - "traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True,\n", - " desc='aerodynamic reference area',\n", - " targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'],\n", - " 'rotate': ['S'], 'climb': ['S']})\n", - "\n", - "traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True,\n", - " desc='zero-lift drag coefficient',\n", - " targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True,\n", - " desc='wing aspect ratio',\n", - " targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('e', val=801, opt=False, units=None, static_target=True,\n", - " desc='Oswald span efficiency factor',\n", - " targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True,\n", - " desc='wingspan',\n", - " targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True,\n", - " desc='height of wing above CG',\n", - " targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True,\n", - " desc='zero-alpha lift coefficient',\n", - " targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True,\n", - " desc='maximum lift coefficient for linear fit',\n", - " targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('alpha_max', val=10.0, opt=False, units='deg', static_target=True,\n", - " desc='angle of attack at maximum lift',\n", - " targets={f'{phase}': ['alpha_max'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})" + " targets={phase: ['h'] for phase in groundroll_phases})" ] }, { @@ -807,7 +763,7 @@ "source": [ "from openmdao.utils.assert_utils import assert_near_equal\n", "\n", - "assert_near_equal(p.get_val('traj.rto.states:r')[-1], 2188.2, tolerance=0.01)" + "assert_near_equal(p.get_val('traj.rto.states:r')[-1], 2197.7, tolerance=0.01)" ] }, { diff --git a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb index f69c3da10..3e97c5536 100644 --- a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb +++ b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb @@ -61,7 +61,7 @@ "In many cases, users are likely to have existing code that performs their calculations.\n", "Rather than forcing them to rewrite this code in the form of OpenMDAO components, the function wrapping components make it easier to get models up and running in OpenMDAO.\n", "\n", - "Also, doing so gives users another option for differenting their models by providing a format that can be used by Google's [jax package](https://github.com/google/jax), which is a python-based automatic differentiation (AD) tool that provides another option, in addition to complex-step, for providing derivatives across the user's model.\n", + "Also, doing so gives users another option for differentiating their models by providing a format that can be used by Google's [jax package](https://github.com/google/jax), which is a python-based automatic differentiation (AD) tool that provides another option, in addition to complex-step, for providing derivatives across the user's model.\n", "\n", "In our case, we have two functions to wrap: the runway ODE and the climb ODE.\n", "\n", @@ -155,10 +155,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### The Function Wrapping Component\n", + "### A component creating function\n", "\n", - "The following component wraps one of the above ODE functions depending on its 'mode'.\n", - "This allows us to declare a single ODE component for use in the problem regardless of whether the given phase models the ground roll or the climb." + "This function accepts num_nodes and returns a component that wraps the above functions." ] }, { @@ -170,69 +169,73 @@ "import openmdao.api as om\n", "import openmdao.func_api as omf\n", "\n", - "\n", - "class BalancedFieldODEComp(om.ExplicitFuncComp):\n", + "def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=False):\n", + " \"\"\"\n", + " Returns the metadata from omf needed to create a new ExplciitFuncComp.\n", + " \"\"\"\n", + " nn = num_nodes\n", + " ode_func = runway_ode if flight_mode == 'runway' else climb_ode\n", + " \n", + " meta = (omf.wrap(ode_func)\n", + " .add_input('rho', val=1.225, desc='atmospheric density at runway', units='kg/m**3')\n", + " .add_input('S', val=124.7, desc='aerodynamic reference area', units='m**2')\n", + " .add_input('CD0', val=0.03, desc='zero-lift drag coefficient', units=None)\n", + " .add_input('CL0', val=0.5, desc='zero-alpha lift coefficient', units=None)\n", + " .add_input('CL_max', val=2.0, desc='maximum lift coefficient for linear fit', units=None)\n", + " .add_input('alpha_max', val=np.radians(10), desc='angle of attack at CL_max', units='rad')\n", + " .add_input('h_w', val=1.0, desc='height of the wing above the CG', units='m')\n", + " .add_input('AR', val=9.45, desc='wing aspect ratio', units=None)\n", + " .add_input('e', val=0.801, desc='Oswald span efficiency factor', units=None)\n", + " .add_input('span', val=35.7, desc='Wingspan', units='m')\n", + " .add_input('T', val=1.0, desc='thrust', units='N')\n", + "\n", + " # Dynamic inputs (can assume a different value at every node)\n", + " .add_input('m', shape=(nn,), desc='aircraft mass', units='kg')\n", + " .add_input('v', shape=(nn,), desc='aircraft true airspeed', units='m/s')\n", + " .add_input('h', shape=(nn,), desc='altitude', units='m')\n", + " .add_input('alpha', shape=(nn,), desc='angle of attack', units='rad')\n", + "\n", + " # Outputs\n", + " .add_output('CL', shape=(nn,), desc='lift coefficient', units=None)\n", + " .add_output('q', shape=(nn,), desc='dynamic pressure', units='Pa')\n", + " .add_output('L', shape=(nn,), desc='lift force', units='N')\n", + " .add_output('D', shape=(nn,), desc='drag force', units='N')\n", + " .add_output('K', val=np.ones(nn), desc='drag-due-to-lift factor', units=None)\n", + " .add_output('F_r', shape=(nn,), desc='runway normal force', units='N')\n", + " .add_output('v_dot', shape=(nn,), desc='rate of change of speed', units='m/s**2',\n", + " tags=['dymos.state_rate_source:v'])\n", + " .add_output('r_dot', shape=(nn,), desc='rate of change of range', units='m/s',\n", + " tags=['dymos.state_rate_source:r'])\n", + " .add_output('W', shape=(nn,), desc='aircraft weight', units='N')\n", + " .add_output('v_stall', shape=(nn,), desc='stall speed', units='m/s')\n", + " .add_output('v_over_v_stall', shape=(nn,), desc='stall speed ratio', units=None)\n", + " )\n", + "\n", + " if flight_mode == 'runway':\n", + " meta.add_input('mu_r', val=0.05, desc='runway friction coefficient', units=None)\n", + " else:\n", + " meta.add_input('gam', shape=(nn,), desc='flight path angle', units='rad')\n", + " meta.add_output('gam_dot', shape=(nn,), desc='rate of change of flight path angle',\n", + " units='rad/s', tags=['dymos.state_rate_source:gam'])\n", + " meta.add_output('h_dot', shape=(nn,), desc='rate of change of altitude', units='m/s',\n", + " tags=['dymos.state_rate_source:h'])\n", + " \n", + " meta.declare_coloring('*', method=grad_method)\n", + " meta.declare_partials(of='*', wrt='*', method=grad_method)\n", " \n", - " def __init__(self, **kwargs):\n", - " nn = kwargs['num_nodes'] \n", - " mode = kwargs.get('mode', 'runway')\n", - " ode_func = runway_ode if mode == 'runway' else climb_ode\n", - " \n", - " meta = (omf.wrap(ode_func)\n", - " .add_input('rho', val=1.225, desc='atmospheric density at runway', units='kg/m**3')\n", - " .add_input('S', val=124.7, desc='aerodynamic reference area', units='m**2')\n", - " .add_input('CD0', val=0.03, desc='zero-lift drag coefficient', units=None)\n", - " .add_input('CL0', val=0.5, desc='zero-alpha lift coefficient', units=None)\n", - " .add_input('CL_max', val=2.0, desc='maximum lift coefficient for linear fit', units=None)\n", - " .add_input('alpha_max', val=np.radians(10), desc='angle of attack at CL_max', units='rad')\n", - " .add_input('h_w', val=1.0, desc='height of the wing above the CG', units='m')\n", - " .add_input('AR', val=9.45, desc='wing aspect ratio', units=None)\n", - " .add_input('e', val=0.801, desc='Oswald span efficiency factor', units=None)\n", - " .add_input('span', val=35.7, desc='Wingspan', units='m')\n", - " .add_input('T', val=1.0, desc='thrust', units='N')\n", - "\n", - " # Dynamic inputs (can assume a different value at every node)\n", - " .add_input('m', shape=(nn,), desc='aircraft mass', units='kg')\n", - " .add_input('v', shape=(nn,), desc='aircraft true airspeed', units='m/s')\n", - " .add_input('h', shape=(nn,), desc='altitude', units='m')\n", - " .add_input('alpha', shape=(nn,), desc='angle of attack', units='rad')\n", - " \n", - " # Outputs\n", - " .add_output('CL', shape=(nn,), desc='lift coefficient', units=None)\n", - " .add_output('q', shape=(nn,), desc='dynamic pressure', units='Pa')\n", - " .add_output('L', shape=(nn,), desc='lift force', units='N')\n", - " .add_output('D', shape=(nn,), desc='drag force', units='N')\n", - " .add_output('K', val=np.ones(nn), desc='drag-due-to-lift factor', units=None)\n", - " .add_output('F_r', shape=(nn,), desc='runway normal force', units='N')\n", - " .add_output('v_dot', shape=(nn,), desc='rate of change of speed', units='m/s**2',\n", - " tags=['dymos.state_rate_source:v'])\n", - " .add_output('r_dot', shape=(nn,), desc='rate of change of range', units='m/s',\n", - " tags=['dymos.state_rate_source:r'])\n", - " .add_output('W', shape=(nn,), desc='aircraft weight', units='N')\n", - " .add_output('v_stall', shape=(nn,), desc='stall speed', units='m/s')\n", - " .add_output('v_over_v_stall', shape=(nn,), desc='stall speed ratio', units=None)\n", - " )\n", - " \n", - " if mode == 'runway':\n", - " meta.add_input('mu_r', val=0.05, desc='runway friction coefficient', units=None)\n", - " else:\n", - " meta.add_input('gam', shape=(nn,), desc='flight path angle', units='rad')\n", - " meta.add_output('gam_dot', shape=(nn,), desc='rate of change of flight path angle',\n", - " units='rad/s', tags=['dymos.state_rate_source:gam'])\n", - " meta.add_output('h_dot', shape=(nn,), desc='rate of change of altitude', units='m/s',\n", - " tags=['dymos.state_rate_source:h'])\n", - "\n", - " grad_method = 'jax'\n", - " meta.declare_coloring('*', method=grad_method)\n", - " meta.declare_partials(of='*', wrt='*', method=grad_method)\n", - " \n", - " # Remove kwargs that aren't options on the component\n", - " kwargs.pop('mode')\n", - " \n", - " super().__init__(meta, **kwargs)\n", + " return om.ExplicitFuncComp(meta)\n", " \n", - " def initialize(self):\n", - " self.options.declare('num_nodes', types=int)" + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The Function Wrapping Component\n", + "\n", + "The following component wraps one of the above ODE functions depending on its 'mode'.\n", + "This allows us to declare a single ODE component for use in the problem regardless of whether the given phase models the ground roll or the climb." ] }, { @@ -278,8 +281,8 @@ " p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'\n", "\n", "# First Phase: Brake release to V1 - both engines operable\n", - "br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", - " ode_init_kwargs={'mode': 'runway'})\n", + "br_to_v1 = dm.Phase(ode_class=wrap_ode_func, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'flight_mode': 'runway'})\n", "br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0)\n", "br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0)\n", "br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0)\n", @@ -287,8 +290,8 @@ "br_to_v1.add_timeseries_output('*')\n", "\n", "# Second Phase: Rejected takeoff at V1 - no engines operable\n", - "rto = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", - " ode_init_kwargs={'mode': 'runway'})\n", + "rto = dm.Phase(ode_class=wrap_ode_func, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'flight_mode': 'runway'})\n", "rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0)\n", "rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", "rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", @@ -296,8 +299,8 @@ "rto.add_timeseries_output('*')\n", "\n", "# Third Phase: V1 to Vr - single engine operable\n", - "v1_to_vr = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", - " ode_init_kwargs={'mode': 'runway'})\n", + "v1_to_vr = dm.Phase(ode_class=wrap_ode_func, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'flight_mode': 'runway'})\n", "v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0)\n", "v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", "v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", @@ -305,8 +308,8 @@ "v1_to_vr.add_timeseries_output('*')\n", "\n", "# Fourth Phase: Rotate - single engine operable\n", - "rotate = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3),\n", - " ode_init_kwargs={'mode': 'runway'})\n", + "rotate = dm.Phase(ode_class=wrap_ode_func, transcription=dm.Radau(num_segments=3),\n", + " ode_init_kwargs={'flight_mode': 'runway'})\n", "rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0)\n", "rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", "rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)\n", @@ -314,8 +317,8 @@ "rotate.add_timeseries_output('*')\n", "\n", "# Fifth Phase: Climb to target speed and altitude at end of runway.\n", - "climb = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=5),\n", - " ode_init_kwargs={'mode': 'climb'})\n", + "climb = dm.Phase(ode_class=wrap_ode_func, transcription=dm.Radau(num_segments=5),\n", + " ode_init_kwargs={'flight_mode': 'climb'})\n", "climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0)\n", "climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)\n", "climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0)\n", @@ -333,12 +336,15 @@ "traj.add_phase('rotate', rotate)\n", "traj.add_phase('climb', climb)\n", "\n", + "all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']\n", + "groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate']\n", + "\n", "# Add parameters common to multiple phases to the trajectory\n", "traj.add_parameter('m', val=174200., opt=False, units='lbm',\n", " desc='aircraft mass',\n", - " targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'],\n", - " 'rotate': ['m'], 'climb': ['m']})\n", + " targets={phase: ['m'] for phase in all_phases})\n", "\n", + "# Handle parameters which change from phase to phase.\n", "traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True,\n", " desc='nominal aircraft thrust',\n", " targets={'br_to_v1': ['T']})\n", @@ -361,58 +367,10 @@ "\n", "traj.add_parameter('h_runway', val=0., opt=False, units='ft',\n", " desc='runway altitude',\n", - " targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'],\n", - " 'rotate': ['h']})\n", - "\n", - "traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True,\n", - " desc='atmospheric density',\n", - " targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'],\n", - " 'rotate': ['rho']})\n", - "\n", - "traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True,\n", - " desc='aerodynamic reference area',\n", - " targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'],\n", - " 'rotate': ['S'], 'climb': ['S']})\n", - "\n", - "traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True,\n", - " desc='zero-lift drag coefficient',\n", - " targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True,\n", - " desc='wing aspect ratio',\n", - " targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('e', val=801, opt=False, units=None, static_target=True,\n", - " desc='Oswald span efficiency factor',\n", - " targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True,\n", - " desc='wingspan',\n", - " targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True,\n", - " desc='height of wing above CG',\n", - " targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True,\n", - " desc='zero-alpha lift coefficient',\n", - " targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True,\n", - " desc='maximum lift coefficient for linear fit',\n", - " targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", - "\n", - "traj.add_parameter('alpha_max', val=10.0, opt=False, units='deg', static_target=True,\n", - " desc='angle of attack at maximum lift',\n", - " targets={f'{phase}': ['alpha_max'] for phase in ['br_to_v1', 'v1_to_vr',\n", - " 'rto', 'rotate' 'climb']})\n", + " targets={phase: ['h'] for phase in groundroll_phases})\n", + "\n", + "# Here we're omitting some constants that are common throughout all phases for the sake of brevity.\n", + "# Their correct defaults are specified in add_input calls to `wrap_ode_func`.\n", "\n", "# Standard \"end of first phase to beginning of second phase\" linkages\n", "# Alpha changes from being a parameter in v1_to_vr to a polynomial control\n", @@ -542,7 +500,7 @@ "source": [ "from openmdao.utils.assert_utils import assert_near_equal\n", "\n", - "assert_near_equal(p.get_val('traj.rto.states:r')[-1], 2188.2, tolerance=0.01)" + "assert_near_equal(p.get_val('traj.rto.states:r')[-1], 2197.7, tolerance=0.01)" ] }, { From 65937a43240b561309f18a18a363960633c284e7 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Fri, 28 Jan 2022 13:00:56 -0500 Subject: [PATCH 05/10] use jax and jit in the ExplicitFuncCompExample, runtime changed from ~5.6 seconds to ~3.6 seconds --- .../examples/balanced_field/balanced_field_funccomp.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb index 3e97c5536..2ad49e115 100644 --- a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb +++ b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb @@ -169,7 +169,7 @@ "import openmdao.api as om\n", "import openmdao.func_api as omf\n", "\n", - "def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=False):\n", + "def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True):\n", " \"\"\"\n", " Returns the metadata from omf needed to create a new ExplciitFuncComp.\n", " \"\"\"\n", @@ -223,7 +223,7 @@ " meta.declare_coloring('*', method=grad_method)\n", " meta.declare_partials(of='*', wrt='*', method=grad_method)\n", " \n", - " return om.ExplicitFuncComp(meta)\n", + " return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit)\n", " \n", " " ] From eed391136a23096e09efabe5c03ac6ba796d05df Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 31 Jan 2022 10:59:46 -0500 Subject: [PATCH 06/10] Display doc build logs in case of failure --- .github/workflows/dymos_tests_workflow.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index b70dd7e1f..a0fb7f212 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.github/workflows/dymos_tests_workflow.yml @@ -258,6 +258,7 @@ jobs: - name: Build docs if: ${{ matrix.DOCS == '1' || matrix.DOCS == '2' }} + id: build_docs shell: bash -l {0} run: | echo "=============================================================" @@ -269,7 +270,7 @@ jobs: python copy_build_artifacts.py - name: Display doc build reports - if: failure() && matrix.DOCS == '1' && steps.build_docs.outcome == 'failure' + if: failure() && { matrix.DOCS == '1' || matrix.DOCS == '2' } && steps.build_docs.outcome == 'failure' run: | for f in /home/runner/work/OpenMDAO/dymos/docs/dymos_book/_build/html/reports/*; do echo "=============================================================" From 539a1f711940e94f89a1d949083a761b2acb7810 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 31 Jan 2022 11:48:55 -0500 Subject: [PATCH 07/10] action fix --- .github/workflows/dymos_tests_workflow.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index a0fb7f212..ae2ffa3b4 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.github/workflows/dymos_tests_workflow.yml @@ -270,7 +270,7 @@ jobs: python copy_build_artifacts.py - name: Display doc build reports - if: failure() && { matrix.DOCS == '1' || matrix.DOCS == '2' } && steps.build_docs.outcome == 'failure' + if: failure() && (matrix.DOCS == '1' || matrix.DOCS == '2') && steps.build_docs.outcome == 'failure' run: | for f in /home/runner/work/OpenMDAO/dymos/docs/dymos_book/_build/html/reports/*; do echo "=============================================================" From 97226fbc5cd8a92d18f84cf5a4700da94a4280a5 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 31 Jan 2022 13:03:47 -0500 Subject: [PATCH 08/10] More action tweaks. Removed 3.6 from test matrix since pyoptsparse no longer supports it. --- .github/workflows/dymos_tests_workflow.yml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index ae2ffa3b4..937f19580 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.github/workflows/dymos_tests_workflow.yml @@ -70,7 +70,7 @@ jobs: DOCS: 0 # oldest supported versions - - PY: 3.6 + - PY: 3.7 NUMPY: 1.16 SCIPY: 1.2 PETSc: 3.10.2 @@ -270,9 +270,10 @@ jobs: python copy_build_artifacts.py - name: Display doc build reports + continue-on-error: True if: failure() && (matrix.DOCS == '1' || matrix.DOCS == '2') && steps.build_docs.outcome == 'failure' run: | - for f in /home/runner/work/OpenMDAO/dymos/docs/dymos_book/_build/html/reports/*; do + for f in /home/runner/work/dymos/dymos/docs/dymos_book/_build/html/reports/*; do echo "=============================================================" echo $f echo "=============================================================" From b500b1e448ca678d6438afeb480c3dcb068cc31c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 31 Jan 2022 13:04:28 -0500 Subject: [PATCH 09/10] bump minimum python to 3.7 --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 2393d6dd7..212669fec 100644 --- a/setup.py +++ b/setup.py @@ -63,7 +63,7 @@ ], license='Apache License', packages=find_packages(), - python_requires=">=3.6", + python_requires=">=3.7", install_requires=[ 'openmdao>=3.13.0', 'numpy>=1.14.1', From e6978e0a45cc8903b3abbb6942bd07268cb4eda6 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 31 Jan 2022 14:27:35 -0500 Subject: [PATCH 10/10] some numerical noise on the mountain car problem --- docs/dymos_book/examples/mountain_car/mountain_car.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/dymos_book/examples/mountain_car/mountain_car.ipynb b/docs/dymos_book/examples/mountain_car/mountain_car.ipynb index 05c3d5546..b450ab0d1 100644 --- a/docs/dymos_book/examples/mountain_car/mountain_car.ipynb +++ b/docs/dymos_book/examples/mountain_car/mountain_car.ipynb @@ -289,7 +289,7 @@ "p.driver.opt_settings['max_iter'] = 500\n", "p.driver.opt_settings['mu_strategy'] = 'adaptive'\n", "p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n", - "p.driver.opt_settings['mu_init'] = 0.01\n", + "p.driver.opt_settings['tol'] = 1.0E-8\n", "p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' # for faster convergence\n", "\n", "p.driver.declare_coloring()" @@ -530,8 +530,8 @@ "source": [ "from openmdao.utils.assert_utils import assert_near_equal\n", "\n", - "assert_near_equal(t[-1, 0], 102.479, tolerance=1.0E-3)\n", - "assert_near_equal(x[-1, 0], 0.5, tolerance=1.0E-3)" + "assert_near_equal(t[-1, 0], 102.479, tolerance=5.0E-3)\n", + "assert_near_equal(x[-1, 0], 0.5, tolerance=5.0E-3)" ] }, {