diff --git a/.github/workflows/dymos_docs_workflow.yml b/.github/workflows/dymos_docs_workflow.yml index 07956a67a..7b89d6eb0 100644 --- a/.github/workflows/dymos_docs_workflow.yml +++ b/.github/workflows/dymos_docs_workflow.yml @@ -31,14 +31,15 @@ jobs: SCIPY: '1.13' PETSc: '3.19' PYOPTSPARSE: 'v2.11.0' - OPENMDAO: 'dev' + OPENMDAO: 'latest' OPTIONAL: '[docs]' JAX: '0.4.28' PUBLISH_DOCS: 1 # make sure the latest versions of things don't break the docs + # sticking with Python 3.12 for now, 3.13 requires NumPy 2.1 which does not work yet with PETSc/pyoptsparse - NAME: latest - PY: 3 + PY: '3.12' NUMPY: 1 SCIPY: 1 PETSc: 3 @@ -268,7 +269,6 @@ jobs: cd $HOME/work/dymos/dymos ghp-import -n -p -f docs/dymos_book/_build/html - - name: Publish docs to openmdao.org if: | github.event_name == 'push' && github.ref == 'refs/heads/master' && diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index 3a1b9cbde..e6b374099 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.github/workflows/dymos_tests_workflow.yml @@ -60,6 +60,12 @@ on: required: false default: false + python313: + type: boolean + description: "Include 'python313' in test matrix" + required: false + default: false + oldest: type: boolean description: "Include 'oldest' in test matrix" @@ -134,8 +140,9 @@ jobs: EXCLUDE: ${{ github.event_name == 'workflow_dispatch' && ! inputs.no_mpi }} # try latest versions + # sticking with Python 3.12 here, 3.13 requires NumPy 2.1 which does not work yet with pyoptsparse - NAME: latest - PY: 3 + PY: '3.12' NUMPY: 1 SCIPY: 1 PETSc: 3.21.0 @@ -146,6 +153,19 @@ jobs: JAX: 'latest' EXCLUDE: ${{ github.event_name == 'workflow_dispatch' && ! inputs.latest }} + # Python 3.13 (requires NumPy 2.1 which does not work yet with pyoptsparse) + - NAME: python313 + PY: '3.13' + NUMPY: 2 + SCIPY: 1 + PETSc: 3 + # PYOPTSPARSE: 'latest' + SNOPT: 7.7 + OPENMDAO: 'dev' + OPTIONAL: '[test]' + JAX: 'latest' + EXCLUDE: ${{ github.event_name == 'workflow_dispatch' && ! inputs.python313 }} + # oldest supported versions - NAME: oldest PY: 3.9 @@ -237,9 +257,9 @@ jobs: if [[ "${{ matrix.OPENMPI }}" && "${{ matrix.MPI4PY }}" ]]; then conda install openmpi=${{ matrix.OPENMPI }} mpi4py=${{ matrix.MPI4PY }} petsc4py=${{ matrix.PETSc }} -q -y elif [[ "${{ matrix.MPI4PY }}" ]]; then - conda install mpi4py=${{ matrix.MPI4PY }} petsc4py=${{ matrix.PETSc }} -q -y + conda install mpich mpi4py=${{ matrix.MPI4PY }} petsc4py=${{ matrix.PETSc }} -q -y else - conda install mpi4py petsc4py=${{ matrix.PETSc }} -q -y + conda install mpich mpi4py petsc4py=${{ matrix.PETSc }} -q -y fi export OMPI_MCA_rmaps_base_oversubscribe=1 export OMPI_MCA_btl=^openib 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 829c43ec1..c89dca539 100644 --- a/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb +++ b/docs/dymos_book/examples/balanced_field/balanced_field_funccomp.ipynb @@ -168,6 +168,7 @@ "source": [ "import openmdao.api as om\n", "import openmdao.func_api as omf\n", + "from dymos.utils.misc import om_version", "\n", "def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True):\n", " \"\"\"\n", @@ -223,7 +224,10 @@ " meta.declare_coloring('*', method=grad_method)\n", " meta.declare_partials(of='*', wrt='*', method=grad_method)\n", " \n", - " return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit)\n", + " if om_version()[0] > (3, 35, 0):\n", + " return om.ExplicitFuncComp(meta, derivs_method=grad_method, use_jit=jax_jit)\n", + " else:\n", + " return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit)\n", " \n", " " ] diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png b/docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png index a431c9811..711521e8d 100644 Binary files a/docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png and b/docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png differ diff --git a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb index 4d0f96408..f34180b4c 100644 --- a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb +++ b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb @@ -100,9 +100,9 @@ "outputs": [], "source": [ "import numpy as np\n", - "from scipy.interpolate import interp1d\n", "\n", "import openmdao.api as om\n", + "from openmdao.components.interp_util.interp import InterpND\n", "\n", "import dymos as dm\n", "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n", @@ -158,7 +158,7 @@ "## The cannonball ODE component\n", "\n", "This component computes the state rates and the kinetic energy of the cannonball.\n", - "By calling the `declare_coloring` method wrt all inputs and using method `'cs'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using complex-step approximation." + "By calling the `declare_coloring` method wrt all inputs and using method `'fd'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using a finite-difference approximation." ] }, { @@ -198,17 +198,17 @@ " self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])\n", " self.add_output('ke', shape=nn, units='J')\n", "\n", - " # Ask OpenMDAO to compute the partial derivatives using complex-step\n", + " # Ask OpenMDAO to compute the partial derivatives using finite-difference\n", " # with a partial coloring algorithm for improved performance, and use\n", " # a graph coloring algorithm to automatically detect the sparsity pattern.\n", - " self.declare_coloring(wrt='*', method='cs')\n", + " self.declare_coloring(wrt='*', method='fd')\n", "\n", " alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]\n", " rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]\n", - " self.rho_interp = interp1d(np.array(alt_data, dtype=complex),\n", - " np.array(rho_data, dtype=complex),\n", - " kind='linear')\n", - "\n", + " self.rho_interp = InterpND(points=np.array(alt_data),\n", + " values=np.array(rho_data),\n", + " method='slinear')\n", + " \n", " def compute(self, inputs, outputs):\n", "\n", " gam = inputs['gam']\n", @@ -220,11 +220,7 @@ "\n", " GRAVITY = 9.80665 # m/s**2\n", "\n", - " # handle complex-step gracefully from the interpolant\n", - " if np.iscomplexobj(h):\n", - " rho = self.rho_interp(inputs['h'])\n", - " else:\n", - " rho = self.rho_interp(inputs['h']).real\n", + " rho = self.rho_interp.interpolate(h)\n", "\n", " q = 0.5*rho*inputs['v']**2\n", " qS = q * S\n", @@ -457,7 +453,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.0" + "version": "3.12.6" } }, "nbformat": 4, diff --git a/dymos/examples/balanced_field/balanced_field_length.py b/dymos/examples/balanced_field/balanced_field_length.py new file mode 100644 index 000000000..a6dd223a8 --- /dev/null +++ b/dymos/examples/balanced_field/balanced_field_length.py @@ -0,0 +1,189 @@ +import openmdao.api as om +import dymos as dm + + +def make_balanced_field_length_problem(ode_class, tx): + """ + Create a balanced field length problem and optionally set default values into it. + + Parameters + ---------- + ode_class : System class + The Dymos ODE System class. + tx_class : Transcription + Transcription to use. + + Returns + ------- + _type_ + _description_ + """ + p = om.Problem() + + p.driver = om.pyOptSparseDriver() + p.driver.declare_coloring() + + # Use IPOPT if available, with fallback to SLSQP + p.driver.options['optimizer'] = 'IPOPT' + p.driver.options['print_results'] = True + + p.driver.opt_settings['print_level'] = 0 + p.driver.opt_settings['mu_strategy'] = 'adaptive' + + p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' + p.driver.opt_settings['mu_init'] = 0.01 + p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' + + # First Phase: Brake release to V1 - both engines operable + br_to_v1 = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) + br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) + br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0) + br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') + br_to_v1.add_timeseries_output('*') + + # Second Phase: Rejected takeoff at V1 - no engines operable + rto = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) + rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + rto.add_parameter('alpha', val=0.0, opt=False, units='deg') + rto.add_timeseries_output('*') + + # Third Phase: V1 to Vr - single engine operable + v1_to_vr = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) + v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') + v1_to_vr.add_timeseries_output('*') + + # Fourth Phase: Rotate - single engine operable + rotate = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'runway'}) + rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) + rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, + val=[0, 10], control_type='polynomial') + rotate.add_timeseries_output('*') + + # Fifth Phase: Climb to target speed and altitude at end of runway. + climb = dm.Phase(ode_class=ode_class, transcription=tx, + ode_init_kwargs={'mode': 'climb'}) + climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) + climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) + climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0) + climb.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) + climb.add_state('gam', fix_initial=True, lower=0, ref=0.05, defect_ref=0.05) + climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10) + climb.add_timeseries_output('*') + + # Instantiate the trajectory and add phases + traj = dm.Trajectory() + p.model.add_subsystem('traj', traj) + traj.add_phase('br_to_v1', br_to_v1) + traj.add_phase('rto', rto) + traj.add_phase('v1_to_vr', v1_to_vr) + traj.add_phase('rotate', rotate) + traj.add_phase('climb', climb) + + all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb'] + groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate'] + + # Add parameters common to multiple phases to the trajectory + traj.add_parameter('m', val=174200., opt=False, units='lbm', + desc='aircraft mass', + targets={phase: ['m'] for phase in all_phases}) + + # Handle parameters which change from phase to phase. + traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, + desc='nominal aircraft thrust', + targets={'br_to_v1': ['T']}) + + traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, + desc='thrust under a single engine', + targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) + + traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, + desc='thrust when engines are shut down for rejected takeoff', + targets={'rto': ['T']}) + + traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, + desc='nominal runway friction coefficient', + targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) + + traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, + desc='runway friction coefficient under braking', + targets={'rto': ['mu_r']}) + + traj.add_parameter('h_runway', val=0., opt=False, units='ft', + desc='runway altitude', + targets={phase: ['h'] for phase in groundroll_phases}) + + # Here we're omitting some constants that are common throughout all phases for the sake of brevity. + # Their correct defaults are specified in add_input calls to `wrap_ode_func`. + + # Standard "end of first phase to beginning of second phase" linkages + # Alpha changes from being a parameter in v1_to_vr to a polynomial control + # in rotate, to a dynamic control in `climb`. + traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) + traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) + traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) + traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) + + # Less common "final value of r must match at ends of two phases". + traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', + phase_b='climb', var_b='r', loc_b='final', + ref=1000) + + # Define the constraints and objective for the optimal control problem + v1_to_vr.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=100) + + rto.add_boundary_constraint('v', loc='final', equals=0., ref=100, linear=True) + + rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) + + climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) + climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) + climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') + climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.25, ref=1.25) + + rto.add_objective('r', loc='final', ref=1000.0) + + # + # Setup the problem and set the initial guess + # + p.setup(check=True) + + br_to_v1.set_time_val(initial=0.0, duration=35.0) + br_to_v1.set_state_val('r', [0, 2500.0]) + br_to_v1.set_state_val('v', [0.0001, 100.0]) + br_to_v1.set_parameter_val('alpha', 0.0, units='deg') + + v1_to_vr.set_time_val(initial=35.0, duration=35.0) + v1_to_vr.set_state_val('r', [2500, 300.0]) + v1_to_vr.set_state_val('v', [100, 110.0]) + v1_to_vr.set_parameter_val('alpha', 0.0, units='deg') + + rto.set_time_val(initial=35.0, duration=1.0) + rto.set_state_val('r', [2500, 5000.0]) + rto.set_state_val('v', [110, 0.0001]) + rto.set_parameter_val('alpha', 0.0, units='deg') + + rotate.set_time_val(initial=35.0, duration=5.0) + rotate.set_state_val('r', [1750, 1800.0]) + rotate.set_state_val('v', [80, 85.0]) + rotate.set_control_val('alpha', 0.0, units='deg') + + climb.set_time_val(initial=30.0, duration=20.0) + climb.set_state_val('r', [5000, 5500.0], units='ft') + climb.set_state_val('v', [160, 170.0], units='kn') + climb.set_state_val('h', [0.0, 35.0], units='ft') + climb.set_state_val('gam', [0.0, 5.0], units='deg') + climb.set_control_val('alpha', 5.0, units='deg') + + return p diff --git a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py index 13b2edd6e..14c2b3caf 100644 --- a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py +++ b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py @@ -2,11 +2,14 @@ import openmdao.api as om import openmdao.func_api as omf +from dymos.utils.misc import om_version from openmdao.utils.assert_utils import assert_near_equal from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse import dymos as dm import numpy as np +from dymos.examples.balanced_field.balanced_field_length import make_balanced_field_length_problem + try: import jax except ImportError: @@ -77,12 +80,12 @@ def climb_ode(rho, S, CD0, CL0, CL_max, alpha_max, h_w, AR, e, span, T, m, v, h, return CL, q, L, D, K, F_r, v_dot, r_dot, W, v_stall, v_over_v_stall, gam_dot, h_dot -def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True): +def wrap_ode_func(num_nodes, mode, grad_method='jax', jax_jit=True): """ Returns the metadata from omf needed to create a new ExplciitFuncComp. """ nn = num_nodes - ode_func = runway_ode if flight_mode == 'runway' else climb_ode + ode_func = runway_ode if mode == 'runway' else climb_ode meta = (omf.wrap(ode_func) .add_input('rho', val=1.225, desc='atmospheric density at runway', units='kg/m**3') @@ -119,7 +122,7 @@ def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True): .add_output('v_over_v_stall', shape=(nn,), desc='stall speed ratio', units=None) ) - if flight_mode == 'runway': + if mode == 'runway': meta.add_input('mu_r', val=0.05, desc='runway friction coefficient', units=None) else: meta.add_input('gam', shape=(nn,), desc='flight path angle', units='rad') @@ -131,7 +134,10 @@ def wrap_ode_func(num_nodes, flight_mode, grad_method='jax', jax_jit=True): meta.declare_coloring('*', method=grad_method) meta.declare_partials(of='*', wrt='*', method=grad_method) - return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit) + if om_version()[0] > (3, 35, 0): + return om.ExplicitFuncComp(meta, derivs_method=grad_method, use_jit=jax_jit) + else: + return om.ExplicitFuncComp(meta, use_jax=grad_method == 'jax', use_jit=jax_jit) @use_tempdirs @@ -140,181 +146,18 @@ class TestBalancedFieldFuncComp(unittest.TestCase): @unittest.skipIf(jax is None, 'requires jax and jaxlib') @require_pyoptsparse('IPOPT') def test_balanced_field_func_comp_radau(self): - self._run_problem(dm.Radau) + p = make_balanced_field_length_problem(ode_class=wrap_ode_func, + tx=dm.Radau(num_segments=3)) + + dm.run_problem(p, run_driver=True, simulate=True) + + assert_near_equal(p.get_val('traj.rto.states:r')[-1], 2197.7, tolerance=0.01) @unittest.skipIf(jax is None, 'requires jax and jaxlib') @require_pyoptsparse('IPOPT') def test_balanced_field_func_comp_gl(self): - self._run_problem(dm.GaussLobatto) - - def _run_problem(self, tx): - p = om.Problem() - - p.driver = om.pyOptSparseDriver() - p.driver.declare_coloring() - - # Use IPOPT if available, with fallback to SLSQP - p.driver.options['optimizer'] = 'IPOPT' - p.driver.options['print_results'] = True - - p.driver.opt_settings['print_level'] = 0 - p.driver.opt_settings['mu_strategy'] = 'adaptive' - - p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' - p.driver.opt_settings['mu_init'] = 0.01 - p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) - br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) - br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0) - br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') - br_to_v1.add_timeseries_output('*') - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - rto.add_parameter('alpha', val=0.0, opt=False, units='deg') - rto.add_timeseries_output('*') - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') - v1_to_vr.add_timeseries_output('*') - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=3), - ode_init_kwargs={'flight_mode': 'runway'}) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) - rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, - val=[0, 10], control_type='polynomial') - rotate.add_timeseries_output('*') - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase(ode_class=wrap_ode_func, transcription=tx(num_segments=5), - ode_init_kwargs={'flight_mode': 'climb'}) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) - climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0) - climb.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - climb.add_state('gam', fix_initial=True, lower=0, ref=0.05, defect_ref=0.05) - climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10) - climb.add_timeseries_output('*') - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - traj.add_phase('br_to_v1', br_to_v1) - traj.add_phase('rto', rto) - traj.add_phase('v1_to_vr', v1_to_vr) - traj.add_phase('rotate', rotate) - traj.add_phase('climb', climb) - - all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb'] - groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate'] - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter('m', val=174200., opt=False, units='lbm', - desc='aircraft mass', - targets={phase: ['m'] for phase in all_phases}) - - # Handle parameters which change from phase to phase. - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, - desc='nominal aircraft thrust', - targets={'br_to_v1': ['T']}) - - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, - desc='thrust under a single engine', - targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, - desc='thrust when engines are shut down for rejected takeoff', - targets={'rto': ['T']}) - - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, - desc='nominal runway friction coefficient', - targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, - desc='runway friction coefficient under braking', - targets={'rto': ['mu_r']}) - - traj.add_parameter('h_runway', val=0., opt=False, units='ft', - desc='runway altitude', - targets={phase: ['h'] for phase in groundroll_phases}) - - # Here we're omitting some constants that are common throughout all phases for the sake of brevity. - # Their correct defaults are specified in add_input calls to `wrap_ode_func`. - - # Standard "end of first phase to beginning of second phase" linkages - # Alpha changes from being a parameter in v1_to_vr to a polynomial control - # in rotate, to a dynamic control in `climb`. - traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) - traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) - - # Less common "final value of r must match at ends of two phases". - traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', - phase_b='climb', var_b='r', loc_b='final', - ref=1000) - - # Define the constraints and objective for the optimal control problem - v1_to_vr.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=100) - - rto.add_boundary_constraint('v', loc='final', equals=0., ref=100, linear=True) - - rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) - - climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) - climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) - climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') - climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.25, ref=1.25) - - rto.add_objective('r', loc='final', ref=1000.0) - - # - # Setup the problem and set the initial guess - # - p.setup(check=True) - - br_to_v1.set_time_val(initial=0.0, duration=35.0) - br_to_v1.set_state_val('r', [0, 2500.0]) - br_to_v1.set_state_val('v', [0.0001, 100.0]) - br_to_v1.set_parameter_val('alpha', 0.0, units='deg') - - v1_to_vr.set_time_val(initial=35.0, duration=35.0) - v1_to_vr.set_state_val('r', [2500, 300.0]) - v1_to_vr.set_state_val('v', [100, 110.0]) - v1_to_vr.set_parameter_val('alpha', 0.0, units='deg') - - rto.set_time_val(initial=35.0, duration=1.0) - rto.set_state_val('r', [2500, 5000.0]) - rto.set_state_val('v', [110, 0.0001]) - rto.set_parameter_val('alpha', 0.0, units='deg') - - rotate.set_time_val(initial=35.0, duration=5.0) - rotate.set_state_val('r', [1750, 1800.0]) - rotate.set_state_val('v', [80, 85.0]) - rotate.set_control_val('alpha', 0.0, units='deg') - - climb.set_time_val(initial=30.0, duration=20.0) - climb.set_state_val('r', [5000, 5500.0], units='ft') - climb.set_state_val('v', [160, 170.0], units='kn') - climb.set_state_val('h', [0.0, 35.0], units='ft') - climb.set_state_val('gam', [0.0, 5.0], units='deg') - climb.set_control_val('alpha', 5.0, units='deg') + p = make_balanced_field_length_problem(ode_class=wrap_ode_func, + tx=dm.GaussLobatto(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=True) diff --git a/dymos/examples/balanced_field/test/test_balanced_field_length.py b/dymos/examples/balanced_field/test/test_balanced_field_length.py index 9fc428257..16f68e1a3 100644 --- a/dymos/examples/balanced_field/test/test_balanced_field_length.py +++ b/dymos/examples/balanced_field/test/test_balanced_field_length.py @@ -10,225 +10,25 @@ import dymos as dm from dymos.examples.balanced_field.balanced_field_ode import BalancedFieldODEComp +from dymos.examples.balanced_field.balanced_field_length import make_balanced_field_length_problem from dymos.utils.misc import om_version @use_tempdirs class TestBalancedFieldLengthRestart(unittest.TestCase): - def _make_problem(self): - p = om.Problem(reports=True) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'IPOPT' - p.driver.options['invalid_desvar_behavior'] = 'ignore' - p.driver.opt_settings['print_level'] = 0 - p.driver.opt_settings['derivative_test'] = 'first-order' - - p.driver.declare_coloring() - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) - br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) - br_to_v1.add_state('v', fix_initial=True, lower=0.0001, ref=100.0, defect_ref=100.0) - br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') - br_to_v1.add_timeseries_output('*') - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rto.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - rto.add_parameter('alpha', val=0.0, opt=False, units='deg') - rto.add_timeseries_output('*') - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - v1_to_vr.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') - v1_to_vr.add_timeseries_output('*') - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) - rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rotate.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, - val=[0, 10], control_type='polynomial') - rotate.add_timeseries_output('*') - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=5), - ode_init_kwargs={'mode': 'climb'}) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) - climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - climb.add_state('h', fix_initial=True, lower=0.0, ref=1.0, defect_ref=1.0) - climb.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0) - climb.add_state('gam', fix_initial=True, lower=0.0, ref=0.05, defect_ref=0.05) - climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10) - climb.add_timeseries_output('*') - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - traj.add_phase('br_to_v1', br_to_v1) - traj.add_phase('rto', rto) - traj.add_phase('v1_to_vr', v1_to_vr) - traj.add_phase('rotate', rotate) - traj.add_phase('climb', climb) - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter('m', val=174200., opt=False, units='lbm', - desc='aircraft mass', - targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'], - 'rotate': ['m'], 'climb': ['m']}) - - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, - desc='nominal aircraft thrust', - targets={'br_to_v1': ['T']}) - - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, - desc='thrust under a single engine', - targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, - desc='thrust when engines are shut down for rejected takeoff', - targets={'rto': ['T']}) - - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, - desc='nominal runway friction coeffcient', - targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, - desc='runway friction coefficient under braking', - targets={'rto': ['mu_r']}) - - traj.add_parameter('h_runway', val=0., opt=False, units='ft', static_target=False, - desc='runway altitude', - targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'], - 'rotate': ['h']}) - - traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True, - desc='atmospheric density', - targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'], - 'rotate': ['rho']}) - - traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True, - desc='aerodynamic reference area', - targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'], - 'rotate': ['S'], 'climb': ['S']}) - - traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True, - desc='zero-lift drag coefficient', - targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True, - desc='wing aspect ratio', - targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('e', val=801, opt=False, units=None, static_target=True, - desc='Oswald span efficiency factor', - targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True, - desc='wingspan', - targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True, - desc='height of wing above CG', - targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True, - desc='zero-alpha lift coefficient', - targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True, - desc='maximum lift coefficient for linear fit', - targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - # Standard "end of first phase to beginning of second phase" linkages - traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) - traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) - - # Less common "final value of r must be the match at ends of two phases". - traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', - phase_b='climb', var_b='r', loc_b='final', - ref=1000) - - # Define the constraints and objective for the optimal control problem - rto.add_boundary_constraint('v', loc='final', upper=0.001, ref=100, linear=True) - - rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) - - climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) - climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) - climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') - climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=1.2) - - rto.add_objective('r', loc='final', ref=1000.0) - - # - # Setup the problem and set the initial guess - # - p.setup(check=True) - - br_to_v1.set_time_val(initial=0.0, duration=35.0) - br_to_v1.set_state_val('r', [0, 2500.0]) - br_to_v1.set_state_val('v', [0.0001, 100.0]) - br_to_v1.set_parameter_val('alpha', 0.0, units='deg') - - v1_to_vr.set_time_val(initial=35.0, duration=35.0) - v1_to_vr.set_state_val('r', [2500, 300.0]) - v1_to_vr.set_state_val('v', [100, 110.0]) - v1_to_vr.set_parameter_val('alpha', 0.0, units='deg') - - rto.set_time_val(initial=35.0, duration=1.0) - rto.set_state_val('r', [2500, 5000.0]) - rto.set_state_val('v', [110, 0.0001]) - rto.set_parameter_val('alpha', 0.0, units='deg') - - rotate.set_time_val(initial=35.0, duration=5.0) - rotate.set_state_val('r', [1750, 1800.0]) - rotate.set_state_val('v', [80, 85.0]) - rotate.set_control_val('alpha', 0.0, units='deg') - - climb.set_time_val(initial=30.0, duration=20.0) - climb.set_state_val('r', [5000, 5500.0], units='ft') - climb.set_state_val('v', [160, 170.0], units='kn') - climb.set_state_val('h', [0.0, 35.0], units='ft') - climb.set_state_val('gam', [0.0, 5.0], units='deg') - climb.set_control_val('alpha', 5.0, units='deg') - - return p - @require_pyoptsparse(optimizer='IPOPT') @unittest.skipUnless(Version(openmdao.__version__) > Version("3.23"), reason='Test requires OpenMDAO 3.23.0 or later.') def test_make_plots(self): - p = self._make_problem() + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.Radau(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=True, make_plots=True) @require_pyoptsparse(optimizer='IPOPT') @unittest.skipUnless(Version(openmdao.__version__) > Version("3.23"), reason='Test requires OpenMDAO 3.23.0 or later.') def test_restart_from_sol(self): - p = self._make_problem() + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.Radau(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=False) sol_db = 'dymos_solution.db' @@ -237,7 +37,7 @@ def test_restart_from_sol(self): sol_results = om.CaseReader(sol_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) dm.run_problem(p, run_driver=True, simulate=True, restart=sol_db) @@ -250,17 +50,17 @@ def test_restart_from_sol(self): sol_results = om.CaseReader(sol_db).get_case('final') sim_results = om.CaseReader(sim_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) - assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) @require_pyoptsparse(optimizer='IPOPT') @unittest.skipUnless(Version(openmdao.__version__) > Version("3.23"), reason='Test requires OpenMDAO 3.23.0 or later.') def test_restart_from_sim(self): - p = self._make_problem() + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.Radau(num_segments=3)) dm.run_problem(p, run_driver=True, simulate=True) sol_db = 'dymos_solution.db' @@ -271,215 +71,39 @@ def test_restart_from_sim(self): sol_results = om.CaseReader(sol_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) dm.run_problem(p, run_driver=True, simulate=True, restart=sim_db) sol_results = om.CaseReader(sol_db).get_case('final') sim_results = om.CaseReader(sim_db).get_case('final') - assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.climb.timeseries.r')[-1], 2197, tolerance=0.01) - assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) - assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 1881, tolerance=0.01) + assert_near_equal(sol_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) + assert_near_equal(sim_results.get_val('traj.rto.timeseries.r')[-1], 2197, tolerance=0.01) @use_tempdirs class TestBalancedFieldLengthDefaultValues(unittest.TestCase): + @require_pyoptsparse(optimizer='IPOPT') def test_default_vals_stick(self): """ Make the balanced field problem without any set_val calls after setup. """ - p = om.Problem() - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0, - initial_val=0.0, duration_val=35.0) - br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0, - val=br_to_v1.interp(ys=[0, 2500.0], nodes='state_input')) - br_to_v1.add_state('v', fix_initial=True, lower=0.0001, ref=100.0, defect_ref=100.0, - val=br_to_v1.interp(ys=[0.0001, 100.0], nodes='state_input')) - br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg') - br_to_v1.add_timeseries_output('*') - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0, - initial_val=35.0, duration_val=35.0) - rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=rto.interp(ys=[2500, 5000.0], nodes='state_input')) - rto.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=rto.interp(ys=[110, 0.0001], nodes='state_input')) - rto.add_parameter('alpha', val=0.0, opt=False, units='deg') - rto.add_timeseries_output('*') - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0, - initial_val=35.0, duration_val=35.0) - v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=v1_to_vr.interp(ys=[2500, 300.0], nodes='state_input')) - v1_to_vr.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=v1_to_vr.interp(ys=[100, 110.0], nodes='state_input')) - v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg') - v1_to_vr.add_timeseries_output('*') - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=3), - ode_init_kwargs={'mode': 'runway'}) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0, - initial_val=70.0, duration_val=5.0) - rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=rotate.interp(ys=[1750, 1800.0], nodes='state_input')) - rotate.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=rotate.interp(ys=[80, 85.0], nodes='state_input')) - rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10, - val=[0, 10], control_type='polynomial') - rotate.add_timeseries_output('*') - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase(ode_class=BalancedFieldODEComp, transcription=dm.Radau(num_segments=5), - ode_init_kwargs={'mode': 'climb'}) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0, - initial_val=75.0, duration_val=10.0) - climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0, - val=climb.interp(ys=[1800, 2500.0], nodes='state_input')) - climb.add_state('h', fix_initial=True, lower=0.0, ref=1.0, defect_ref=1.0, val=0.0) - climb.add_state('v', fix_initial=False, lower=0.0001, ref=100.0, defect_ref=100.0, - val=climb.interp(ys=[85, 90], nodes='state_input')) - climb.add_state('gam', fix_initial=True, lower=0.0, ref=0.05, defect_ref=0.05, - val=climb.interp(ys=[0, 0.05], nodes='state_input')) - climb.add_control('alpha', opt=True, units='deg', lower=-10, upper=15, ref=10, - val=climb.interp(ys=[0.01, 0.01], nodes='control_input')) - climb.add_timeseries_output('*') - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - traj.add_phase('br_to_v1', br_to_v1) - traj.add_phase('rto', rto) - traj.add_phase('v1_to_vr', v1_to_vr) - traj.add_phase('rotate', rotate) - traj.add_phase('climb', climb) - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter('m', val=174200., opt=False, units='lbm', - desc='aircraft mass', - targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'], - 'rotate': ['m'], 'climb': ['m']}) - - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, - desc='nominal aircraft thrust', targets={'br_to_v1': ['T']}) - - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, - desc='thrust under a single engine', - targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, - desc='thrust when engines are shut down for rejected takeoff', - targets={'rto': ['T']}) - - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, - desc='nominal runway friction coeffcient', - targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, - desc='runway friction coefficient under braking', - targets={'rto': ['mu_r']}) - - traj.add_parameter('h_runway', val=0., opt=False, units='ft', static_target=False, - desc='runway altitude', - targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'], - 'rotate': ['h']}) - - traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True, - desc='atmospheric density', - targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'], - 'rotate': ['rho']}) - - traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True, - desc='aerodynamic reference area', - targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'], - 'rotate': ['S'], 'climb': ['S']}) - - traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True, - desc='zero-lift drag coefficient', - targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True, - desc='wing aspect ratio', - targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('e', val=801, opt=False, units=None, static_target=True, - desc='Oswald span efficiency factor', - targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True, - desc='wingspan', - targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True, - desc='height of wing above CG', - targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True, - desc='zero-alpha lift coefficient', - targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True, - desc='maximum lift coefficient for linear fit', - targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr', - 'rto', 'rotate', 'climb']}) - - # Standard "end of first phase to beginning of second phase" linkages - traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v']) - traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha']) - traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v']) - - # Less common "final value of r must be the match at ends of two phases". - traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final', - phase_b='climb', var_b='r', loc_b='final', - ref=1000) - - # Define the constraints and objective for the optimal control problem - rto.add_boundary_constraint('v', loc='final', upper=0.001, ref=100, linear=True) - - rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000) - - climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True) - climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True) - climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg') - climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=1.2) - - rto.add_objective('r', loc='final', ref=1000.0) - - # - # Setup the problem and set the initial guess - # - p.setup(check=True) + p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.GaussLobatto(num_segments=3)) p.run_model() - assert_near_equal(p.get_val('traj.rotate.t_initial'), 70) + assert_near_equal(p.get_val('traj.rotate.t_initial'), 35) assert_near_equal(p.get_val('traj.rotate.t_duration'), 5) - assert_near_equal(p.get_val('traj.rotate.controls:alpha'), np.array([[0, 10]]).T) - assert_near_equal(p.get_val('traj.climb.controls:alpha'), - p.model.traj.phases.climb.interp('', [0.01, 0.01], nodes='control_input')) - assert_near_equal(p.get_val('traj.climb.states:gam'), - p.model.traj.phases.climb.interp(ys=[0.0, 0.05], nodes='state_input')) - assert_near_equal(p.get_val('traj.climb.states:h'), - p.model.traj.phases.climb.interp(ys=[0.0, 0.0], nodes='state_input')) + assert_near_equal(p.get_val('traj.rotate.controls:alpha'), np.array([[0, 0]]).T) + assert_near_equal(p.get_val('traj.climb.controls:alpha', units='deg'), + p.model.traj.phases.climb.interp('', [5, 5], nodes='control_input')) + assert_near_equal(p.get_val('traj.climb.states:gam', units='deg'), + p.model.traj.phases.climb.interp(ys=[0.0, 5.0], nodes='state_input')) + assert_near_equal(p.get_val('traj.climb.states:h', units='ft'), + p.model.traj.phases.climb.interp(ys=[0.0, 35.0], nodes='state_input')) assert_near_equal(p.get_val('traj.v1_to_vr.parameters:alpha'), 0.0) diff --git a/dymos/examples/cannonball/cannonball_ode.py b/dymos/examples/cannonball/cannonball_ode.py index 7e9bded31..b9b97d1f7 100644 --- a/dymos/examples/cannonball/cannonball_ode.py +++ b/dymos/examples/cannonball/cannonball_ode.py @@ -1,14 +1,16 @@ import numpy as np -from scipy.interpolate import interp1d import openmdao.api as om +from openmdao.components.interp_util.interp import InterpND + from dymos.models.atmosphere.atmos_1976 import USatm1976Data english_to_metric_rho = om.unit_conversion('slug/ft**3', 'kg/m**3')[0] english_to_metric_h = om.unit_conversion('ft', 'm')[0] -rho_interp = interp1d(np.array(USatm1976Data.alt*english_to_metric_h, dtype=complex), - np.array(USatm1976Data.rho*english_to_metric_rho, dtype=complex), kind='linear') +rho_interp = InterpND(points=np.array(USatm1976Data.alt*english_to_metric_h), + values=np.array(USatm1976Data.rho*english_to_metric_rho), + method='slinear').interpolate class CannonballODE(om.ExplicitComponent): @@ -56,11 +58,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = rho_interp(inputs['h']) - else: - rho = rho_interp(inputs['h']).real + rho = rho_interp(h) q = 0.5*rho*inputs['v']**2 qS = q * S diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py index 194dd79b5..ddb0f9b55 100644 --- a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py +++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py @@ -24,9 +24,8 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='IPOPT') def test_two_phase_cannonball_for_docs(self): - from scipy.interpolate import interp1d - import openmdao.api as om + from openmdao.components.interp_util.interp import InterpND from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -107,16 +106,16 @@ def setup(self): self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) self.add_output('ke', shape=nn, units='J') - # Ask OpenMDAO to compute the partial derivatives using complex-step + # Ask OpenMDAO to compute the partial derivatives using finite-difference # with a partial coloring algorithm for improved performance, and use # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='cs') + self.declare_coloring(wrt='*', method='fd') alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = interp1d(np.array(alt_data, dtype=complex), - np.array(rho_data, dtype=complex), - kind='linear') + self.rho_interp = InterpND(points=np.array(alt_data), + values=np.array(rho_data), + method='slinear') def compute(self, inputs, outputs): @@ -129,13 +128,9 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = self.rho_interp(inputs['h']) - else: - rho = self.rho_interp(inputs['h']).real + rho = self.rho_interp.interpolate(h) - q = 0.5*rho*inputs['v']**2 + q = 0.5*rho*v**2 qS = q * S D = qS * CD cgam = np.cos(gam) diff --git a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py index 96909213e..81ebd9008 100644 --- a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py +++ b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py @@ -9,9 +9,9 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='SLSQP') def test_two_phase_cannonball_for_docs(self): import numpy as np - from scipy.interpolate import interp1d import openmdao.api as om + from openmdao.components.interp_util.interp import InterpND from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -92,16 +92,16 @@ def setup(self): self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) self.add_output('ke', shape=nn, units='J') - # Ask OpenMDAO to compute the partial derivatives using complex-step + # Ask OpenMDAO to compute the partial derivatives using finite-difference # with a partial coloring algorithm for improved performance, and use # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='cs') + self.declare_coloring(wrt='*', method='fd') alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = interp1d(np.array(alt_data, dtype=complex), - np.array(rho_data, dtype=complex), - kind='linear') + self.rho_interp = InterpND(points=np.array(alt_data), + values=np.array(rho_data), + method='slinear') def compute(self, inputs, outputs): @@ -114,11 +114,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = self.rho_interp(inputs['h']) - else: - rho = self.rho_interp(inputs['h']).real + rho = self.rho_interp.interpolate(inputs['h']) q = 0.5*rho*inputs['v']**2 qS = q * S diff --git a/dymos/examples/cannonball/test/test_load_case_missing_phase.py b/dymos/examples/cannonball/test/test_load_case_missing_phase.py index 0aa20ec74..a5aedb161 100644 --- a/dymos/examples/cannonball/test/test_load_case_missing_phase.py +++ b/dymos/examples/cannonball/test/test_load_case_missing_phase.py @@ -184,7 +184,7 @@ def test_load_case_missing_phase(self): 848.54403903, 864.60323344, 895.67083544, 926.90913594, 934.25949414, 946.77982935, 954.88116543, 955.37021704]) - assert_near_equal(p.get_val('traj.ascent.states:h').ravel(), h_expected, tolerance=1.0E-6) + assert_near_equal(p.get_val('traj.ascent.states:h').ravel(), h_expected, tolerance=1.0E-5) if __name__ == '__main__': # pragma: no cover diff --git a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py index 1f8e3b87a..48a8d0850 100644 --- a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py +++ b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py @@ -3,9 +3,9 @@ from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse import numpy as np -from scipy.interpolate import interp1d import openmdao.api as om +from openmdao.components.interp_util.interp import InterpND from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -83,16 +83,16 @@ def setup(self): self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) self.add_output('ke', shape=nn, units='J') - # Ask OpenMDAO to compute the partial derivatives using complex-step + # Ask OpenMDAO to compute the partial derivatives using finite-difference # with a partial coloring algorithm for improved performance, and use # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='cs') + self.declare_coloring(wrt='*', method='fd') alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = interp1d(np.array(alt_data, dtype=complex), - np.array(rho_data, dtype=complex), - kind='linear', bounds_error=False, fill_value='extrapolate') + self.rho_interp = InterpND(points=np.array(alt_data), + values=np.array(rho_data), + method='slinear', extrapolate=True) def compute(self, inputs, outputs): @@ -105,11 +105,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - # handle complex-step gracefully from the interpolant - if np.iscomplexobj(h): - rho = self.rho_interp(inputs['h']) - else: - rho = self.rho_interp(inputs['h']).real + rho = self.rho_interp.interpolate(inputs['h']) q = 0.5 * rho * inputs['v'] ** 2 qS = q * S diff --git a/pyproject.toml b/pyproject.toml index a7c048d95..d75b1c4c9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ all = [ "dymos[docs]", "dymos[notebooks]", "dymos[test]", + "dymos[numba]", ] docs = [ "bokeh", @@ -62,6 +63,8 @@ test = [ "playwright>=1.20", "pycodestyle", "testflo>=1.3.6", +] +numba = [ "numba", ]