Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated dymos to handle the new OpenMDAO distributed I/O approach #597

Merged
merged 4 commits into from
May 19, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
185 changes: 95 additions & 90 deletions dymos/examples/vanderpol/doc/test_doc_vanderpol.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,9 @@
import os
import unittest
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.style.use('ggplot')

from dymos.utils.doc_utils import save_for_docs
from openmdao.utils.testing_utils import use_tempdirs
from openmdao.utils.mpi import MPI


@use_tempdirs
Expand All @@ -18,116 +15,124 @@ def tearDown(self):

@save_for_docs
def test_vanderpol_for_docs_simulation(self):
from dymos.examples.plotting import plot_results
import dymos as dm
from dymos.examples.vanderpol.vanderpol_dymos import vanderpol

# Create the Dymos problem instance
p = vanderpol(transcription='gauss-lobatto', num_segments=75)

# Run the problem (simulate only)
p.run_model()

# check validity by using scipy.integrate.solve_ivp to integrate the solution
exp_out = p.model.traj.simulate()

# Display the results
plot_results([('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x1',
'time (s)',
'x1 (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x0',
'time (s)',
'x0 (V/s)'),
('traj.phase0.timeseries.states:x0',
'traj.phase0.timeseries.states:x1',
'x0 vs x1',
'x0 vs x1'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.controls:u',
'time (s)',
'control u'),
],
title='Van Der Pol Simulation',
p_sol=p, p_sim=exp_out)

plt.show()
dm.run_problem(p, run_driver=False, simulate=True, make_plots=True)

@save_for_docs
def test_vanderpol_for_docs_optimize(self):
import dymos as dm
from dymos.examples.plotting import plot_results
from dymos.examples.vanderpol.vanderpol_dymos import vanderpol

# Create the Dymos problem instance
p = vanderpol(transcription='gauss-lobatto', num_segments=75,
transcription_order=3, compressed=True, optimizer='SLSQP')

# Find optimal control solution to stop oscillation
dm.run_problem(p)

# check validity by using scipy.integrate.solve_ivp to integrate the solution
exp_out = p.model.traj.simulate()

# Display the results
plot_results([('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x1',
'time (s)',
'x1 (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x0',
'time (s)',
'x0 (V/s)'),
('traj.phase0.timeseries.states:x0',
'traj.phase0.timeseries.states:x1',
'x0 vs x1',
'x0 vs x1'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.controls:u',
'time (s)',
'control u'),
],
title='Van Der Pol Optimization',
p_sol=p, p_sim=exp_out)

plt.show()
dm.run_problem(p, simulate=True, make_plots=True)

@save_for_docs
def test_vanderpol_for_docs_optimize_refine(self):
import dymos as dm
from dymos.examples.plotting import plot_results
from dymos.examples.vanderpol.vanderpol_dymos import vanderpol
from openmdao.utils.assert_utils import assert_near_equal

# Create the Dymos problem instance
p = vanderpol(transcription='gauss-lobatto', num_segments=15,
transcription_order=3, compressed=True, optimizer='SLSQP')

# Enable grid refinement and find optimal control solution to stop oscillation
p.model.traj.phases.phase0.set_refine_options(refine=True)
dm.run_problem(p, refine_iteration_limit=10)

# check validity by using scipy.integrate.solve_ivp to integrate the solution
exp_out = p.model.traj.simulate()

# Display the results
plot_results([('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x1',
'time (s)',
'x1 (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x0',
'time (s)',
'x0 (V/s)'),
('traj.phase0.timeseries.states:x0',
'traj.phase0.timeseries.states:x1',
'x0 vs x1',
'x0 vs x1'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.controls:u',
'time (s)',
'control u'),
],
title='Van Der Pol Optimization with Grid Refinement',
p_sol=p, p_sim=exp_out)

plt.show()

dm.run_problem(p, refine_iteration_limit=10, simulate=True, make_plots=True)

assert_near_equal(p.get_val('traj.phase0.states:x0')[-1, ...], 0.0)
assert_near_equal(p.get_val('traj.phase0.states:x1')[-1, ...], 0.0)
assert_near_equal(p.get_val('traj.phase0.states:J')[-1, ...], 5.2808, tolerance=1.0E-3)
assert_near_equal(p.get_val('traj.phase0.controls:u')[-1, ...], 0.0, tolerance=1.0E-3)


@unittest.skipUnless(MPI, "MPI is required.")
@save_for_docs
@use_tempdirs
class TestVanderpolDelayMPI(unittest.TestCase):
N_PROCS = 2

def test_vanderpol_delay_mpi(self):
import openmdao.api as om
import dymos as dm
from dymos.examples.vanderpol.vanderpol_ode import VanderpolODE
from openmdao.utils.assert_utils import assert_near_equal

DELAY = 0.005

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

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

# define a Trajectory object and add to model
traj = dm.Trajectory()
p.model.add_subsystem('traj', subsys=traj)

t = dm.Radau(num_segments=30, order=3)

# define a Phase as specified above and add to Phase
phase = dm.Phase(ode_class=VanderpolODE, transcription=t,
ode_init_kwargs={'delay': DELAY, 'distrib': True})
traj.add_phase(name='phase0', phase=phase)

t_final = 15
phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=t_final, units='s')

# set the State time options
phase.add_state('x0', fix_initial=False, fix_final=False,
rate_source='x0dot',
units='V/s',
targets='x0') # target required because x0 is an input
phase.add_state('x1', fix_initial=False, fix_final=False,
rate_source='x1dot',
units='V',
targets='x1') # target required because x1 is an input
phase.add_state('J', fix_initial=False, fix_final=False,
rate_source='Jdot',
units=None)

