Skip to content

Commit

Permalink
Fixed an issue where applying DirectSolver to StateIndependentsComp w…
Browse files Browse the repository at this point in the history 
…as breaking when used with other linear solvers under MPI. (#1020)

* Removed DirectSolver from StateIndependentsComp in the pseudospectral transcriptions. Phase now has default_linear_solver and default_nonlinear_solver options like trajectory.

* Don't constraint type of solver applied to Trajectory so users can use other solvers potentially unknown to us

* When Trajectory.phases is a ParallelGroup, always add a NonlinearBlockJac if there are direct connections.

* Removed default_(non)linear_solver options.
Added auto_solvers options to Trajectory and Phases, allowing users to disable the behavior if desired.
Added parallel_phases option to Trajectory, which toggles whether the top level phase container is a parallel group.

* Made consistent configure_solvers behavior in Birkhoff and ExplicitShooting
  • Loading branch information
@robfalck
robfalck authored Nov 29, 2023
1 parent 10c2df0 commit 37ac4e5
Show file tree
Hide file tree
Showing 13 changed files with 262 additions and 294 deletions.
77 changes: 27 additions & 50 deletions docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true,
"tags": [
"active-ipynb",
"remove-input",
Expand Down Expand Up @@ -76,6 +77,7 @@
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true,
"tags": [
"output_scroll"
]
Expand Down Expand Up @@ -148,7 +150,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"class CannonballODE(om.ExplicitComponent):\n",
Expand Down Expand Up @@ -232,7 +236,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"def initial_guess(t_dur, gam0=np.pi/3):\n",
Expand Down Expand Up @@ -265,10 +271,12 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"p = om.Problem(model=om.Group())\n",
"p = om.Problem(name='cannonball_implicit_duration', model=om.Group())\n",
"\n",
"p.driver = om.pyOptSparseDriver()\n",
"p.driver.options['optimizer'] = 'IPOPT'\n",
Expand Down Expand Up @@ -321,6 +329,14 @@
"# A nonlinear solver is used to find the duration of required to satisfy the condition.\n",
"phase.set_duration_balance('h', val=0.0)\n",
"\n",
"# Prior to setup, this phase has the default NonlinearRunOnce and LinearRunOnce solvers.\n",
"# Because of the implicit behavior introduced by set_duration_balance, it will automatically\n",
"# use a NewtonSolver with an Armijo Goldstein linesearch.\n",
"# In this case that linesearch causes extrapolation of the atmospheric model into invalid regions.\n",
"# To account for this, we explicitly assign a NewtonSolver with a different linesearch.\n",
"phase.nonlinear_solver = om.NewtonSolver(iprint=0, solve_subsystems=True)\n",
"phase.nonlinear_solver.linesearch = None\n",
"\n",
"phase.add_objective('r', loc='final', scaler=-1.0)\n",
"\n",
"p.model.connect('size_comp.mass', 'traj.phase.parameters:m')\n",
Expand Down Expand Up @@ -350,7 +366,7 @@
"#####################################################\n",
"# Run the optimization and final explicit simulation\n",
"#####################################################\n",
"dm.run_problem(p)"
"dm.run_problem(p, simulate=True, make_plots=True)"
]
},
{
Expand All @@ -363,60 +379,21 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"exp_out = traj.simulate()\n",
"\n",
"rad = p.get_val('radius', units='m')[0]\n",
"print(f'optimal radius: {rad} m ')\n",
"mass = p.get_val('size_comp.mass', units='kg')[0]\n",
"print(f'cannonball mass: {mass} kg ')\n",
"angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]\n",
"print(f'launch angle: {angle} deg')\n",
"max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]\n",
"print(f'maximum range: {max_range} m')\n",
"\n",
"fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))\n",
"\n",
"time_imp = p.get_val('traj.phase.timeseries.time')\n",
"\n",
"time_exp = exp_out.get_val('traj.phase.timeseries.time')\n",
"\n",
"r_imp = p.get_val('traj.phase.timeseries.r')\n",
"\n",
"r_exp = exp_out.get_val('traj.phase.timeseries.r')\n",
"\n",
"h_imp = p.get_val('traj.phase.timeseries.h')\n",
"\n",
"h_exp = exp_out.get_val('traj.phase.timeseries.h')\n",
"\n",
"axes.plot(r_imp, h_imp, 'bo')\n",
"\n",
"axes.plot(r_exp, h_exp, 'r--')\n",
"\n",
"axes.set_xlabel('range (m)')\n",
"axes.set_ylabel('altitude (m)')\n",
"\n",
"fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))\n",
"states = ['r', 'h', 'v', 'gam']\n",
"for i, state in enumerate(states):\n",
" x_imp = p.get_val(f'traj.phase.timeseries.{state}')\n",
"\n",
" x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')\n",
"\n",
" axes[i].set_ylabel(state)\n",
"\n",
" axes[i].plot(time_imp, x_imp, 'bo')\n",
" axes[i].plot(time_exp, x_exp, 'r--')\n",
"from IPython.display import IFrame\n",
"\n",
"plt.show()"
"IFrame(src='reports/cannonball_implicit_duration/traj_results_report.html', width=700, height=1000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true,
"tags": [
"remove-input",
"remove-output"
Expand Down
101 changes: 100 additions & 1 deletion dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,20 @@
from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse

try:
from petsc4py import PETSc
except ImportError:
PETSc = None

from openmdao.utils.general_utils import set_pyoptsparse_opt

OPT, OPTIMIZER = set_pyoptsparse_opt('SLSQP')


def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
compressed=True, optimizer='SLSQP', force_alloc_complex=False,
solve_segments=False, y_bounds=None):
solve_segments=False, y_bounds=None, phase_nonlinear_solver=None,
phase_linear_solver=None):
p = om.Problem(model=om.Group())

