From c8a2be800d42deda19e9720fedbdc384ca26a34a Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Fri, 14 May 2021 16:33:17 -0400 Subject: [PATCH 1/4] Updated vanderpol MPI case to follow the new OpenMDAO way of doing serial/distributed I/O rather than entire components. Trajectory no longer assigns solvers to phases. Phases will now only assign a solver if it is needed by solve_segments and if the exisiting one is NonlinearRunOnce or LinearRunOnce. This should allow users to override the default solver before setup, but also make turning solve_segments on and off easy in most cases. --- .../vanderpol/doc/test_doc_vanderpol.py | 173 ++++++++------- dymos/examples/vanderpol/vanderpol_dymos.py | 40 ++-- dymos/examples/vanderpol/vanderpol_ode.py | 199 +++--------------- dymos/trajectory/trajectory.py | 4 +- .../pseudospectral/pseudospectral_base.py | 15 +- 5 files changed, 137 insertions(+), 294 deletions(-) diff --git a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py index 5eab79c14..b4c95a7dc 100644 --- a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py +++ b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py @@ -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 @@ -18,84 +15,28 @@ 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 # Create the Dymos problem instance @@ -104,30 +45,82 @@ def test_vanderpol_for_docs_optimize_refine(self): # 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) + + +@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 + + 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) diff --git a/dymos/examples/vanderpol/vanderpol_dymos.py b/dymos/examples/vanderpol/vanderpol_dymos.py index 6263b3326..c9acfa75a 100644 --- a/dymos/examples/vanderpol/vanderpol_dymos.py +++ b/dymos/examples/vanderpol/vanderpol_dymos.py @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/dymos/examples/vanderpol/vanderpol_ode.py b/dymos/examples/vanderpol/vanderpol_ode.py index fc84663a9..791912ac2 100644 --- a/dymos/examples/vanderpol/vanderpol_ode.py +++ b/dymos/examples/vanderpol/vanderpol_ode.py @@ -4,137 +4,7 @@ from openmdao.utils.array_utils import evenly_distrib_idxs -class vanderpol_ode(om.ExplicitComponent): - """ODE for optimal control of a Van der Pol oscillator - - objective J: - minimize integral of (x0**2 + x1**2 + u**2) for 0.0 <= t <= 15 - - subject to: - x0dot = (1 - x1^2) * x0 - x1 + u - x1dot = x0 - -0.75 <= u <= 1.0 - - initial conditions: - x0(0) = 1.0 x1(0) = 1.0 - - final conditions: - x0(15) = 0.0 x1(15) = 0.0 - """ - - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - # inputs: 2 states and a control - self.add_input('x0', val=np.ones(nn), desc='derivative of Output', units='V/s') - self.add_input('x1', val=np.ones(nn), desc='Output', units='V') - self.add_input('u', val=np.ones(nn), desc='control', units=None) - - # outputs: derivative of states - # the objective function will be treated as a state for computation, so its derivative is an output - self.add_output('x0dot', val=np.ones(nn), desc='second derivative of Output', units='V/s**2') - self.add_output('x1dot', val=np.ones(nn), desc='derivative of Output', units='V/s') - self.add_output('Jdot', val=np.ones(nn), desc='derivative of objective', units='1.0/s') - - # partials - r = c = np.arange(nn) - - self.declare_partials(of='x0dot', wrt='x0', rows=r, cols=c) - self.declare_partials(of='x0dot', wrt='x1', rows=r, cols=c) - self.declare_partials(of='x0dot', wrt='u', rows=r, cols=c, val=1.0) - - self.declare_partials(of='x1dot', wrt='x0', rows=r, cols=c, val=1.0) - self.declare_partials(of='x1dot', wrt='x1', rows=r, cols=c, val=0.0) - self.declare_partials(of='x1dot', wrt='u', rows=r, cols=c, val=0.0) - - self.declare_partials(of='Jdot', wrt='x0', rows=r, cols=c) - self.declare_partials(of='Jdot', wrt='x1', rows=r, cols=c) - self.declare_partials(of='Jdot', wrt='u', rows=r, cols=c) - - def compute(self, inputs, outputs): - x0 = inputs['x0'] - x1 = inputs['x1'] - u = inputs['u'] - - outputs['x0dot'] = (1.0 - x1**2) * x0 - x1 + u - outputs['x1dot'] = x0 - outputs['Jdot'] = x0**2 + x1**2 + u**2 - - def compute_partials(self, inputs, jacobian): - # partials declared with 'val' above do not need to be computed - x0 = inputs['x0'] - x1 = inputs['x1'] - u = inputs['u'] - - jacobian['x0dot', 'x0'] = 1.0 - x1 * x1 - jacobian['x0dot', 'x1'] = -2.0 * x1 * x0 - 1.0 - - jacobian['Jdot', 'x0'] = 2.0 * x0 - jacobian['Jdot', 'x1'] = 2.0 * x1 - jacobian['Jdot', 'u'] = 2.0 * u - - -class vanderpol_ode_group(om.Group): - """Group containing distributed vanderpol_ode pass through and calculation""" - - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - self.add_subsystem(name='vanderpol_ode_passthrough', - subsys=vanderpol_ode_passthrough(num_nodes=nn), - promotes_inputs=['x0', 'x1', 'u']) - - self.add_subsystem(name='vanderpol_ode_delay', - subsys=vanderpol_ode_delay(num_nodes=nn), - promotes_outputs=['x0dot', 'x1dot', 'Jdot']) - - # connect collect_comp (pass through) output to distributed ODE input - self.connect('vanderpol_ode_passthrough.x0pass', 'vanderpol_ode_delay.x0') - self.connect('vanderpol_ode_passthrough.x1pass', 'vanderpol_ode_delay.x1') - self.connect('vanderpol_ode_passthrough.upass', 'vanderpol_ode_delay.u') - - -class vanderpol_ode_passthrough(om.ExplicitComponent): - """Pass through component that just copies control and state from input to output - - if you just use a plain old passthrough (non-distributed with a full sized input and output) component to connect - to the distributed output, the framework will do the MPI allgathering for you""" - - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] # total number of inputs and outputs over all processes - - # inputs: 2 states and a control - self.add_input('x0', val=np.ones(nn), desc='derivative of Output', units='V/s') - self.add_input('x1', val=np.ones(nn), desc='Output', units='V') - self.add_input('u', val=np.ones(nn), desc='control', units=None) - - # outputs: same as inputs - self.add_output('x0pass', val=np.ones(nn), desc='derivative of Output', units='V/s') - self.add_output('x1pass', val=np.ones(nn), desc='Output', units='V') - self.add_output('upass', val=np.ones(nn), desc='control', units=None) - - # partials - row_col = np.arange(nn) - self.declare_partials(of='x0pass', wrt='x0', rows=row_col, cols=row_col, val=1.0) - self.declare_partials(of='x1pass', wrt='x1', rows=row_col, cols=row_col, val=1.0) - self.declare_partials(of='upass', wrt='u', rows=row_col, cols=row_col, val=1.0) - - def compute(self, inputs, outputs): - outputs['x0pass'] = inputs['x0'] - outputs['x1pass'] = inputs['x1'] - outputs['upass'] = inputs['u'] - - -class vanderpol_ode_delay(om.ExplicitComponent): +class VanderpolODE(om.ExplicitComponent): """intentionally slow version of vanderpol_ode for effects of demonstrating distributed component calculations MPI can run this component in multiple processes, distributing the calculation of derivatives. @@ -142,14 +12,13 @@ class vanderpol_ode_delay(om.ExplicitComponent): """ def __init__(self, *args, **kwargs): - self.delay_time = 0.005 self.progress_prints = False super().__init__(*args, **kwargs) def initialize(self): self.options.declare('num_nodes', types=int) - self.options['distributed'] = True - self.options.declare('size', types=int, default=1, desc="Size of input and output vectors.") + self.options.declare('distrib', types=bool, default=False) + self.options.declare('delay', types=(float,), default=0.0) def setup(self): nn = self.options['num_nodes'] @@ -157,68 +26,58 @@ def setup(self): rank = comm.rank sizes, offsets = evenly_distrib_idxs(comm.size, nn) # (#cpus, #inputs) -> (size array, offset array) - start = offsets[rank] + self.start_idx = offsets[rank] self.io_size = sizes[rank] # number of inputs and outputs managed by this distributed process - end = start + self.io_size - - if self.progress_prints: - print('in vanderpol_ode_delay.setup', self.io_size, self.comm.rank) + self.end_idx = self.start_idx + self.io_size # inputs: 2 states and a control - self.add_input('x0', val=np.ones(self.io_size), desc='derivative of Output', units='V/s', - src_indices=np.arange(start, end, dtype=int)) - self.add_input('x1', val=np.ones(self.io_size), desc='Output', units='V', - src_indices=np.arange(start, end, dtype=int)) - self.add_input('u', val=np.ones(self.io_size), desc='control', units=None, - src_indices=np.arange(start, end, dtype=int)) + self.add_input('x0', val=np.ones(nn), desc='derivative of Output', units='V/s') + self.add_input('x1', val=np.ones(nn), desc='Output', units='V') + self.add_input('u', val=np.ones(nn), desc='control', units=None) # outputs: derivative of states # the objective function will be treated as a state for computation, so its derivative is an output - self.add_output('x0dot', val=np.ones(self.io_size), desc='second derivative of Output', units='V/s**2') - self.add_output('x1dot', val=np.ones(self.io_size), desc='derivative of Output', units='V/s') - self.add_output('Jdot', val=np.ones(self.io_size), desc='derivative of objective', units='1.0/s') - - # partials - r = c = np.arange(self.io_size) + self.add_output('x0dot', val=np.ones(self.io_size), desc='second derivative of Output', + units='V/s**2', distributed=self.options['distrib']) + self.add_output('x1dot', val=np.ones(self.io_size), desc='derivative of Output', + units='V/s', distributed=self.options['distrib']) + self.add_output('Jdot', val=np.ones(self.io_size), desc='derivative of objective', + units='1.0/s', distributed=self.options['distrib']) + + # self.declare_coloring(method='cs') + # # partials + r = np.arange(self.io_size, dtype=int) + c = r + self.start_idx self.declare_partials(of='x0dot', wrt='x0', rows=r, cols=c) self.declare_partials(of='x0dot', wrt='x1', rows=r, cols=c) self.declare_partials(of='x0dot', wrt='u', rows=r, cols=c, val=1.0) self.declare_partials(of='x1dot', wrt='x0', rows=r, cols=c, val=1.0) - self.declare_partials(of='x1dot', wrt='x1', rows=r, cols=c, val=0.0) - self.declare_partials(of='x1dot', wrt='u', rows=r, cols=c, val=0.0) self.declare_partials(of='Jdot', wrt='x0', rows=r, cols=c) self.declare_partials(of='Jdot', wrt='x1', rows=r, cols=c) self.declare_partials(of='Jdot', wrt='u', rows=r, cols=c) def compute(self, inputs, outputs): - if self.progress_prints: - sizes = (len(inputs['x0']), len(inputs['x1']), len(inputs['u'])) - print('in vanderpol_ode_delay.compute', sizes, self.comm.rank) + # introduce slowness proportional to size of computation + time.sleep(self.options['delay'] * self.io_size) - time.sleep(self.delay_time * self.io_size) # introduce slowness proportional to size of computation - - x0 = inputs['x0'] - x1 = inputs['x1'] - u = inputs['u'] + # The inputs contain the entire vector, be each rank will only operate on a portion of it. + x0 = inputs['x0'][self.start_idx:self.end_idx] + x1 = inputs['x1'][self.start_idx:self.end_idx] + u = inputs['u'][self.start_idx:self.end_idx] outputs['x0dot'] = (1.0 - x1**2) * x0 - x1 + u outputs['x1dot'] = x0 outputs['Jdot'] = x0**2 + x1**2 + u**2 def compute_partials(self, inputs, jacobian): - if self.progress_prints: - sizes = (len(inputs['x0']), len(inputs['x1']), len(inputs['u'])) - print('in vanderpol_ode_delay.compute_partials', sizes) - - time.sleep(self.delay_time * self.io_size) # introduce slowness proportional to size of computation + time.sleep(self.options['delay'] * self.io_size) - # partials declared with 'val' above do not need to be computed - x0 = inputs['x0'] - x1 = inputs['x1'] - u = inputs['u'] + x0 = inputs['x0'][self.start_idx:self.end_idx] + x1 = inputs['x1'][self.start_idx:self.end_idx] + u = inputs['u'][self.start_idx:self.end_idx] jacobian['x0dot', 'x0'] = 1.0 - x1 * x1 jacobian['x0dot', 'x1'] = -2.0 * x1 * x0 - 1.0 diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py index 202c07550..1d1f2d532 100644 --- a/dymos/trajectory/trajectory.py +++ b/dymos/trajectory/trajectory.py @@ -267,9 +267,7 @@ def setup(self): phases_group = self.add_subsystem('phases', subsys=om.ParallelGroup()) for name, phs in self._phases.items(): - g = phases_group.add_subsystem(name, phs, **self._phase_add_kwargs[name]) - # DirectSolvers were moved down into the phases for use with MPI - g.linear_solver = om.DirectSolver() + phases_group.add_subsystem(name, phs, **self._phase_add_kwargs[name]) if self._linkages: self._setup_linkages() diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py index 514d20f09..dd71275a4 100644 --- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py +++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py @@ -441,12 +441,15 @@ def configure_solvers(self, phase): The phase object to which this transcription instance applies. """ if self.any_solved_segs or self.any_connected_opt_segs: - newton = phase.nonlinear_solver = om.NewtonSolver() - newton.options['solve_subsystems'] = True - newton.options['maxiter'] = 100 - newton.options['iprint'] = -1 - newton.linesearch = om.BoundsEnforceLS() - phase.linear_solver = om.DirectSolver() + # Only override the solvers if the user hasn't set them to something else. + if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce): + newton = phase.nonlinear_solver = om.NewtonSolver() + newton.options['solve_subsystems'] = True + newton.options['maxiter'] = 100 + newton.options['iprint'] = 2 + newton.linesearch = om.BoundsEnforceLS() + if isinstance(phase.linear_solver, om.LinearRunOnce): + phase.linear_solver = om.DirectSolver() def setup_timeseries_outputs(self, phase): """ From c9111fca978261d97932c32b09a127288dd54667 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Fri, 14 May 2021 17:59:28 -0400 Subject: [PATCH 2/4] fixed existing vanderpol optimization mpi test --- dymos/examples/vanderpol/test/test_vanderpol.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dymos/examples/vanderpol/test/test_vanderpol.py b/dymos/examples/vanderpol/test/test_vanderpol.py index 55cf26de7..f85474d44 100644 --- a/dymos/examples/vanderpol/test/test_vanderpol.py +++ b/dymos/examples/vanderpol/test/test_vanderpol.py @@ -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, ...]) From 3d6937e3e1f869b21d433fb2358cdf39cae97b81 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 17 May 2021 08:15:38 -0400 Subject: [PATCH 3/4] added some asserts to the vanderpol doc tests. removed unneccessary prints from the timeseries plots. --- dymos/examples/vanderpol/doc/test_doc_vanderpol.py | 12 ++++++++++++ dymos/visualization/timeseries_plots.py | 4 ++-- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py index b4c95a7dc..2ba749dfc 100644 --- a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py +++ b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py @@ -38,6 +38,7 @@ def test_vanderpol_for_docs_optimize(self): def test_vanderpol_for_docs_optimize_refine(self): import dymos as dm 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, @@ -48,6 +49,11 @@ def test_vanderpol_for_docs_optimize_refine(self): 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 @@ -59,6 +65,7 @@ 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 @@ -124,3 +131,8 @@ def test_vanderpol_delay_mpi(self): 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) diff --git a/dymos/visualization/timeseries_plots.py b/dymos/visualization/timeseries_plots.py index 398b8f95d..2da58b497 100644 --- a/dymos/visualization/timeseries_plots.py +++ b/dymos/visualization/timeseries_plots.py @@ -295,13 +295,13 @@ def timeseries_plots(solution_recorder_filename, simulation_record_file=None, pl cr = om.CaseReader(solution_recorder_filename) # get outputs from the solution - solution_cases = cr.list_cases('problem') + solution_cases = cr.list_cases('problem', out_stream=None) last_solution_case = cr.get_case(solution_cases[-1]) # If plotting simulation results, get the values for those variables if simulation_record_file: cr_simulate = om.CaseReader(simulation_record_file) - system_simulation_cases = cr_simulate.list_cases('problem') + system_simulation_cases = cr_simulate.list_cases('problem', out_stream=None) last_simulation_case = cr_simulate.get_case(system_simulation_cases[-1]) else: last_simulation_case = None From 1273551423df1f2e40d469ed10899b3cf421dc02 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 17 May 2021 09:40:24 -0400 Subject: [PATCH 4/4] cleaned up numerous unnecessary imports in run_problem --- dymos/run_problem.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/dymos/run_problem.py b/dymos/run_problem.py index ae51bf09b..60eced5de 100755 --- a/dymos/run_problem.py +++ b/dymos/run_problem.py @@ -1,20 +1,11 @@ import warnings -from .grid_refinement.ph_adaptive.ph_adaptive import PHAdaptive -from .grid_refinement.hp_adaptive.hp_adaptive import HPAdaptive -from .grid_refinement.write_iteration import write_error, write_refine_iter -from .grid_refinement.refinement import _refine_iter -from .phase.phase import Phase - import openmdao.api as om -import dymos as dm -import numpy as np from dymos.trajectory.trajectory import Trajectory -from dymos.load_case import load_case, find_phases +from dymos.load_case import load_case from dymos.visualization.timeseries_plots import timeseries_plots -from dymos.grid_refinement.error_estimation import check_error -import os -import sys + +from .grid_refinement.refinement import _refine_iter def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_driver=True,