Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/OpenMDAO/dymos into times…
Browse files Browse the repository at this point in the history
…eries_multi_traj_fix
  • Loading branch information
robfalck committed Nov 1, 2024
2 parents a0daa2c + 8338af1 commit e19d474
Show file tree
Hide file tree
Showing 14 changed files with 304 additions and 640 deletions.
6 changes: 3 additions & 3 deletions .github/workflows/dymos_docs_workflow.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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' &&
Expand Down
26 changes: 23 additions & 3 deletions .github/workflows/dymos_tests_workflow.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
" "
]
Expand Down
Binary file modified docs/dymos_book/examples/brachistochrone/brachistochrone_fbd.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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."
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -457,7 +453,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.0"
"version": "3.12.6"
}
},
"nbformat": 4,
Expand Down
189 changes: 189 additions & 0 deletions dymos/examples/balanced_field/balanced_field_length.py
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit e19d474

Please sign in to comment.