p.driver = om.pyOptSparseDriver()
Expand Down Expand Up @@ -44,6 +50,12 @@ def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_o

traj = dm.Trajectory()
phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)

if phase_nonlinear_solver is not None:
phase.nonlinear_solver = phase_nonlinear_solver
if phase_linear_solver is not None:
phase.linear_solver = phase_linear_solver

p.model.add_subsystem('traj0', traj)
traj.add_phase('phase0', phase)

Expand Down Expand Up @@ -461,3 +473,90 @@ def test_brachistochrone_bounded_solve_segments(self):

warnings.simplefilter('always')
self.assertIn(expected_warning, [str(w.message) for w in ctx])

@unittest.skipIf(PETSc is None, 'PETSc is not available.')
def test_brachistochrone_solve_segments_radau_krylov_solver(self, optimizer='SLSQP'):
#
# Initialize the Problem and the optimization driver
#
p = om.Problem(model=om.Group())

p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = optimizer
if optimizer == 'SNOPT':
p.driver.opt_settings['iSumm'] = 6
p.driver.opt_settings['Verify level'] = 3
elif optimizer == 'IPOPT':
p.driver.opt_settings['mu_init'] = 1e-3
p.driver.opt_settings['max_iter'] = 500
p.driver.opt_settings['print_level'] = 0
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' # for faster convergence
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['mu_strategy'] = 'monotone'
p.driver.declare_coloring(tol=1.0E-12)

#
# Create a trajectory and add a phase to it
#
traj = p.model.add_subsystem('traj0', dm.Trajectory())

# -------------------------------
# forward solver-based shooting
# -------------------------------
phase0 = dm.Phase(ode_class=BrachistochroneODE,
transcription=dm.GaussLobatto(num_segments=10, compressed=True,
solve_segments='forward'))
phase = traj.add_phase('phase0', phase0)

#
# Set the variables
#
phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

phase.add_state('x', fix_initial=True, fix_final=False)

phase.add_state('y', fix_initial=True, fix_final=False)