# define the control
phase.add_control(name='u', units=None, lower=-0.75, upper=1.0, continuity=True,
rate_continuity=True,
targets='u') # target required because u is an input

# add constraints
phase.add_boundary_constraint('x0', loc='initial', equals=1.0)
phase.add_boundary_constraint('x1', loc='initial', equals=1.0)
phase.add_boundary_constraint('J', loc='initial', equals=0.0)

phase.add_boundary_constraint('x0', loc='final', equals=0.0)
phase.add_boundary_constraint('x1', loc='final', equals=0.0)

# define objective to minimize
phase.add_objective('J', loc='final')

# setup the problem
p.setup(check=True)

p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = t_final

# add a linearly interpolated initial guess for the state and control curves
p['traj.phase0.states:x0'] = phase.interp('x0', [1, 0])
p['traj.phase0.states:x1'] = phase.interp('x1', [1, 0])
p['traj.phase0.states:J'] = phase.interp('J', [0, 1])
p['traj.phase0.controls:u'] = phase.interp('u', [-0.75, -0.75])

dm.run_problem(p, run_driver=True, simulate=False)

assert_near_equal(p.get_val('traj.phase0.states:x0')[-1, ...], 0.0)
assert_near_equal(p.get_val('traj.phase0.states:x1')[-1, ...], 0.0)
assert_near_equal(p.get_val('traj.phase0.states:J')[-1, ...], 5.2808, tolerance=1.0E-3)
assert_near_equal(p.get_val('traj.phase0.controls:u')[-1, ...], 0.0, tolerance=1.0E-3)
4 changes: 2 additions & 2 deletions dymos/examples/vanderpol/test/test_vanderpol.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,8 +54,8 @@ def test_vanderpol_optimal_mpi(self):
dymos/examples/vanderpol/test/test_vanderpol.py TestVanderpolExampleMPI.test_vanderpol_optimal_mpi
(using varying values for n should give the same answer)
"""
p = vanderpol(transcription='gauss-lobatto', num_segments=75, delay=True,
use_pyoptsparse=True, optimizer='IPOPT')
p = vanderpol(transcription='gauss-lobatto', num_segments=75, delay=0.005,
use_pyoptsparse=True, distrib=True, optimizer='IPOPT')
p.run_driver() # find optimal control solution to stop oscillation

print('Objective function minimized to', p.get_val('traj.phase0.states:J')[-1, ...])
Expand Down
40 changes: 15 additions & 25 deletions dymos/examples/vanderpol/vanderpol_dymos.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
import openmdao.api as om
import dymos as dm
from dymos.examples.vanderpol.vanderpol_ode import vanderpol_ode, vanderpol_ode_group
from dymos.examples.vanderpol.vanderpol_ode import VanderpolODE


def vanderpol(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
compressed=True, optimizer='SLSQP', use_pyoptsparse=False, delay=None):
def vanderpol(transcription='gauss-lobatto', num_segments=40, transcription_order=3,
compressed=True, optimizer='SLSQP', use_pyoptsparse=False, delay=0.0, distrib=True,
solve_segments=False):
"""Dymos problem definition for optimal control of a Van der Pol oscillator"""

# define the OpenMDAO problem
Expand All @@ -30,20 +31,20 @@ def vanderpol(transcription='gauss-lobatto', num_segments=8, transcription_order
if transcription == 'gauss-lobatto':
t = dm.GaussLobatto(num_segments=num_segments,
order=transcription_order,
compressed=compressed)
compressed=compressed,
solve_segments=solve_segments)
elif transcription == 'radau-ps':
t = dm.Radau(num_segments=num_segments,
order=transcription_order,
compressed=compressed)
compressed=compressed,
solve_segments=solve_segments)

# define a Phase as specified above and add to Phase
if not delay:
phase = dm.Phase(ode_class=vanderpol_ode, transcription=t)
else:
phase = dm.Phase(ode_class=vanderpol_ode_group, transcription=t) # distributed component group
phase = dm.Phase(ode_class=VanderpolODE, transcription=t,
ode_init_kwargs={'delay': delay, 'distrib': distrib})
traj.add_phase(name='phase0', phase=phase)

t_final = 15.0
t_final = 15
phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=t_final, units='s')

# set the State time options
Expand Down Expand Up @@ -77,11 +78,6 @@ def vanderpol(transcription='gauss-lobatto', num_segments=8, transcription_order
# setup the problem
p.setup(check=True)

# TODO - Dymos API will soon provide a way to specify this.
# the linear solver used to compute derivatives is not working on MPI, so switch to LinearRunOnce
for phase in traj._phases.values():
phase.linear_solver = om.LinearRunOnce()

p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = t_final

Expand All @@ -91,18 +87,12 @@ def vanderpol(transcription='gauss-lobatto', num_segments=8, transcription_order
p['traj.phase0.states:J'] = phase.interp('J', [0, 1])
p['traj.phase0.controls:u'] = phase.interp('u', [-0.75, -0.75])

p.final_setup()

# debugging helpers:
# om.n2(p) # show n2 diagram
#
# with np.printoptions(linewidth=1024): # display partials for manual checking
# p.check_partials(compact_print=True)

return p


if __name__ == '__main__':
# just set up the problem, test it elsewhere
p = vanderpol(transcription='gauss-lobatto', num_segments=75, transcription_order=3,
compressed=True, optimizer='SLSQP')
p = vanderpol(transcription='radau-ps', num_segments=30, transcription_order=3,
compressed=True, optimizer='SLSQP', delay=0.005, distrib=True, use_pyoptsparse=True)

dm.run_problem(p, run_driver=True, simulate=False)
Loading