Skip to content

Commit

Permalink
Fixed an issue where a DirectSolver was automatically added to a phas…
Browse files Browse the repository at this point in the history
…e if any state has `input_initial=True`. (#982)

- Having a connected initial state value is not itself enough to require a linear solver.
- Added an test to verify input_initial=True works without a direct solver.
- Fixed some logic that was causing errant warnings due in the handling of static_target parameter option.
  • Loading branch information
robfalck authored Sep 6, 2023
1 parent a129652 commit e2d882e
Show file tree
Hide file tree
Showing 4 changed files with 127 additions and 5 deletions.
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
import unittest
import matplotlib.pyplot as plt

from openmdao.utils.testing_utils import use_tempdirs
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal

import dymos as dm
from dymos.examples.cannonball.size_comp import CannonballSizeComp
from dymos.examples.cannonball.cannonball_ode import CannonballODE


plt.switch_backend('Agg')


@use_tempdirs
class TestCannonballConnectedInitialStates(unittest.TestCase):

def test_cannonball_connected_initial_states(self):

p = om.Problem(model=om.Group())

p.driver = om.ScipyOptimizeDriver()
p.driver.options['optimizer'] = 'SLSQP'
p.driver.declare_coloring()

p.model.add_subsystem('size_comp', CannonballSizeComp(),
promotes_inputs=['radius', 'dens'])
p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
p.model.add_design_var('radius', lower=0.01, upper=0.10,
ref0=0.01, ref=0.10, units='m')

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

transcription = dm.Radau(num_segments=5, order=3, compressed=True)
ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription)

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

ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
duration_ref=100, units='s')

ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), duration_ref=100, units='s')
ascent.add_state('r', fix_initial=True, fix_final=False, rate_source='r_dot', units='m')
ascent.add_state('h', fix_initial=True, fix_final=False, units='m', rate_source='h_dot')
ascent.add_state('gam', fix_initial=False, fix_final=True, units='rad', rate_source='gam_dot')
ascent.add_state('v', fix_initial=False, fix_final=False, units='m/s', rate_source='v_dot')

ascent.add_parameter('S', units='m**2', static_targets=True)
ascent.add_parameter('m', units='kg', static_targets=True)

# Limit the muzzle energy
ascent.add_boundary_constraint('ke = 0.5 * m * v**2', loc='initial',
upper=400000, lower=0, ref=100000)

# Second Phase (descent)
transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)

traj.add_phase('descent', descent)

descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
duration_ref=100, units='s')
descent.add_state('r', units='m', rate_source='r_dot', input_initial=True)
descent.add_state('h', units='m', rate_source='h_dot', input_initial=True, fix_final=True)
descent.add_state('gam', units='rad', rate_source='gam_dot', input_initial=True, fix_final=False)
descent.add_state('v', units='m/s', rate_source='v_dot', input_initial=True, fix_final=False)

descent.add_parameter('S', units='m**2', static_targets=True)
descent.add_parameter('m', units='kg', static_targets=True)

descent.add_objective('r', loc='final', scaler=-1.0)

traj.add_parameter('CD',
targets={'ascent': ['CD'], 'descent': ['CD']},
val=0.5, units=None, opt=False, static_targets=True)

traj.add_parameter('m', units='kg', val=1.0,
targets={'ascent': 'm', 'descent': 'm'}, static_targets=True)

traj.add_parameter('S', units='m**2', val=0.005, static_targets=True)

# Link Phases (link time and all state variables)
traj.link_phases(phases=['ascent', 'descent'], vars=['time'])
traj.link_phases(phases=['ascent', 'descent'], vars=['r', 'h', 'v', 'gam'], connected=True)

p.model.connect('size_comp.mass', 'traj.parameters:m')
p.model.connect('size_comp.S', 'traj.parameters:S')

p.model.linear_solver = om.DirectSolver()

# Finish Problem Setup
p.setup()

# Set constants and initial guesses
p.set_val('radius', 0.05, units='m')
p.set_val('dens', 7.87, units='g/cm**3')

p.set_val('traj.parameters:CD', 0.5)

p.set_val('traj.ascent.t_initial', 0.0)
p.set_val('traj.ascent.t_duration', 10.0)

p.set_val('traj.ascent.states:r', ascent.interp(ys=[0, 100], nodes='state_input'))
p.set_val('traj.ascent.states:h', ascent.interp(ys=[0, 100], nodes='state_input'))
p.set_val('traj.ascent.states:v', ascent.interp(ys=[200, 150], nodes='state_input'))
p.set_val('traj.ascent.states:gam', ascent.interp(ys=[25, 0], nodes='state_input'), units='deg')

p.set_val('traj.descent.t_initial', 10.0)
p.set_val('traj.descent.t_duration', 10.0)

p.set_val('traj.descent.states:r', descent.interp(ys=[100, 200], nodes='state_input'))
p.set_val('traj.descent.states:h', descent.interp(ys=[100, 0], nodes='state_input'))
p.set_val('traj.descent.states:v', descent.interp(ys=[150, 200], nodes='state_input'))
p.set_val('traj.descent.states:gam', descent.interp(ys=[0, -45], nodes='state_input'), units='deg')

# Run the optimization and final explicit simulation
dm.run_problem(p)

assert_near_equal(p.get_val('traj.descent.states:r')[-1],
3183.25, tolerance=1.0E-2)


if __name__ == '__main__':
unittest.main()
2 changes: 1 addition & 1 deletion dymos/phase/options.py
Original file line number Diff line number Diff line change
Expand Up @@ -268,7 +268,7 @@ def __init__(self, read_only=False):
'each node in the ODE.'
'If _unspecified, attempt to determine through introspection.',
deprecation='Use option `static_targets` to specify whether all targets\n'
'are static (static_targegts=True), none are static (static_targets=False),\n'
'are static (static_targets=True), none are static (static_targets=False),\n'
'static_targets are determined via introspection (static_targets=_unspecified),\n'
'or give an explicit sequence of the static targets.')

Expand Down
3 changes: 0 additions & 3 deletions dymos/phase/phase.py
Original file line number Diff line number Diff line change
Expand Up @@ -1050,15 +1050,12 @@ def set_parameter_options(self, name, val=_unspecified, units=_unspecified, opt=
self.parameter_options[name]['shape'] = shape

if dynamic is not _unspecified:
self.parameter_options[name]['static_target'] = not dynamic
self.parameter_options[name]['static_targets'] = not dynamic

if static_target is not _unspecified:
self.parameter_options[name]['static_target'] = static_target
self.parameter_options[name]['static_targets'] = static_target

if static_targets is not _unspecified:
self.parameter_options[name]['static_target'] = static_targets
self.parameter_options[name]['static_targets'] = static_targets

if dynamic is not _unspecified and static_target is not _unspecified:
Expand Down
2 changes: 1 addition & 1 deletion dymos/transcriptions/pseudospectral/pseudospectral_base.py
Original file line number Diff line number Diff line change
Expand Up @@ -514,7 +514,7 @@ def configure_solvers(self, phase):
newton.linesearch = om.BoundsEnforceLS()

# even though you don't need a nl_solver for connections, you still ln_solver since its implicit
if self.any_solved_segs or self.any_connected_opt_segs or self._implicit_duration:
if self.any_solved_segs or self._implicit_duration:
if isinstance(phase.linear_solver, om.LinearRunOnce):
phase.linear_solver = om.DirectSolver()

Expand Down

0 comments on commit e2d882e

Please sign in to comment.