Skip to content

Commit

Permalink
Added better logic to detect different units across linkage constrain…
Browse files Browse the repository at this point in the history 
…ts. (#913)

* Added better logic to detect different units across linkage constraints.
  • Loading branch information
@robfalck
robfalck authored Mar 22, 2023
1 parent 2822281 commit 08b625e
Show file tree
Hide file tree
Showing 2 changed files with 121 additions and 6 deletions.
112 changes: 110 additions & 2 deletions dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,12 +121,120 @@ def test_two_burn_orbit_raise_gl_radau_gl_changing_units_error(self):
with self.assertRaises(ValueError) as e:
p.setup(check=True, force_alloc_complex=True)

expected_exception = 'traj: Linkage units were not specified but the units of var_a (DU/TU**2) and var_b ' \
'(km/s**2) are not the same. Units for this linkage constraint must ' \
expected_exception = 'traj: Linkage units were not specified but the units of burn1.accel (DU/TU**2) and ' \
'burn2.accel (km/s**2) are not equivalent. Units for this linkage constraint must ' \
'be specified explicitly.'

self.assertEqual(expected_exception, str(e.exception))

def test_two_burn_orbit_raise_gl_radau_gl_equivalent_units_error(self):
import openmdao.api as om

import dymos as dm
from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE

traj = dm.Trajectory()
p = om.Problem(model=om.Group())
p.model.add_subsystem('traj', traj)

p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = 'IPOPT'

p.driver.declare_coloring()

traj.add_parameter('c', opt=False, val=1.5, units='DU/TU',
targets={'burn1': ['c'], 'coast': ['c'], 'burn2': ['c']})

# First Phase (burn)

burn1 = dm.Phase(ode_class=FiniteBurnODE,
transcription=dm.GaussLobatto(num_segments=10, order=3, compressed=True))

burn1 = traj.add_phase('burn1', burn1)

burn1.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='TU')
burn1.add_state('r', fix_initial=True, fix_final=False, defect_scaler=100.0,
rate_source='r_dot', targets=['r'], units='nmi')
burn1.add_state('theta', fix_initial=True, fix_final=False, defect_scaler=100.0,
rate_source='theta_dot', targets=['theta'], units='rad')
burn1.add_state('vr', fix_initial=True, fix_final=False, defect_scaler=100.0,
rate_source='vr_dot', targets=['vr'], units='DU/TU')
burn1.add_state('vt', fix_initial=True, fix_final=False, defect_scaler=100.0,
rate_source='vt_dot', targets=['vt'], units='DU/TU')
burn1.add_state('accel', fix_initial=True, fix_final=False,
rate_source='at_dot', targets=['accel'], units='DU/TU**2')
burn1.add_state('deltav', fix_initial=True, fix_final=False,
rate_source='deltav_dot', units='DU/TU')
burn1.add_control('u1', rate_continuity=True, rate2_continuity=True, units='deg',
scaler=0.01,
rate_continuity_scaler=0.001, rate2_continuity_scaler=0.001,
lower=-30, upper=30, targets=['u1'])

# Second Phase (Coast)
coast = dm.Phase(ode_class=FiniteBurnODE,
transcription=dm.Radau(num_segments=10, order=3))

traj.add_phase('coast', coast)

coast.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 10), duration_ref=10, units='TU')

coast.add_state('r', fix_initial=False, fix_final=False, defect_scaler=100.0,
rate_source='r_dot', targets=['r'], units='nmi')
coast.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=100.0,
units='rad', rate_source='theta_dot', targets=['theta'])
coast.add_state('vr', fix_initial=False, fix_final=False, defect_scaler=100.0,
rate_source='vr_dot', targets=['vr'], units='DU/TU')
coast.add_state('vt', fix_initial=False, fix_final=False, defect_scaler=100.0,
rate_source='vt_dot', targets=['vt'], units='DU/TU')
coast.add_state('accel', fix_initial=True, fix_final=False, ref=1.0E-12, defect_ref=1.0E-12,
rate_source='at_dot', targets=['accel'], units='DU/TU**2')
coast.add_state('deltav', fix_initial=False, fix_final=False,
rate_source='deltav_dot', units='DU/TU')
coast.add_parameter('u1', targets=['u1'], opt=False, val=0.0, units='deg')