phase.add_state('v', fix_initial=True, fix_final=False)

phase.add_control('theta', continuity=True, rate_continuity=True,
units='deg', lower=0.01, upper=179.9)

phase.add_parameter('g', units='m/s**2', val=9.80665)

#
# Minimize time at the end of the phase
#
phase.add_objective('time', loc='final', scaler=10)
phase.add_boundary_constraint('x', loc='final', equals=10)
phase.add_boundary_constraint('y', loc='final', equals=5)

# ---------------------------------------------
# change linear solver for shooting
# ---------------------------------------------
# linear solver
phase0.linear_solver = om.PETScKrylov(assemble_jac=False, iprint=-1)
phase0.linear_solver.precon = om.LinearRunOnce(iprint=-1)

#
# Setup the Problem
#
p.setup()

#
# Set the initial values
#
p['traj0.phase0.t_initial'] = 0.0
p['traj0.phase0.t_duration'] = 2.0

p.set_val('traj0.phase0.states:x', phase.interp('x', ys=[0, 10]))
p.set_val('traj0.phase0.states:y', phase.interp('y', ys=[10, 5]))
p.set_val('traj0.phase0.states:v', phase.interp('v', ys=[0, 9.9]))
p.set_val('traj0.phase0.controls:theta', phase.interp('theta', ys=[5, 100.5]))

#
# Solve for the optimal trajectory
#
dm.run_problem(p, run_driver=True)

self.assert_results(p)
61 changes: 7 additions & 54 deletions dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,10 @@
import unittest

import matplotlib.pyplot as plt
import numpy as np

from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse


plt.switch_backend('Agg')


def initial_guess(t_dur, gam0=np.pi/3):
t = np.linspace(0, t_dur, num=100)
g = -9.80665
Expand Down Expand Up @@ -206,6 +202,12 @@ def compute(self, inputs, outputs):
# converge to the initial point.
phase.set_duration_balance('h', val=0.0)

# In this problem, the default ArmijoGoldsteinLS has issues with extrapolating
# the states and causes the optimization to fail.
# Using the default linesearch or BoundsEnforceLS work better here.
phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, iprint=2)
phase.nonlinear_solver.linesearch = None

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

p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
Expand Down Expand Up @@ -235,60 +237,11 @@ def compute(self, inputs, outputs):
#####################################################
# Run the optimization and final explicit simulation
#####################################################
dm.run_problem(p)
dm.run_problem(p, simulate=True)

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

exp_out = traj.simulate()

#############################################
# Plot the results
#############################################
rad = p.get_val('radius', units='m')[0]
print(f'optimal radius: {rad} m ')
mass = p.get_val('size_comp.mass', units='kg')[0]
print(f'cannonball mass: {mass} kg ')
angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]
print(f'launch angle: {angle} deg')
max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]
print(f'maximum range: {max_range} m')

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))

time_imp = p.get_val('traj.phase.timeseries.time')

time_exp = exp_out.get_val('traj.phase.timeseries.time')

r_imp = p.get_val('traj.phase.timeseries.r')

r_exp = exp_out.get_val('traj.phase.timeseries.r')

h_imp = p.get_val('traj.phase.timeseries.h')

h_exp = exp_out.get_val('traj.phase.timeseries.h')

axes.plot(r_imp, h_imp, 'bo')

axes.plot(r_exp, h_exp, 'b--')

axes.set_xlabel('range (m)')
axes.set_ylabel('altitude (m)')

fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
states = ['r', 'h', 'v', 'gam']
for i, state in enumerate(states):
x_imp = p.get_val(f'traj.phase.timeseries.{state}')

x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')

axes[i].set_ylabel(state)

axes[i].plot(time_imp, x_imp, 'bo')
axes[i].plot(time_exp, x_exp, 'b--')

plt.show()


if __name__ == '__main__': # pragma: no cover
unittest.main()
Loading

0 comments on commit 37ac4e5

Please sign in to comment.