# Third Phase (burn)

burn2 = dm.Phase(ode_class=FiniteBurnODE,
transcription=dm.GaussLobatto(num_segments=10, order=3, compressed=True))

traj.add_phase('burn2', burn2)

burn2.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 10), initial_ref=10, units='TU')
burn2.add_state('r', fix_initial=False, fix_final=True,
rate_source='r_dot', targets=['r'], units='NM')
burn2.add_state('theta', fix_initial=False, fix_final=False,
rate_source='theta_dot', targets=['theta'], units='rad')
burn2.add_state('vr', fix_initial=False, fix_final=True,
rate_source='vr_dot', targets=['vr'], units='DU/TU')
burn2.add_state('vt', fix_initial=False, fix_final=True,
rate_source='vt_dot', targets=['vt'], units='DU/TU')
burn2.add_state('accel', fix_initial=False, fix_final=False, defect_ref=1.0E-6,
rate_source='at_dot', targets=['accel'], units='DU/TU**2')
burn2.add_state('deltav', fix_initial=False, fix_final=False,
rate_source='deltav_dot', units='DU/TU')
burn2.add_control('u1', targets=['u1'], rate_continuity=True, rate2_continuity=True,
units='deg', scaler=0.01, lower=-30, upper=30)

burn2.add_objective('deltav', loc='final', scaler=1.0)

burn1.add_timeseries_output('pos_x')
coast.add_timeseries_output('pos_x')
burn2.add_timeseries_output('pos_x')

burn1.add_timeseries_output('pos_y')
coast.add_timeseries_output('pos_y')
burn2.add_timeseries_output('pos_y')

# Link Phases
traj.link_phases(phases=['burn1', 'coast', 'burn2'],
vars=['time', 'r', 'vr', 'vt', 'deltav'])
traj.link_phases(phases=['burn1', 'burn2'], vars=['accel'])

# Finish Problem Setup
p.model.linear_solver = om.DirectSolver()

p.setup(check=True, force_alloc_complex=True)

def test_two_burn_orbit_raise_gl_radau_gl_constrained(self):
import numpy as np

Expand Down
15 changes: 11 additions & 4 deletions dymos/trajectory/trajectory.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

from openmdao.utils.om_warnings import warn_deprecation
from openmdao.utils.graph_utils import get_sccs_topo
from openmdao.utils.units import unit_conversion

import numpy as np
import networkx as nx
Expand Down Expand Up @@ -596,10 +597,16 @@ def _update_linkage_options_configure(self, linkage_options):
if linkage_options['units_b'] is _unspecified:
linkage_options['units_b'] = units['b']

if not linkage_options['connected'] and (linkage_options['units_a'] != linkage_options['units_b']) and \
linkage_options['units'] is _unspecified:
raise ValueError(f'{info_str}Linkage units were not specified but the units of '
f'var_a ({units["a"]}) and var_b ({units["b"]}) are not the same. '
if units['a'] is not None and units['b'] is not None:
conversion_scaler, conversion_offset = unit_conversion(units['a'], units['b'])
else:
conversion_scaler, conversion_offset = (1.0, 0.0)

if not linkage_options['connected'] \
and linkage_options['units'] is _unspecified \
and (abs(conversion_scaler - 1.0) > 1.0E-15 or abs(conversion_offset) > 1.0E-15):
raise ValueError(f'{info_str}Linkage units were not specified but the units of {phase_name_a}.{var_a} '
f'({units["a"]}) and {phase_name_b}.{var_b} ({units["b"]}) are not equivalent. '
f'Units for this linkage constraint must be specified explicitly.')
else:
linkage_options['units'] = units['a']
Expand Down

0 comments on commit 08b625e

Please sign in to comment.