diff --git a/benchmark/benchmark_balanced_field.py b/benchmark/benchmark_balanced_field.py index cfc4b2f26..332103315 100644 --- a/benchmark/benchmark_balanced_field.py +++ b/benchmark/benchmark_balanced_field.py @@ -20,7 +20,7 @@ def _run_balanced_field_length_problem(tx=dm.GaussLobatto, timeseries=True, sim= p.driver.options['optimizer'] = optimizer p.driver.options['print_results'] = False if optimizer == 'IPOPT': - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 # First Phase: Brake release to V1 - both engines operable br_to_v1 = dm.Phase(ode_class=BalancedFieldODEComp, diff --git a/benchmark/benchmark_racecar.py b/benchmark/benchmark_racecar.py index d30fe1124..b1b5636c4 100644 --- a/benchmark/benchmark_racecar.py +++ b/benchmark/benchmark_racecar.py @@ -129,7 +129,7 @@ def _run_racecar_problem(transcription, timeseries=False): p.driver.opt_settings['compl_inf_tol'] = 1e-3 p.driver.opt_settings['acceptable_iter'] = 0 p.driver.opt_settings['tol'] = 1e-3 - p.driver.opt_settings['print_level'] = 5 + 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' diff --git a/docs/dymos_book/faq/downstream_analysis.ipynb b/docs/dymos_book/faq/downstream_analysis.ipynb index ca515d5a4..7c0dae322 100644 --- a/docs/dymos_book/faq/downstream_analysis.ipynb +++ b/docs/dymos_book/faq/downstream_analysis.ipynb @@ -47,7 +47,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# How do connect the outputs of Dymos to a downstream analysis?\n", + "# How do I connect the outputs of Dymos to a downstream analysis?\n", "\n", "One of the design goals of Dymos is to allow the trajectory to be a part of a larger multidisciplinary optimization problem.\n", "Sometimes, you may want to take the results from the Dymos trajectory and feed them to some downstream analysis.\n", diff --git a/docs/dymos_book/faq/upstream_analysis.ipynb b/docs/dymos_book/faq/upstream_analysis.ipynb index 9a50d9545..204fe61ce 100644 --- a/docs/dymos_book/faq/upstream_analysis.ipynb +++ b/docs/dymos_book/faq/upstream_analysis.ipynb @@ -47,7 +47,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# How do connect the outputs of an upstream analysis as inputs to Dymos?\n", + "# How do I connect the outputs of an upstream analysis as inputs to Dymos?\n", "\n", "One of the design goals of Dymos is to allow the trajectory to be a part of a larger multidisciplinary optimization problem.\n", "Sometimes, you may want to take the results of some upstream design and use them as an input to a Dymos trajectory.\n", diff --git a/docs/dymos_book/features/trajectories/trajectories.ipynb b/docs/dymos_book/features/trajectories/trajectories.ipynb index 73add00b6..5fd96cda9 100644 --- a/docs/dymos_book/features/trajectories/trajectories.ipynb +++ b/docs/dymos_book/features/trajectories/trajectories.ipynb @@ -224,7 +224,7 @@ "If we implemented these as parameters within each phase individually, we would need some constraints to ensure that they held the same value within each phase.\n", "To avoid this complexity, Dymos Trajectory objects support their own Parameters.\n", "\n", - "Like their Phase-based counterparts, Trajectory parameters produce may be design variables for the problem or used as inputs to the trajectory from external sources.\n", + "Like their Phase-based counterparts, Trajectory parameters may be design variables for the problem or they may be inputs to the trajectory from external sources.\n", "\n", "When using Trajectory parameters, their values are connected to each phase as an Input Parameter within the Phase.\n", "Because ODEs in different phases may have different names for parameters (e.g. 'mass', 'm', 'm_total', etc) Dymos allows the user to specify the targeted ODE parameters on a phase-by-phase basis using the `targets` and `target_params` option.\n", diff --git a/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb b/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb index bdf9ac319..d1c35b618 100644 --- a/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb +++ b/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb @@ -156,7 +156,7 @@ "\n", "## How does the number and order of segments affect the solution?\n", "\n", - "Obviously, a single third-order polynomial won't be able to fit highly oscillator behavior.\n", + "Obviously, a single third-order polynomial won't be able to fit highly oscillatory behavior.\n", "In this case, our guess of using four segments (equally spaced) in the phase wasn't quite sufficient.\n", "Let's try increasing that number to ten third-order segments." ] diff --git a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py index 3eb1f6968..37d55b62a 100644 --- a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py +++ b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py @@ -29,7 +29,7 @@ def test_balanced_field_length_for_docs(self): p.driver.options['optimizer'] = optimizer p.driver.options['print_results'] = False if optimizer == 'IPOPT': - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['derivative_test'] = 'first-order' # First Phase: Brake release to V1 - both engines operable diff --git a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py index 2811178ec..04aeb4a6f 100644 --- a/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py +++ b/dymos/examples/balanced_field/test/test_balanced_field_func_comp.py @@ -158,7 +158,7 @@ def _run_problem(self, tx): p.driver.options['optimizer'] = 'IPOPT' p.driver.options['print_results'] = True - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['mu_strategy'] = 'adaptive' p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' diff --git a/dymos/examples/balanced_field/test/test_balanced_field_length.py b/dymos/examples/balanced_field/test/test_balanced_field_length.py index 6a23e3010..7e373c76a 100644 --- a/dymos/examples/balanced_field/test/test_balanced_field_length.py +++ b/dymos/examples/balanced_field/test/test_balanced_field_length.py @@ -18,7 +18,7 @@ def _make_problem(self): p.driver = om.pyOptSparseDriver() p.driver.options['optimizer'] = 'IPOPT' - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['derivative_test'] = 'first-order' p.driver.declare_coloring() diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py index 6c13227ab..c3dbc8c14 100644 --- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py +++ b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py @@ -310,7 +310,7 @@ def test_brachistochrone_for_docs_coloring_demo_solve_segments(self): # p = om.Problem(model=om.Group()) p.driver = om.pyOptSparseDriver(optimizer='IPOPT') - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 # p.driver.declare_coloring() # diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py index 275423a57..c38057867 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py @@ -27,7 +27,7 @@ def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_o elif optimizer == 'IPOPT': p.driver.opt_settings['mu_init'] = 1e-3 p.driver.opt_settings['max_iter'] = 500 - p.driver.opt_settings['print_level'] = 5 + 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' diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py index a5eaf23c1..8b363b183 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py @@ -1,308 +1,312 @@ """Unit Tests for the code that does automatic report generation""" import unittest import pathlib -import sys import os -from io import StringIO +from packaging.version import Version import openmdao.api as om -from openmdao.test_suite.components.paraboloid import Paraboloid -from openmdao.test_suite.components.sellar_feature import SellarMDA import openmdao.core.problem -from openmdao.core.constants import _UNDEFINED -from openmdao.utils.assert_utils import assert_warning -from openmdao.utils.general_utils import set_pyoptsparse_opt -from openmdao.utils.reports_system import set_default_reports_dir, _reports_dir, register_report, \ - list_reports, clear_reports, run_n2_report, setup_default_reports, report_function +import openmdao.utils.reports_system as reports_system from openmdao.utils.testing_utils import use_tempdirs -from openmdao.utils.mpi import MPI from openmdao.utils.tests.test_hooks import hooks_active from openmdao.visualization.n2_viewer.n2_viewer import _default_n2_filename from openmdao.visualization.scaling_viewer.scaling_report import _default_scaling_filename +from openmdao import __version__ as openmdao_version -try: - from openmdao.vectors.petsc_vector import PETScVector -except ImportError: - PETScVector = None - -OPT, OPTIMIZER = set_pyoptsparse_opt('SLSQP') - -if OPTIMIZER: - from openmdao.drivers.pyoptsparse_driver import pyOptSparseDriver import dymos as dm from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE -@use_tempdirs -class TestSubproblemReportToggle(unittest.TestCase): +def setup_model_radau(do_reports): + p = om.Problem(model=om.Group()) - def setUp(self): - self.n2_filename = _default_n2_filename - self.scaling_filename = _default_scaling_filename + p.driver = om.ScipyOptimizeDriver() + p.driver.declare_coloring(tol=1.0E-12) - # set things to a known initial state for all the test runs - openmdao.core.problem._problem_names = [] # need to reset these to simulate separate runs - os.environ.pop('OPENMDAO_REPORTS', None) - os.environ.pop('OPENMDAO_REPORTS_DIR', None) - # We need to remove the TESTFLO_RUNNING environment variable for these tests to run. - # The reports code checks to see if TESTFLO_RUNNING is set and will not do anything if set - # But we need to remember whether it was set so we can restore it - self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None) - clear_reports() - set_default_reports_dir(_reports_dir) + t = dm.Radau(num_segments=10, order=3) - self.count = 0 + traj = dm.Trajectory() + phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t) + p.model.add_subsystem('traj0', traj) + traj.add_phase('phase0', phase) - def tearDown(self): - # restore what was there before running the test - if self.testflo_running is not None: - os.environ['TESTFLO_RUNNING'] = self.testflo_running + # p.model.add_subsystem('traj0', traj) - @hooks_active - def test_no_sim_reports(self): - setup_default_reports() + phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - p = om.Problem(model=om.Group()) + phase.add_state('x', fix_initial=True, fix_final=False) + phase.add_state('y', fix_initial=True, fix_final=False) - p.driver = om.ScipyOptimizeDriver() - p.driver.declare_coloring(tol=1.0E-12) + # Note that by omitting the targets here Dymos will automatically attempt to connect + # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found. + phase.add_state('v', fix_initial=True, fix_final=False) - t = dm.Radau(num_segments=10, order=3) + phase.add_control('theta', continuity=True, rate_continuity=True, + units='deg', lower=0.01, upper=179.9) - traj = dm.Trajectory() - phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t) - p.model.add_subsystem('traj0', traj) - traj.add_phase('phase0', phase) + phase.add_parameter('g', targets=['g'], units='m/s**2') - # p.model.add_subsystem('traj0', traj) + phase.add_boundary_constraint('x', loc='final', equals=10) + phase.add_boundary_constraint('y', loc='final', equals=5) + # Minimize time at the end of the phase + phase.add_objective('time_phase', loc='final', scaler=10) - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) + p.setup() - phase.add_state('x', fix_initial=True, fix_final=False) - phase.add_state('y', fix_initial=True, fix_final=False) + phase.set_simulate_options(method='RK23') - # Note that by omitting the targets here Dymos will automatically attempt to connect - # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found. - phase.add_state('v', fix_initial=True, fix_final=False) + p['traj0.phase0.t_initial'] = 0.0 + p['traj0.phase0.t_duration'] = 2.0 - phase.add_control('theta', - continuity=True, rate_continuity=True, - units='deg', lower=0.01, upper=179.9) + p['traj0.phase0.states:x'] = phase.interp('x', [0, 10]) + p['traj0.phase0.states:y'] = phase.interp('y', [10, 5]) + p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9]) + p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100]) + p['traj0.phase0.parameters:g'] = 9.80665 - phase.add_parameter('g', targets=['g'], units='m/s**2') + if do_reports: + dm.run_problem(p, run_driver=True, simulate=True, simulate_kwargs={'reports': True}) + else: + dm.run_problem(p, run_driver=True, simulate=True) - phase.add_boundary_constraint('x', loc='final', equals=10) - phase.add_boundary_constraint('y', loc='final', equals=5) - # Minimize time at the end of the phase - phase.add_objective('time_phase', loc='final', scaler=10) + return p - p.setup() - phase.set_simulate_options(method='RK23') +def setup_model_shooting(do_reports): + prob = om.Problem() - p['traj0.phase0.t_initial'] = 0.0 - p['traj0.phase0.t_duration'] = 2.0 + prob.driver = om.ScipyOptimizeDriver() + prob.driver.declare_coloring(tol=1.0E-12) - p['traj0.phase0.states:x'] = phase.interp('x', [0, 10]) - p['traj0.phase0.states:y'] = phase.interp('y', [10, 5]) - p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9]) - p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100]) - p['traj0.phase0.parameters:g'] = 9.80665 + tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto', + method='rk4', order=5, + num_steps_per_segment=5, + compressed=False, + subprob_reports=do_reports) - dm.run_problem(p, run_driver=True, simulate=True) + phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) - problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name) - report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] + phase.set_time_options(units='s', fix_initial=True, duration_bounds=(1.0, 10.0)) - # Test that a report subdir was made - self.assertEqual(len(report_subdirs), 1) + # automatically discover states + phase.set_state_options('x', fix_initial=True) + phase.set_state_options('y', fix_initial=True) + phase.set_state_options('v', fix_initial=True) - path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename) - self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') - path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename) - self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') + phase.add_parameter('g', val=1.0, units='m/s**2', opt=True, lower=1, upper=9.80665) + phase.add_control('theta', val=45.0, units='deg', opt=True, lower=1.0E-6, upper=179.9, + ref=90., rate2_continuity=True) - @hooks_active - def test_make_sim_reports(self): - setup_default_reports() + phase.add_boundary_constraint('x', loc='final', equals=10.0) + phase.add_boundary_constraint('y', loc='final', equals=5.0) - p = om.Problem(model=om.Group()) + prob.model.add_subsystem('phase0', phase) - p.driver = om.ScipyOptimizeDriver() - p.driver.declare_coloring(tol=1.0E-12) + phase.add_objective('time', loc='final') - t = dm.Radau(num_segments=10, order=3) + prob.setup(force_alloc_complex=True) - traj = dm.Trajectory() - phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t) - p.model.add_subsystem('traj0', traj) - traj.add_phase('phase0', phase) + prob.set_val('phase0.t_initial', 0.0) + prob.set_val('phase0.t_duration', 2) + prob.set_val('phase0.states:x', 0.0) + prob.set_val('phase0.states:y', 10.0) + prob.set_val('phase0.states:v', 1.0E-6) + prob.set_val('phase0.parameters:g', 1.0, units='m/s**2') + prob.set_val('phase0.controls:theta', phase.interp('theta', ys=[0.01, 90]), units='deg') - # p.model.add_subsystem('traj0', traj) + dm.run_problem(prob, run_driver=True, simulate=False) - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) + return prob - phase.add_state('x', fix_initial=True, fix_final=False) - phase.add_state('y', fix_initial=True, fix_final=False) - # Note that by omitting the targets here Dymos will automatically attempt to connect - # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found. - phase.add_state('v', fix_initial=True, fix_final=False) +# reports API between 3.18 and 3.19, so handle it here in order to be able to test against older +# versions of openmdao +if Version(openmdao_version) > Version("3.18"): + from openmdao.utils.reports_system import get_reports_dir, clear_reports - phase.add_control('theta', - continuity=True, rate_continuity=True, - units='deg', lower=0.01, upper=179.9) + @use_tempdirs + class TestSubproblemReportToggle(unittest.TestCase): - phase.add_parameter('g', targets=['g'], units='m/s**2') + def setUp(self): + self.n2_filename = _default_n2_filename + self.scaling_filename = _default_scaling_filename - phase.add_boundary_constraint('x', loc='final', equals=10) - phase.add_boundary_constraint('y', loc='final', equals=5) - # Minimize time at the end of the phase - phase.add_objective('time_phase', loc='final', scaler=10) + # set things to a known initial state for all the test runs + openmdao.core.problem._problem_names = [] # need to reset these to simulate separate runs + os.environ.pop('OPENMDAO_REPORTS', None) + os.environ.pop('OPENMDAO_REPORTS_DIR', None) + # We need to remove the TESTFLO_RUNNING environment variable for these tests to run. + # The reports code checks to see if TESTFLO_RUNNING is set and will not do anything if set + # But we need to remember whether it was set so we can restore it + self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None) + clear_reports() - p.setup() + self.count = 0 - phase.set_simulate_options(method='RK23') + def tearDown(self): + # restore what was there before running the test + if self.testflo_running is not None: + os.environ['TESTFLO_RUNNING'] = self.testflo_running - p['traj0.phase0.t_initial'] = 0.0 - p['traj0.phase0.t_duration'] = 2.0 + @hooks_active + def test_no_sim_reports(self): + p = setup_model_radau(do_reports=False) - p['traj0.phase0.states:x'] = phase.interp('x', [0, 10]) - p['traj0.phase0.states:y'] = phase.interp('y', [10, 5]) - p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9]) - p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100]) - p['traj0.phase0.parameters:g'] = 9.80665 + report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()] - dm.run_problem(p, run_driver=True, simulate=True, simulate_kwargs={'reports': True}) + # Test that a report subdir was made + self.assertEqual(len(report_subdirs), 1) - problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name) - report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] + path = pathlib.Path(report_subdirs[0]).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + path = pathlib.Path(report_subdirs[0]).joinpath(self.scaling_filename) + self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') - # Test that a report subdir was made - # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation. - self.assertEqual(len(report_subdirs), 12) + @hooks_active + def test_make_sim_reports(self): + p = setup_model_radau(do_reports=True) - for subdir in report_subdirs: - path = pathlib.Path(subdir).joinpath(self.n2_filename) - self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()] - @hooks_active - def test_explicitshooting_no_subprob_reports(self): - setup_default_reports() + # Test that a report subdir was made + # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation. + self.assertEqual(len(report_subdirs), 12) - prob = om.Problem() + for subdir in report_subdirs: + path = pathlib.Path(subdir).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') - prob.driver = om.ScipyOptimizeDriver() - prob.driver.declare_coloring(tol=1.0E-12) + @hooks_active + def test_explicitshooting_no_subprob_reports(self): + p = setup_model_shooting(do_reports=False) - tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto', - method='rk4', order=5, - num_steps_per_segment=5, - compressed=False) + report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()] - phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) + # Test that a report subdir was made + self.assertEqual(len(report_subdirs), 1) - phase.set_time_options(units='s', fix_initial=True, duration_bounds=(1.0, 10.0)) + path = pathlib.Path(report_subdirs[0]).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + path = pathlib.Path(report_subdirs[0]).joinpath(self.scaling_filename) + self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') - # automatically discover states - phase.set_state_options('x', fix_initial=True) - phase.set_state_options('y', fix_initial=True) - phase.set_state_options('v', fix_initial=True) + @hooks_active + def test_explicitshooting_make_subprob_reports(self): + p = setup_model_shooting(do_reports=True) - phase.add_parameter('g', val=1.0, units='m/s**2', opt=True, lower=1, upper=9.80665) - phase.add_control('theta', val=45.0, units='deg', opt=True, lower=1.0E-6, upper=179.9, - ref=90., rate2_continuity=True) + report_subdirs = [e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()] - phase.add_boundary_constraint('x', loc='final', equals=10.0) - phase.add_boundary_constraint('y', loc='final', equals=5.0) + # Test that a report subdir was made + # There is the nominal problem, a subproblem for integration, and a subproblem for the derivatives. + self.assertEqual(len(report_subdirs), 3) - prob.model.add_subsystem('phase0', phase) + path = pathlib.Path(report_subdirs[0]).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + path = pathlib.Path(report_subdirs[0]).joinpath(self.scaling_filename) + self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') - phase.add_objective('time', loc='final') + for subdir in report_subdirs: + path = pathlib.Path(subdir).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') - prob.setup(force_alloc_complex=True) +else: # old OM versions before reports API changed... + from openmdao.utils.reports_system import set_default_reports_dir, _reports_dir, clear_reports, \ + setup_default_reports - prob.set_val('phase0.t_initial', 0.0) - prob.set_val('phase0.t_duration', 2) - prob.set_val('phase0.states:x', 0.0) - prob.set_val('phase0.states:y', 10.0) - prob.set_val('phase0.states:v', 1.0E-6) - prob.set_val('phase0.parameters:g', 1.0, units='m/s**2') - prob.set_val('phase0.controls:theta', phase.interp('theta', ys=[0.01, 90]), units='deg') + @use_tempdirs + class TestSubproblemReportToggle(unittest.TestCase): - dm.run_problem(prob, run_driver=True, simulate=False) + def setUp(self): + self.n2_filename = _default_n2_filename + self.scaling_filename = _default_scaling_filename - problem_reports_dir = pathlib.Path(_reports_dir).joinpath(prob._name) - report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] + # set things to a known initial state for all the test runs + openmdao.core.problem._problem_names = [] # need to reset these to simulate separate runs + os.environ.pop('OPENMDAO_REPORTS', None) + os.environ.pop('OPENMDAO_REPORTS_DIR', None) + # We need to remove the TESTFLO_RUNNING environment variable for these tests to run. + # The reports code checks to see if TESTFLO_RUNNING is set and will not do anything if set + # But we need to remember whether it was set so we can restore it + self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None) + clear_reports() + set_default_reports_dir(_reports_dir) - # Test that a report subdir was made - self.assertEqual(len(report_subdirs), 1) + self.count = 0 - path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename) - self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') - path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename) - self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') + def tearDown(self): + # restore what was there before running the test + if self.testflo_running is not None: + os.environ['TESTFLO_RUNNING'] = self.testflo_running - @hooks_active - def test_explicitshooting_make_subprob_reports(self): - setup_default_reports() + @hooks_active + def test_no_sim_reports(self): + setup_default_reports() - prob = om.Problem() + p = setup_model_radau(do_reports=False) - prob.driver = om.ScipyOptimizeDriver() - prob.driver.declare_coloring(tol=1.0E-12) + problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name) + report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] - tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto', - method='rk4', order=5, - num_steps_per_segment=5, - compressed=False, - subprob_reports=True) + # Test that a report subdir was made + self.assertEqual(len(report_subdirs), 1) - phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) + path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename) + self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') - phase.set_time_options(units='s', fix_initial=True, duration_bounds=(1.0, 10.0)) + @hooks_active + def test_make_sim_reports(self): + setup_default_reports() - # automatically discover states - phase.set_state_options('x', fix_initial=True) - phase.set_state_options('y', fix_initial=True) - phase.set_state_options('v', fix_initial=True) + p = setup_model_radau(do_reports=True) - phase.add_parameter('g', val=1.0, units='m/s**2', opt=True, lower=1, upper=9.80665) - phase.add_control('theta', val=45.0, units='deg', opt=True, lower=1.0E-6, upper=179.9, - ref=90., rate2_continuity=True) + report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] - phase.add_boundary_constraint('x', loc='final', equals=10.0) - phase.add_boundary_constraint('y', loc='final', equals=5.0) + # Test that a report subdir was made + # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation. + self.assertEqual(len(report_subdirs), 12) - prob.model.add_subsystem('phase0', phase) + for subdir in report_subdirs: + path = pathlib.Path(subdir).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') - phase.add_objective('time', loc='final') + @hooks_active + def test_explicitshooting_no_subprob_reports(self): + setup_default_reports() - prob.setup(force_alloc_complex=True) + p = setup_model_shooting(do_reports=False) - prob.set_val('phase0.t_initial', 0.0) - prob.set_val('phase0.t_duration', 2) - prob.set_val('phase0.states:x', 0.0) - prob.set_val('phase0.states:y', 10.0) - prob.set_val('phase0.states:v', 1.0E-6) - prob.set_val('phase0.parameters:g', 1.0, units='m/s**2') - prob.set_val('phase0.controls:theta', phase.interp('theta', ys=[0.01, 90]), units='deg') + problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name) + report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] - dm.run_problem(prob, run_driver=True, simulate=False) + # Test that a report subdir was made + self.assertEqual(len(report_subdirs), 1) - problem_reports_dir = pathlib.Path(_reports_dir).joinpath(prob._name) - report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] + path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename) + self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') - # Test that a report subdir was made - # There is the nominal problem, a subproblem for integration, and a subproblem for the derivatives. - self.assertEqual(len(report_subdirs), 3) + @hooks_active + def test_explicitshooting_make_subprob_reports(self): + setup_default_reports() - path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename) - self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') - path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename) - self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') + p = setup_model_shooting(do_reports=True) - for subdir in report_subdirs: - path = pathlib.Path(subdir).joinpath(self.n2_filename) + problem_reports_dir = pathlib.Path(_reports_dir).joinpath(p._name) + report_subdirs = [e for e in pathlib.Path(_reports_dir).iterdir() if e.is_dir()] + + # Test that a report subdir was made + # There is the nominal problem, a subproblem for integration, and a subproblem for the derivatives. + self.assertEqual(len(report_subdirs), 3) + + path = pathlib.Path(problem_reports_dir).joinpath(self.n2_filename) self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') + path = pathlib.Path(problem_reports_dir).joinpath(self.scaling_filename) + self.assertTrue(path.is_file(), f'The scaling report file, {str(path)}, was not found') + + for subdir in report_subdirs: + path = pathlib.Path(subdir).joinpath(self.n2_filename) + self.assertTrue(path.is_file(), f'The N2 report file, {str(path)} was not found') diff --git a/dymos/examples/brachistochrone/test/test_duplicate_constraints.py b/dymos/examples/brachistochrone/test/test_duplicate_constraints.py index 4e18cba97..4f8e6fcb0 100644 --- a/dymos/examples/brachistochrone/test/test_duplicate_constraints.py +++ b/dymos/examples/brachistochrone/test/test_duplicate_constraints.py @@ -614,4 +614,4 @@ def test_duplicate_timeseries_output(self): warnings.simplefilter("error") phase.add_path_constraint('check', constraint_name='bounded_check', upper=100, lower=-100) - self.assertIn('bounded_check', [op['output_name'] for op in phase._timeseries['timeseries']['outputs']]) + self.assertIn('bounded_check', [op['output_name'] for op in phase._timeseries['timeseries']['outputs'].values()]) diff --git a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py index 3df30fdea..db8879165 100644 --- a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py +++ b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py @@ -75,7 +75,7 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F rate_source='deltav_dot', units='DU/TU') coast.add_parameter('u1', opt=False, val=0.0, units='deg') - coast.add_parameter('c', val=0.0, units='DU/TU') + coast.add_parameter('c', val=1.0, units='DU/TU') # Third Phase (burn) burn2 = dm.Phase(ode_class=FiniteBurnODE, transcription=t[transcription]) @@ -156,6 +156,13 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F traj.link_phases(phases=['burn1', 'burn2'], vars=['accel']) + if connected: + # If running connected and under MPI the phases subsystem requires a Nonlinear Block Jacobi solver. + # This is not the most efficient way to actually solve this problem but it demonstrates access + # to the traj.phases subsystem before setup. + traj.phases.nonlinear_solver = om.NonlinearBlockJac() + traj.phases.linear_solver = om.PETScKrylov() + return traj @@ -255,6 +262,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP p.set_val('traj.coast.states:vr', val=coast.interp('vr', [0.3285, 0])) p.set_val('traj.coast.states:vt', val=coast.interp('vt', [0.97, 1])) p.set_val('traj.coast.states:accel', val=coast.interp('accel', [0, 0])) + p.set_val('traj.burn1.controls:u1', val=burn1.interp('u1', [0, 0])) if burn2 in p.model.traj.phases._subsystems_myproc: if connected: @@ -267,6 +275,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP p.set_val('traj.burn2.states:vt', val=burn2.interp('vt', [np.sqrt(1 / r_target), 1])) p.set_val('traj.burn2.states:deltav', val=burn2.interp('deltav', [0.2, 0.1])) p.set_val('traj.burn2.states:accel', val=burn2.interp('accel', [0., 0.1])) + p.set_val('traj.burn2.controls:u1', val=burn2.interp('u1', [0, 0])) else: p.set_val('traj.burn2.t_initial', val=5.25) @@ -280,6 +289,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP p.set_val('traj.burn2.states:accel', val=burn2.interp('accel', [0.1, 0])) p.set_val('traj.burn2.controls:u1', val=burn2.interp('u1', [0, 0])) - dm.run_problem(p, run_driver=run_driver, simulate=simulate, restart=restart) + if run_driver or simulate: + dm.run_problem(p, run_driver=run_driver, simulate=simulate, restart=restart) return p diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py index 6f1f4ca39..516ec16c9 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py @@ -129,7 +129,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP elif optimizer == 'IPOPT': p.driver.opt_settings['hessian_approximation'] = 'limited-memory' p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['linear_solver'] = 'mumps' p.driver.opt_settings['mu_strategy'] = 'adaptive' # p.driver.opt_settings['derivative_test'] = 'first-order' @@ -216,7 +216,7 @@ def test_mpl_plots(self): 'states_accel.png', 'controls_u1.png', 'states_vr.png', 'pos_x.png', 'states_vt.png', 'pos_y.png', 'parameters_u1.png', 'states_theta.png', 'control_rates_u1_rate2.png', 'state_rates_vt.png', 'time_phase.png', - 'parameters_c.png', 'state_rates_theta.png', 'state_rates_vr.png'} + 'parameters_c.png', 'state_rates_theta.png', 'state_rates_vr.png', 'dt_dstau.png'} self.assertSetEqual(expected_files, set(os.listdir('plots'))) diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py index e6143e454..51a99fd27 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py @@ -17,7 +17,7 @@ class TestExampleTwoBurnOrbitRaiseMPI(unittest.TestCase): N_PROCS = 3 - def test_ex_two_burn_orbit_raise(self): + def test_ex_two_burn_orbit_raise_mpi(self): optimizer = 'IPOPT' p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3, @@ -28,6 +28,17 @@ def test_ex_two_burn_orbit_raise(self): assert_near_equal(p.get_val('traj.burn2.states:deltav')[-1], 0.3995, tolerance=2.0E-3) + def test_ex_two_burn_orbit_raise_connected_mpi(self): + optimizer = 'IPOPT' + + p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3, + compressed=False, optimizer=optimizer, simulate=False, + connected=True, show_output=False) + + if p.model.traj.phases.burn2 in p.model.traj.phases._subsystems_myproc: + assert_near_equal(p.get_val('traj.burn2.states:deltav')[0], 0.3995, + tolerance=2.0E-3) + if __name__ == '__main__': # pragma: no cover unittest.main() diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py index 356a5239f..29815fb6c 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py @@ -38,8 +38,8 @@ def test_ex_two_burn_orbit_raise_connected(self): sim_case2 = om.CaseReader('dymos_simulation.db').get_case('final') # Verify that the second case has the same inputs and outputs - assert_cases_equal(case1, p, tol=1.0E-9) - assert_cases_equal(sim_case1, sim_case2, tol=1.0E-9) + assert_cases_equal(case1, p, tol=1.0E-8) + assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8) def test_restart_from_solution_radau(self): optimizer = 'IPOPT' @@ -94,8 +94,8 @@ def test_ex_two_burn_orbit_raise_connected(self): sim_case2 = om.CaseReader('dymos_simulation.db').get_case('final') # Verify that the second case has the same inputs and outputs - assert_cases_equal(case1, p, tol=1.0E-9) - assert_cases_equal(sim_case1, sim_case2, tol=1.0E-9) + assert_cases_equal(case1, p, tol=1.0E-8) + assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8) def test_restart_from_solution_radau_to_connected(self): optimizer = 'IPOPT' diff --git a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py index 9a69a655c..6c4298bdd 100644 --- a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py +++ b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py @@ -42,7 +42,7 @@ def make_problem(self, transcription=GaussLobatto, optimizer='SLSQP', numseg=30) elif optimizer == 'IPOPT': p.driver.opt_settings['hessian_approximation'] = 'limited-memory' # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling' - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['linear_solver'] = 'mumps' p.driver.declare_coloring() else: diff --git a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py index 94dfe8161..52c8b2e70 100644 --- a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py +++ b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py @@ -1,5 +1,6 @@ import unittest - +import numpy as np +from numpy.polynomial import Polynomial as P import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -7,9 +8,8 @@ from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse -@require_pyoptsparse(optimizer='SLSQP') def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', - transcription_order=3, force_alloc_complex=False): + transcription_order=3, force_alloc_complex=False, add_rate=False): p = om.Problem(model=om.Group()) @@ -25,6 +25,12 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', p.driver.opt_settings['Function precision'] = 1.0E-12 p.driver.opt_settings['Linesearch tolerance'] = 0.1 p.driver.opt_settings['Major step limit'] = 0.5 + elif optimizer == 'IPOPT': + p.driver.opt_settings['tol'] = 1.0E-5 + p.driver.opt_settings['print_level'] = 0 + p.driver.opt_settings['mu_strategy'] = 'monotone' + p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' + p.driver.opt_settings['mu_init'] = 0.01 t = {'gauss-lobatto': dm.GaussLobatto(num_segments=num_seg, order=transcription_order), 'radau-ps': dm.Radau(num_segments=num_seg, order=transcription_order)} @@ -86,6 +92,10 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', phase.add_timeseries_output(['aero.*', 'prop.thrust', 'prop.m_dot'], units={'aero.f_lift': 'lbf', 'prop.thrust': 'lbf'}) + # test adding rate as timeseries output + if add_rate: + phase.add_timeseries_rate_output('aero.mach') + p.model.linear_solver = om.DirectSolver() p.setup(check=True, force_alloc_complex=force_alloc_complex) @@ -108,48 +118,111 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', @use_tempdirs class TestMinTimeClimb(unittest.TestCase): - def test_results_gauss_lobatto(self): - p = min_time_climb(optimizer='SLSQP', num_seg=12, transcription_order=3, - transcription='gauss-lobatto') + def _test_results(self, p): + """ Verify the results of the optimization. """ + # Verify that ODE output mach is added to the timeseries + assert_near_equal(p.get_val('traj.phase0.timeseries.mach')[-1], 1.0, tolerance=1.0E-2) # Check that time matches to within 1% of an externally verified solution. assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 321.0, tolerance=0.02) - # Verify that ODE output mach is added to the timeseries - assert_near_equal(p.get_val('traj.phase0.timeseries.mach')[-1], 1.0, tolerance=1.0E-2) - - # verify all wildcard timeseries exist - output_dict = dict(p.model.list_outputs(units=True)) - ts = [k for k, v in output_dict.items() if 'timeseries' in k] + def _test_wilcard_outputs(self, p): + """ Test that all wilcard outputs are provided. """ + output_dict = dict(p.model.list_outputs(units=True, out_stream=None)) + ts = {k: v for k, v in output_dict.items() if 'timeseries.' in k} for c in ['mach', 'CD0', 'kappa', 'CLa', 'CL', 'CD', 'q', 'f_lift', 'f_drag', 'thrust', 'm_dot']: assert(any([True for t in ts if 'timeseries.' + c in t])) - # verify time series units + + def _test_timeseries_units(self, p): + """ Test that the units from the timeseries are correct. """ + output_dict = dict(p.model.list_outputs(units=True, out_stream=None)) assert(output_dict['traj.phases.phase0.timeseries.thrust']['units'] == 'lbf') # no wildcard, from units dict assert(output_dict['traj.phases.phase0.timeseries.m_dot']['units'] == 'kg/s') # no wildcard, from ODE assert(output_dict['traj.phases.phase0.timeseries.f_drag']['units'] == 'N') # wildcard, from ODE assert(output_dict['traj.phases.phase0.timeseries.f_lift']['units'] == 'lbf') # wildcard, from units dict + def _test_mach_rate(self, p, plot=False): + """ Test that the mach rate is provided by the timeseries and is accurate. """ + # Verify correct timeseries output of mach_rate + output_dict = dict(p.model.list_outputs(units=True, out_stream=None)) + ts = {k: v for k, v in output_dict.items() if 'timeseries.' in k} + self.assertTrue('traj.phases.phase0.timeseries.mach_rate' in ts) + + time = p['traj.phases.phase0.timeseries.time'][:, 0] + mach = p['traj.phases.phase0.timeseries.mach'][:, 0] + mach_rate = p['traj.phases.phase0.timeseries.mach_rate'][:, 0] + + # Fit a numpy polynomial segment by segment to mach vs time, and compare the derivatives to mach_rate + gd = p.model.traj.phases.phase0.options['transcription'].grid_data + seg_idxs = gd.subset_segment_indices['all'] + num_seg = seg_idxs.shape[0] + + if plot: + import matplotlib.pyplot as plt + import matplotlib.cm as cm + fig, axes = plt.subplots(2, 1, sharex=True) + color = iter(cm.viridis(np.linspace(0, 1, num_seg))) + axes[0].set_ylabel('Mach') + axes[1].set_ylabel('Mach rate (1/s)') + axes[1].set_xlabel('time (s)') + + for i in range(num_seg): + time_seg = time[seg_idxs[i, 0]: seg_idxs[i, 1]] + mach_seg = mach[seg_idxs[i, 0]: seg_idxs[i, 1]] + order = len(time_seg) - 1 + p = P.fit(time_seg, mach_seg, order) + deriv = p.deriv(1) + + time_nodes = time[seg_idxs[i, 0]: seg_idxs[i, 1]] + mach_nodes = mach[seg_idxs[i, 0]: seg_idxs[i, 1]] + mach_rate_nodes = mach_rate[seg_idxs[i, 0]: seg_idxs[i, 1]] + + if plot: + c = next(color) + t_plot = np.linspace(time_seg[0], time_seg[-1], 20) + m_plot = p(t_plot) + m_rate_plot = deriv(t_plot) + axes[0].plot(t_plot, m_plot, c=c) + axes[1].plot(t_plot, m_rate_plot, '--', c=c) + + axes[0].plot(time[seg_idxs[i, 0]: seg_idxs[i, 1]], mach[seg_idxs[i, 0]: seg_idxs[i, 1]], 'o', c=c) + axes[1].plot(time[seg_idxs[i, 0]: seg_idxs[i, 1]], mach_rate[seg_idxs[i, 0]: seg_idxs[i, 1]], 'o', c=c) + + assert_near_equal(mach_nodes, p(time_nodes), tolerance=1.0E-9) + assert_near_equal(mach_rate_nodes, deriv(time_nodes), tolerance=1.0E-9) + + if plot: + plt.show() + + @require_pyoptsparse(optimizer='SLSQP') + def test_results_gauss_lobatto(self): + NUM_SEG = 12 + ORDER = 3 + p = min_time_climb(optimizer='IPOPT', num_seg=NUM_SEG, transcription_order=ORDER, + transcription='gauss-lobatto', add_rate=True) + + self._test_results(p) + + self._test_wilcard_outputs(p) + + self._test_timeseries_units(p) + + self._test_mach_rate(p) + + @require_pyoptsparse(optimizer='SLSQP') def test_results_radau(self): - p = min_time_climb(optimizer='SLSQP', num_seg=12, transcription_order=3, - transcription='radau-ps') + NUM_SEG = 12 + ORDER = 3 + p = min_time_climb(optimizer='SLSQP', num_seg=NUM_SEG, transcription_order=ORDER, + transcription='radau-ps', add_rate=True) - # Check that time matches to within 1% of an externally verified solution. - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 321.0, tolerance=0.02) + self._test_results(p) - # Verify that ODE output mach is added to the timeseries - assert_near_equal(p.get_val('traj.phase0.timeseries.mach')[-1], 1.0, tolerance=1.0E-2) + self._test_wilcard_outputs(p) - # verify all wildcard timeseries exist - output_dict = dict(p.model.list_outputs(units=True)) - ts = [k for k, v in output_dict.items() if 'timeseries' in k] - for c in ['mach', 'CD0', 'kappa', 'CLa', 'CL', 'CD', 'q', 'f_lift', 'f_drag', 'thrust', 'm_dot']: - assert(any([True for t in ts if 'timeseries.' + c in t])) - # verify time series units - assert(output_dict['traj.phases.phase0.timeseries.thrust']['units'] == 'lbf') # no wildcard, from units dict - assert(output_dict['traj.phases.phase0.timeseries.m_dot']['units'] == 'kg/s') # no wildcard, from ODE - assert(output_dict['traj.phases.phase0.timeseries.f_drag']['units'] == 'N') # wildcard, from ODE - assert(output_dict['traj.phases.phase0.timeseries.f_lift']['units'] == 'lbf') # wildcard, from units dict + self._test_timeseries_units(p) + self._test_mach_rate(p) if __name__ == '__main__': # pragma: no cover unittest.main() diff --git a/dymos/examples/racecar/doc/test_doc_racecar.py b/dymos/examples/racecar/doc/test_doc_racecar.py index 32ad1af3f..a7261c24b 100644 --- a/dymos/examples/racecar/doc/test_doc_racecar.py +++ b/dymos/examples/racecar/doc/test_doc_racecar.py @@ -137,7 +137,7 @@ def test_racecar_for_docs(self): p.driver.opt_settings['compl_inf_tol'] = 1e-3 p.driver.opt_settings['acceptable_iter'] = 0 p.driver.opt_settings['tol'] = 1e-3 - p.driver.opt_settings['print_level'] = 5 + 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' diff --git a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py index 2a79813a3..992f34740 100644 --- a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py +++ b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py @@ -34,7 +34,7 @@ def make_problem(self, transcription=Radau, optimizer='SLSQP', numseg=30): p.driver.opt_settings['Verify level'] = 3 elif optimizer == 'IPOPT': p.driver.opt_settings['max_iter'] = 500 - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['tol'] = 1.0E-6 p.driver.opt_settings['acceptable_tol'] = 1.0E-5 traj = p.model.add_subsystem('traj', Trajectory()) diff --git a/dymos/examples/robot_arm/test/test_robot_arm.py b/dymos/examples/robot_arm/test/test_robot_arm.py index 84438e6dd..69a561189 100644 --- a/dymos/examples/robot_arm/test/test_robot_arm.py +++ b/dymos/examples/robot_arm/test/test_robot_arm.py @@ -38,7 +38,7 @@ def make_problem(self, transcription=Radau, optimizer='SLSQP', numseg=30): elif optimizer == 'IPOPT': p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' p.driver.opt_settings['max_iter'] = 500 - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 traj = p.model.add_subsystem('traj', Trajectory()) phase = traj.add_phase('phase', Phase(ode_class=RobotArmODE, diff --git a/dymos/examples/shuttle_reentry/test/test_reentry.py b/dymos/examples/shuttle_reentry/test/test_reentry.py index 2b0a3b7bc..76d9a0e6f 100644 --- a/dymos/examples/shuttle_reentry/test/test_reentry.py +++ b/dymos/examples/shuttle_reentry/test/test_reentry.py @@ -26,7 +26,7 @@ def make_problem(self, constrained=True, transcription=GaussLobatto, optimizer=' if optimizer == 'IPOPT': p.driver.opt_settings['max_iter'] = 500 p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas' - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' p.driver.opt_settings['tol'] = 1.0E-7 p.driver.opt_settings['mu_strategy'] = 'monotone' @@ -152,7 +152,7 @@ def test_reentry_mixed_controls(self): p.driver.declare_coloring(tol=1.0E-12) p.driver.options['optimizer'] = 'IPOPT' p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas' - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' p.driver.opt_settings['mu_strategy'] = 'monotone' diff --git a/dymos/examples/vanderpol/vanderpol_dymos.py b/dymos/examples/vanderpol/vanderpol_dymos.py index c9acfa75a..34a5e2572 100644 --- a/dymos/examples/vanderpol/vanderpol_dymos.py +++ b/dymos/examples/vanderpol/vanderpol_dymos.py @@ -20,7 +20,7 @@ def vanderpol(transcription='gauss-lobatto', num_segments=40, transcription_orde if optimizer == 'SNOPT': p.driver.opt_settings['iSumm'] = 6 # show detailed SNOPT output elif optimizer == 'IPOPT': - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.declare_coloring() # define a Trajectory object and add to model diff --git a/dymos/examples/water_rocket/doc/test_doc_water_rocket.py b/dymos/examples/water_rocket/doc/test_doc_water_rocket.py index db405f857..41fefd39e 100644 --- a/dymos/examples/water_rocket/doc/test_doc_water_rocket.py +++ b/dymos/examples/water_rocket/doc/test_doc_water_rocket.py @@ -25,7 +25,7 @@ def test_water_rocket_height_for_docs(self): traj = p.model.add_subsystem('traj', traj) p.driver = om.pyOptSparseDriver(optimizer='IPOPT', print_results=False) - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['max_iter'] = 1000 p.driver.opt_settings['mu_strategy'] = 'monotone' p.driver.declare_coloring(tol=1.0E-12) @@ -64,7 +64,7 @@ def test_water_rocket_range_for_docs(self): traj = p.model.add_subsystem('traj', traj) p.driver = om.pyOptSparseDriver(optimizer='IPOPT') - p.driver.opt_settings['print_level'] = 5 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['max_iter'] = 1000 p.driver.opt_settings['mu_strategy'] = 'monotone' p.driver.declare_coloring(tol=1.0E-12) diff --git a/dymos/phase/options.py b/dymos/phase/options.py index 7efab2bc3..a7cbd373a 100644 --- a/dymos/phase/options.py +++ b/dymos/phase/options.py @@ -685,3 +685,6 @@ def __init__(self, read_only=False): self.declare(name='units', default=None, allow_none=True, desc='Units to be used for the timeseries output, or None to leave the units unchanged.') + + self.declare(name='is_rate', default=False, allow_none=False, + desc='If True this is a rate.') diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index 1f622b851..86f1eb5a9 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -75,7 +75,7 @@ def __init__(self, from_phase=None, **kwargs): self._path_constraints = [] self._timeseries = {'timeseries': {'transcription': None, 'subset': 'all', - 'outputs': []}} + 'outputs': {}}} self._objectives = {} else: self.time_options.update(from_phase.time_options) @@ -613,7 +613,7 @@ def add_polynomial_control(self, name, order, desc=_unspecified, val=_unspecifie name : str Name of the controllable parameter in the ODE. order : int - The order of the interpolating polynomial used to represent the control valeu in + The order of the interpolating polynomial used to represent the control value in phase tau space. desc : str A description of the polynomial control. @@ -1050,8 +1050,7 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None, # Automatically add the requested variable to the timeseries outputs if it's an ODE output. var_type = self.classify_var(name) if var_type == 'ode': - current_output_names = [op['output_name'] for op in self._timeseries['timeseries']['outputs']] - if constraint_name not in current_output_names: + if constraint_name not in self._timeseries['timeseries']['outputs']: self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape) def add_path_constraint(self, name, constraint_name=None, units=None, shape=None, indices=None, @@ -1135,8 +1134,7 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None # Automatically add the requested variable to the timeseries outputs if it's an ODE output. var_type = self.classify_var(name) if var_type == 'ode': - current_output_names = [op['output_name'] for op in self._timeseries['timeseries']['outputs']] - if constraint_name not in current_output_names: + if constraint_name not in self._timeseries['timeseries']['outputs']: self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape) def add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, @@ -1149,7 +1147,61 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap name : str, or list of str The name(s) of the variable to be used as a timeseries output. Must be one of 'time', 'time_phase', one of the states, controls, control rates, or parameters, - in the phase, or the path to an output variable in the ODE. + in the phase, the path to an output variable in the ODE, or a glob pattern + matching some outputs in the ODE. + output_name : str or None or list or dict + The name of the variable as listed in the phase timeseries outputs. By + default this is the last element in `name` when split by dots. The user may + override the constraint name if splitting the path causes name collisions. + units : str or None or _unspecified + The units to express the timeseries output. If None, use the + units associated with the target. If provided, must be compatible with + the target units. + If a list of names is provided, units can be a matching list or dictionary. + shape : tuple or _unspecified + The shape of the timeseries output variable. This must be provided (if not scalar) + since Dymos doesn't necessarily know the shape of ODE outputs until setup time. + timeseries : str or None + The name of the timeseries to which the output is being added. + """ + if type(name) is list: + for i, name_i in enumerate(name): + if type(units) is dict: # accept dict for units when using array of name + unit = units.get(name_i, None) + elif type(units) is list: # allow matching list for units + unit = units[i] + else: + unit = units + + oname = self._add_timeseries_output(name_i, output_name=output_name, + units=unit, + shape=shape, + timeseries=timeseries, + rate=False) + + # Handle specific units for wildcard names. + if oname is not None and '*' in name_i: + self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units + + else: + self._add_timeseries_output(name, output_name=output_name, + units=units, + shape=shape, + timeseries=timeseries, + rate=False) + + def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, + timeseries='timeseries'): + r""" + Add the rate of a variable to the timeseries outputs of the phase. + + Parameters + ---------- + name : str, or list of str + The name(s) of the variable to be used as a timeseries output. Must be one of + 'time', 'time_phase', one of the states, controls, control rates, or parameters, + in the phase, the path to an output variable in the ODE, or a glob pattern + matching some outputs in the ODE. output_name : str or None or list or dict The name of the variable as listed in the phase timeseries outputs. By default this is the last element in `name` when split by dots. The user may @@ -1174,27 +1226,30 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap else: unit = units - self._add_timeseries_output(name_i, output_name=output_name, - units=unit, - shape=shape, - timeseries=timeseries) + oname = self._add_timeseries_output(name_i, output_name=output_name, + units=unit, + shape=shape, + timeseries=timeseries, + rate=True) # Handle specific units for wildcard names. - if '*' in name_i: - self._timeseries[timeseries]['outputs'][-1]['wildcard_units'] = units + if oname is not None and '*' in name_i: + self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units else: self._add_timeseries_output(name, output_name=output_name, units=units, shape=shape, - timeseries=timeseries) + timeseries=timeseries, + rate=True) def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, - timeseries='timeseries'): + timeseries='timeseries', rate=False): r""" - Add a single variable to the timeseries outputs of the phase. + Add a single variable or rate to the timeseries outputs of the phase. - This is called by add_timeseries_output for each variable that is added. + This is called by add_timeseries_output or add_timeseries_rate_output for each variable or rate + that is added. Parameters ---------- @@ -1205,7 +1260,8 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha output_name : str or None The name of the variable as listed in the phase timeseries outputs. By default this is the last element in `name` when split by dots. The user may - override the constraint name if splitting the path causes name collisions. + override the constraint name if splitting the path causes name collisions. If rate + is True, the rate name will be this name + _rate. units : str or None The units to express the timeseries output. If None, use the units associated with the target. If provided, must be compatible with @@ -1215,27 +1271,39 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha since Dymos doesn't necessarily know the shape of ODE outputs until setup time. timeseries : str or None The name of the timeseries to which the output is being added. - """ - if output_name is None: - output_name = name.split('.')[-1] + rate : bool + If True, add the rate of change of the named variable to the timeseries outputs of the + phase. The rate variable will be named f'{name}_rate'. Defaults to False. + Returns + ------- + str or None + Name of output that was added to the timeseries or None if nothing was added. + """ if timeseries not in self._timeseries: raise ValueError(f'Timeseries {timeseries} does not exist in phase {self.pathname}') - current_output_names = [op['output_name'] for op in self._timeseries[timeseries]['outputs']] + if output_name is None: + output_name = name.split('.')[-1] + + if rate: + output_name = output_name + '_rate' - if output_name not in current_output_names: + if output_name in self._timeseries[timeseries]['outputs']: + om.issue_warning(f'Output name `{output_name}` is already in timeseries `{timeseries}`. ' + f'New output ignored.') + else: ts_output = TimeseriesOutputOptionsDictionary() ts_output['name'] = name ts_output['output_name'] = output_name ts_output['wildcard_units'] = {} ts_output['units'] = units ts_output['shape'] = shape + ts_output['is_rate'] = rate - self._timeseries[timeseries]['outputs'].append(ts_output) - else: - om.issue_warning(f'Output name `{output_name}` is already in timeseries `{timeseries}`. ' - f'New output ignored.') + self._timeseries[timeseries]['outputs'][output_name] = ts_output + + return output_name def add_timeseries(self, name, transcription, subset='all'): r""" @@ -1254,7 +1322,7 @@ def add_timeseries(self, name, transcription, subset='all'): """ self._timeseries[name] = {'transcription': transcription, 'subset': subset, - 'outputs': []} + 'outputs': {}} def add_objective(self, name, loc='final', index=None, shape=(1,), ref=None, ref0=None, adder=None, scaler=None, parallel_deriv_color=None): @@ -1585,13 +1653,11 @@ def configure_state_discovery(self): """ transcription = self.options['transcription'] state_options = self.state_options - out_meta = get_promoted_vars(transcription._get_ode(self), 'output') + out_meta = get_promoted_vars(transcription._get_ode(self), 'output', metadata_keys=('tags',)) - for name, meta in out_meta.items(): - tags = meta['tags'] - prom_name = meta['prom_name'] + for prom_name, meta in out_meta.items(): state = None - for tag in sorted(tags): + for tag in sorted(meta['tags']): # Declared as rate_source. if tag.startswith('dymos.state_rate_source:') or tag.startswith('state_rate_source:'): @@ -1607,7 +1673,7 @@ def configure_state_discovery(self): if state_options[state]['rate_source'] is not None: if state_options[state]['rate_source'] != prom_name: raise ValueError(f"rate_source has been declared twice for state " - f"'{state}' which is tagged on '{name}'.") + f"'{state}' which is tagged on '{prom_name}'.") state_options[state]['rate_source'] = prom_name @@ -2275,7 +2341,6 @@ def _indices_in_constraints(self, name, loc): all_flat_idxs : set A C-order flattened set of indices that apply to the constraint. """ - s = {'initial': 'initial boundary', 'final': 'final boundary', 'path': 'path'} cons = {'initial': self._initial_boundary_constraints, 'final': self._final_boundary_constraints, 'path': self._path_constraints} @@ -2285,12 +2350,15 @@ def _indices_in_constraints(self, name, loc): for con in cons[loc]: if con['name'] != name: continue + flat_idxs = get_constraint_flat_idxs(con) duplicate_idxs = all_flat_idxs.intersection(flat_idxs) if duplicate_idxs: + s = {'initial': 'initial boundary', 'final': 'final boundary', 'path': 'path'} raise ValueError(f'Duplicate constraint in phase {self.pathname}. ' f'The following indices of `{name}` are used in ' f'multiple {s[loc]} constraints:\n{duplicate_idxs}') + all_flat_idxs.update(flat_idxs) return all_flat_idxs diff --git a/dymos/test/test_run_problem.py b/dymos/test/test_run_problem.py index 779747211..5f4067bc8 100755 --- a/dymos/test/test_run_problem.py +++ b/dymos/test/test_run_problem.py @@ -35,7 +35,7 @@ def test_run_HS_problem_radau(self): elif optimizer == 'IPOPT': p.driver.opt_settings['hessian_approximation'] = 'limited-memory' # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling' - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['max_iter'] = 200 p.driver.opt_settings['linear_solver'] = 'mumps' @@ -102,7 +102,7 @@ def test_run_HS_problem_gl(self): elif optimizer == 'IPOPT': p.driver.opt_settings['hessian_approximation'] = 'limited-memory' # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling' - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['linear_solver'] = 'mumps' traj = p.model.add_subsystem('traj', dm.Trajectory()) diff --git a/dymos/trajectory/test/test_trajectory_parameters.py b/dymos/trajectory/test/test_trajectory_parameters.py index 89446dc9b..76487855d 100644 --- a/dymos/trajectory/test/test_trajectory_parameters.py +++ b/dymos/trajectory/test/test_trajectory_parameters.py @@ -175,7 +175,7 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP elif optimizer == 'IPOPT': p.driver.opt_settings['hessian_approximation'] = 'limited-memory' p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' - p.driver.opt_settings['print_level'] = 4 + p.driver.opt_settings['print_level'] = 0 p.driver.opt_settings['linear_solver'] = 'mumps' p.driver.opt_settings['mu_strategy'] = 'adaptive' # p.driver.opt_settings['derivative_test'] = 'first-order' diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py index c1bbb1a20..de661501c 100644 --- a/dymos/trajectory/trajectory.py +++ b/dymos/trajectory/trajectory.py @@ -35,14 +35,27 @@ class Trajectory(om.Group): ---------- **kwargs : dict Dictionary of optional arguments. + + Attributes + ---------- + parameter_options : dict + A dictionary of parameter names and their associated TrajectoryParameterOptionsDictionary + phases : om.Group or om.ParallelGroup + The Group which contains phases for this Trajectory. + + _linkages : OrderedDict + A dictionary containing phase linkage information for the Trajectory. + _phases : dict + A dictionary of phase names as keys with the Phase objects being their associated values. """ def __init__(self, **kwargs): super(Trajectory, self).__init__(**kwargs) + self._linkages = {} + self._phases = {} + self.parameter_options = {} - self._linkages = OrderedDict() - self._phases = OrderedDict() - self._phase_add_kwargs = {} + self.phases = om.ParallelGroup() def initialize(self): """ @@ -71,8 +84,7 @@ def add_phase(self, name, phase, **kwargs): PhaseBase The Phase object added to the trajectory. """ - self._phases[name] = phase - self._phase_add_kwargs[name] = kwargs + self._phases[name] = self.phases.add_subsystem(name, phase, **kwargs) return phase def set_parameter_options(self, name, units=_unspecified, val=_unspecified, desc=_unspecified, opt=False, @@ -324,10 +336,8 @@ def setup(self): if self.parameter_options: self._setup_parameters() - phases_group = self.add_subsystem('phases', subsys=om.ParallelGroup()) - - for name, phs in self._phases.items(): - phases_group.add_subsystem(name, phs, **self._phase_add_kwargs[name]) + # This will override the existing phases attribute with the same thing. + self.add_subsystem('phases', subsys=self.phases) if self._linkages: self._setup_linkages() diff --git a/dymos/transcriptions/common/timeseries_output_comp.py b/dymos/transcriptions/common/timeseries_output_comp.py index 1afb28894..8a0fbdee6 100644 --- a/dymos/transcriptions/common/timeseries_output_comp.py +++ b/dymos/transcriptions/common/timeseries_output_comp.py @@ -28,8 +28,6 @@ def initialize(self): """ Declare component options. """ - self._timeseries_outputs = [] - self._vars = {} self.options.declare('input_grid_data', diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py index 804d0b854..eb5c462d7 100644 --- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py +++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py @@ -570,7 +570,7 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None): units = phase.state_options[var]['units'] linear = loc == 'initial' obj_path = f'integrator.states_out:{var}' - elif var_type in 'indep_control': + elif var_type == 'indep_control': shape = phase.control_options[var]['shape'] units = phase.control_options[var]['units'] linear = True @@ -580,7 +580,7 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None): units = phase.control_options[var]['units'] linear = False obj_path = f'control_values:{var}' - elif var_type in 'indep_polynomial_control': + elif var_type == 'indep_polynomial_control': shape = phase.polynomial_control_options[var]['shape'] units = phase.polynomial_control_options[var]['units'] linear = True diff --git a/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py b/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py index 666dd05c8..778fc0d6a 100644 --- a/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py +++ b/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py @@ -31,7 +31,6 @@ def __init__(self, **kwargs): # from another variable has potentially different units self._units = {} self._conversion_factors = {} - self._vars = {} self._no_check_partials = not dymos_options['include_check_partials'] @@ -50,7 +49,7 @@ def setup(self): self.input_num_nodes = igd.subset_num_nodes['segment_ends'] self.output_num_nodes = self.input_num_nodes - def _add_output_configure(self, name, units, shape, desc='', src=None): + def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False): """ Add a single timeseries output. @@ -70,12 +69,18 @@ def _add_output_configure(self, name, units, shape, desc='', src=None): description of the timeseries output variable. src : str The src path of the variables input, used to prevent redundant inputs. + rate : bool + If True, timeseries output is a rate. Returns ------- bool True if a new input was added for the output, or False if it reuses an existing input. """ + if rate: + raise NotImplementedError("Timeseries output rates are not currently supported for " + "ExplicitShooting transcriptions.") + input_num_nodes = self.input_num_nodes output_num_nodes = self.output_num_nodes added_source = False @@ -102,7 +107,7 @@ def _add_output_configure(self, name, units, shape, desc='', src=None): shape=(output_num_nodes,) + shape, units=units, desc=desc) - self._vars[name] = (input_name, output_name, shape) + self._vars[name] = (input_name, output_name, shape, rate) size = np.prod(shape) rs = cs = np.arange(output_num_nodes * size, dtype=int) @@ -114,7 +119,7 @@ def _add_output_configure(self, name, units, shape, desc='', src=None): offset = 0 else: scale, offset = unit_conversion(input_units, units) - self._conversion_factors[output_name] = scale, offset + self._conversion_factors[output_name] = scale, offset self.declare_partials(of=output_name, wrt=input_name, @@ -133,6 +138,9 @@ def compute(self, inputs, outputs): outputs : `Vector` `Vector` containing outputs. """ - for (input_name, output_name, _) in self._vars.values(): - scale, offset = self._conversion_factors[output_name] - outputs[output_name] = scale * (inputs[input_name] + offset) + for (input_name, output_name, _, _) in self._vars.values(): + if output_name in self._conversion_factors: + scale, offset = self._conversion_factors[output_name] + outputs[output_name] = scale * (inputs[input_name] + offset) + else: + outputs[output_name] = inputs[input_name] diff --git a/dymos/transcriptions/explicit_shooting/rk_integration_comp.py b/dymos/transcriptions/explicit_shooting/rk_integration_comp.py index 099c5e6eb..9ef7a749c 100644 --- a/dymos/transcriptions/explicit_shooting/rk_integration_comp.py +++ b/dymos/transcriptions/explicit_shooting/rk_integration_comp.py @@ -577,22 +577,22 @@ def _configure_timeseries_outputs(self): self._timeseries_idxs_in_y = {} self._filtered_timeseries_outputs = {} - for ts_name, ts_opts in self.timeseries_options.items(): - patterns = [output['name'] for output in ts_opts['outputs']] + for ts_opts in self.timeseries_options.values(): + patterns = [output['name'] for output in ts_opts['outputs'].values()] matching_outputs = filter_outputs(patterns, ode_outputs) explicit_requests = set([var for var in patterns if '*' not in var]) - unmatched_requests = sorted(list(set(explicit_requests) - set(matching_outputs.keys()))) + unmatched_requests = explicit_requests.difference(matching_outputs) if unmatched_requests: om.issue_warning(msg='The following timeseries outputs were requested but ' - f'not found in the ODE: {", ".join(unmatched_requests)}', + f'not found in the ODE: {", ".join(sorted(unmatched_requests))}', category=om.OpenMDAOWarning) for var, var_meta in matching_outputs.items(): if var in explicit_requests: - ts_output = next((output for output in ts_opts['outputs'] if output['name'] == var)) + ts_output = next((output for output in ts_opts['outputs'].values() if output['name'] == var)) # var explicitly matched output_name = ts_output['output_name'] if ts_output['output_name'] else ts_output['name'].split('.')[-1] units = ts_output['units'] or var_meta.get('units', None) @@ -1545,12 +1545,12 @@ def compute(self, inputs, outputs): outputs['time_phase'] = self._t[idxs, ...] - inputs['t_initial'] # Extract the state values - for state_name, options in self.state_options.items(): + for state_name in self.state_options: of = self._state_output_names[state_name] outputs[of] = self._x[idxs, self.state_idxs[state_name]] # Extract the control values and rates - for control_name, options in self.control_options.items(): + for control_name in self.control_options: oname = self._control_output_names[control_name] rate_name = self._control_rate_names[control_name] rate2_name = self._control_rate2_names[control_name] @@ -1559,7 +1559,7 @@ def compute(self, inputs, outputs): outputs[rate2_name] = self._y[idxs, self._control_rate2_idxs_in_y[control_name]] # Extract the control values and rates - for control_name, options in self.polynomial_control_options.items(): + for control_name in self.polynomial_control_options: oname = self._polynomial_control_output_names[control_name] rate_name = self._polynomial_control_rate_names[control_name] rate2_name = self._polynomial_control_rate2_names[control_name] @@ -1568,7 +1568,7 @@ def compute(self, inputs, outputs): outputs[rate2_name] = self._y[idxs, self._polynomial_control_rate2_idxs_in_y[control_name]] # Extract the timeseries outputs - for name, options in self._filtered_timeseries_outputs.items(): + for name in self._filtered_timeseries_outputs: oname = self._timeseries_output_names[name] outputs[oname] = self._y[idxs, self._timeseries_idxs_in_y[name]] diff --git a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py index 77b40702c..782efcde5 100644 --- a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py +++ b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py @@ -40,6 +40,14 @@ def __init__(self, **kwargs): # Flag to set if no multiplication by the interpolation matrix is necessary self._no_interp = False + def _declare_options(self): + """ + Declare options before kwargs are processed in the init method. + """ + super()._declare_options() + self.options.declare('time_units', default=None, allow_none=True, types=str, + desc='Units of time') + def setup(self): """ Define the independent variables as output variables. @@ -62,7 +70,8 @@ def setup(self): # To do this, find the nodes in the output grid which fall in each segment of the input # grid. Then build a Lagrange interpolating polynomial for that segment L_blocks = [] - output_nodes_ptau = ogd.node_ptau[ogd.subset_node_indices[output_subset]].tolist() + D_blocks = [] + output_nodes_ptau = ogd.node_ptau[ogd.subset_node_indices[output_subset]] for iseg in range(igd.num_segments): i1, i2 = igd.segment_indices[iseg] @@ -75,32 +84,31 @@ def setup(self): else: ptau_hi = igd.segment_ends[iseg+1] if iseg < igd.num_segments - 1: - idxs_in_iseg = np.where(output_nodes_ptau <= ptau_hi)[0] + optau_segi = output_nodes_ptau[output_nodes_ptau <= ptau_hi] else: - idxs_in_iseg = np.arange(len(output_nodes_ptau)) - optau_segi = np.asarray(output_nodes_ptau)[idxs_in_iseg] + optau_segi = output_nodes_ptau + # Remove the captured nodes so we don't accidentally include them again - output_nodes_ptau = output_nodes_ptau[len(idxs_in_iseg):] + output_nodes_ptau = output_nodes_ptau[len(optau_segi):] # # Now get the output nodes which fall in iseg in iseg's segment tau space. ostau_segi = 2.0 * (optau_segi - iptau_segi[0]) / (iptau_segi[-1] - iptau_segi[0]) - 1 # Create the interpolation matrix and add it to the blocks - L, _ = lagrange_matrices(istau_segi, ostau_segi) + L, D = lagrange_matrices(istau_segi, ostau_segi) L_blocks.append(L) + D_blocks.append(D) if _USE_SPARSE: self.interpolation_matrix = sp.block_diag(L_blocks, format='csr') + self.differentiation_matrix = sp.block_diag(D_blocks, format='csr') else: self.interpolation_matrix = block_diag(*L_blocks) + self.differentiation_matrix = block_diag(*D_blocks) - for (name, kwargs) in self._timeseries_outputs: - units = kwargs['units'] - desc = kwargs['units'] - shape = kwargs['shape'] - self._add_output_configure(name, units, shape, desc) + self.add_input('dt_dstau', shape=(self.input_num_nodes,), units=self.options['time_units']) - def _add_output_configure(self, name, units, shape, desc='', src=None): + def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False): """ Add a single timeseries output. @@ -120,6 +128,8 @@ def _add_output_configure(self, name, units, shape, desc='', src=None): description of the timeseries output variable. src : str The src path of the variables input, used to prevent redundant inputs. + rate : bool + If True, timeseries output is a rate. Returns ------- @@ -148,40 +158,40 @@ def _add_output_configure(self, name, units, shape, desc='', src=None): added_source = True output_name = name - self.add_output(output_name, - shape=(output_num_nodes,) + shape, - units=units, desc=desc) + self.add_output(output_name, shape=(output_num_nodes,) + shape, units=units, desc=desc) - self._vars[name] = (input_name, output_name, shape) + self._vars[name] = (input_name, output_name, shape, rate) size = np.prod(shape) - val_jac = np.zeros((output_num_nodes, size, input_num_nodes, size)) + jac = np.zeros((output_num_nodes, size, input_num_nodes, size)) + + if rate: + mat = self.differentiation_matrix + else: + mat = self.interpolation_matrix for i in range(size): if _USE_SPARSE: - val_jac[:, i, :, i] = self.interpolation_matrix.toarray() + jac[:, i, :, i] = mat.toarray() else: - val_jac[:, i, :, i] = self.interpolation_matrix - - val_jac = val_jac.reshape((output_num_nodes * size, input_num_nodes * size), - order='C') + jac[:, i, :, i] = mat - val_jac_rows, val_jac_cols = np.where(val_jac != 0) + jac = jac.reshape((output_num_nodes * size, input_num_nodes * size), order='C') - rs, cs = val_jac_rows, val_jac_cols + jac_rows, jac_cols = np.nonzero(jac) # There's a chance that the input for this output was pulled from another variable with # different units, so account for that with a conversion. - if None in {input_units, units}: - scale = 1.0 - offset = 0 + if input_units is None or units is None: + self.declare_partials(of=output_name, wrt=input_name, + rows=jac_rows, cols=jac_cols, val=jac[jac_rows, jac_cols]) else: scale, offset = unit_conversion(input_units, units) - self._conversion_factors[output_name] = scale, offset + self._conversion_factors[output_name] = scale, offset - self.declare_partials(of=output_name, - wrt=input_name, - rows=rs, cols=cs, val=scale * val_jac[rs, cs]) + self.declare_partials(of=output_name, wrt=input_name, + rows=jac_rows, cols=jac_cols, + val=scale * jac[jac_rows, jac_cols]) return added_source @@ -196,9 +206,10 @@ def compute(self, inputs, outputs): outputs : `Vector` `Vector` containing outputs. """ - for (input_name, output_name, _) in self._vars.values(): - scale, offset = self._conversion_factors[output_name] + # convert dt_dstau to a column vector + dt_dstau = inputs['dt_dstau'][:, np.newaxis] + for (input_name, output_name, _, is_rate) in self._vars.values(): if self._no_interp: interp_vals = inputs[input_name] else: @@ -208,4 +219,12 @@ def compute(self, inputs, outputs): inp = inp.swapaxes(0, 1) interp_vals = self.interpolation_matrix.dot(inp) - outputs[output_name] = scale * (interp_vals + offset) + + if is_rate: + interp_vals = self.differentiation_matrix.dot(interp_vals) / dt_dstau + + if output_name in self._conversion_factors: + scale, offset = self._conversion_factors[output_name] + outputs[output_name] = scale * (interp_vals + offset) + else: + outputs[output_name] = interp_vals diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py index 58aa9dfa2..bc250ab99 100644 --- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py +++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py @@ -3,7 +3,7 @@ import numpy as np import openmdao.api as om -from openmdao.utils.general_utils import simple_warning +from openmdao.utils.om_warnings import issue_warning from .pseudospectral_base import PseudospectralBase from .components import GaussLobattoInterleaveComp @@ -518,12 +518,12 @@ def configure_timeseries_outputs(self, phase): # phase.promotes(timeseries_name, inputs=[(tgt_name, prom_name)], # src_indices=(src_idxs,), flat_src_indices=True) - for ts_output in phase._timeseries[timeseries_name]['outputs']: + for output_name, ts_output in phase._timeseries[timeseries_name]['outputs'].items(): var = ts_output['name'] - output_name = ts_output['output_name'] units = ts_output['units'] wildcard_units = ts_output['wildcard_units'] shape = ts_output['shape'] + is_rate = ts_output['is_rate'] if '*' in var: # match outputs from the ODE matches = filter(list(ode_outputs.keys()), var) @@ -543,10 +543,10 @@ def configure_timeseries_outputs(self, phase): for output_name, var_list in output_name_groups.items(): if len(var_list) > 1: var_list_as_string = ', '.join(var_list) - simple_warning(f"The timeseries variable name {output_name} is " - f"duplicated in these variables: {var_list_as_string}. " - "Disambiguate by using the add_timeseries_output " - "output_name option.") + issue_warning(f"The timeseries variable name {output_name} is " + f"duplicated in these variables: {var_list_as_string}. " + "Disambiguate by using the add_timeseries_output " + "output_name option.") else: matches = [var] @@ -585,7 +585,8 @@ def configure_timeseries_outputs(self, phase): timeseries_input_added = timeseries_comp._add_output_configure(output_name, units, shape, src=f'interleave_comp.' - f'all_values:{output_name}') + f'all_values:{output_name}', + rate=is_rate) interleave_comp = phase._get_subsystem('interleave_comp') src_added = interleave_comp.add_var(output_name, shape, units, diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py index 32a48711f..fad215da0 100644 --- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py +++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py @@ -45,13 +45,13 @@ def setup_time(self, phase): phase : dymos.Phase The phase object to which this transcription instance applies. """ - time_units = phase.time_options['units'] grid_data = self.grid_data super(PseudospectralBase, self).setup_time(phase) time_comp = TimeComp(num_nodes=grid_data.num_nodes, node_ptau=grid_data.node_ptau, - node_dptau_dstau=grid_data.node_dptau_dstau, units=time_units, + node_dptau_dstau=grid_data.node_dptau_dstau, + units=phase.time_options['units'], initial_val=phase.time_options['initial_val'], duration_val=phase.time_options['duration_val']) @@ -252,14 +252,11 @@ def setup_ode(self, phase): The phase object to which this transcription instance applies. """ grid_data = self.grid_data - transcription = grid_data.transcription - time_units = phase.time_options['units'] - phase.add_subsystem('state_interp', subsys=StateInterpComp(grid_data=grid_data, state_options=phase.state_options, - time_units=time_units, - transcription=transcription)) + time_units=phase.time_options['units'], + transcription=grid_data.transcription)) def configure_ode(self, phase): """ @@ -292,14 +289,10 @@ def setup_defects(self, phase): phase : dymos.Phase The phase object to which this transcription instance applies. """ - grid_data = self.grid_data - - time_units = phase.time_options['units'] - phase.add_subsystem('collocation_constraint', - CollocationComp(grid_data=grid_data, + CollocationComp(grid_data=self.grid_data, state_options=phase.state_options, - time_units=time_units)) + time_units=phase.time_options['units'])) def _configure_solve_segments(self, state_name, options, phase): """ @@ -498,9 +491,12 @@ def setup_timeseries_outputs(self, phase): timeseries_comp = PseudospectralTimeseriesOutputComp(input_grid_data=gd, output_grid_data=ogd, - output_subset=options['subset']) + output_subset=options['subset'], + time_units=phase.time_options['units']) phase.add_subsystem(name, subsys=timeseries_comp) + phase.connect('dt_dstau', (f'{name}.dt_dstau'), flat_src_indices=True) + def _get_objective_src(self, var, loc, phase, ode_outputs=None): """ Return the path to the variable that will be used as the objective. @@ -557,7 +553,7 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None): else: linear = False constraint_path = f'states:{var}' - elif var_type in 'indep_control': + elif var_type == 'indep_control': shape = phase.control_options[var]['shape'] units = phase.control_options[var]['units'] linear = True @@ -567,7 +563,7 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None): units = phase.control_options[var]['units'] linear = False constraint_path = f'control_values:{var}' - elif var_type in 'indep_polynomial_control': + elif var_type == 'indep_polynomial_control': shape = phase.polynomial_control_options[var]['shape'] units = phase.polynomial_control_options[var]['units'] linear = True diff --git a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py index baca78a70..4f6012805 100644 --- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py +++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py @@ -4,7 +4,7 @@ import numpy as np import openmdao.api as om -from openmdao.utils.general_utils import simple_warning +from openmdao.utils.om_warnings import issue_warning from .pseudospectral_base import PseudospectralBase from ..common import RadauPSContinuityComp @@ -388,12 +388,12 @@ def configure_timeseries_outputs(self, phase): src_indices=(src_idxs,), flat_src_indices=True) - for ts_output in phase._timeseries[timeseries_name]['outputs']: + for output_name, ts_output in phase._timeseries[timeseries_name]['outputs'].items(): var = ts_output['name'] - output_name = ts_output['output_name'] units = ts_output['units'] wildcard_units = ts_output['wildcard_units'] shape = ts_output['shape'] + is_rate = ts_output['is_rate'] if '*' in var: # match outputs from the ODE matches = filter(list(ode_outputs.keys()), var) @@ -413,10 +413,10 @@ def configure_timeseries_outputs(self, phase): for output_name, var_list in output_name_groups.items(): if len(var_list) > 1: var_list_as_string = ', '.join(var_list) - simple_warning(f"The timeseries variable name {output_name} is " - f"duplicated in these variables: {var_list_as_string}. " - "Disambiguate by using the add_timeseries_output " - "output_name option.") + issue_warning(f"The timeseries variable name {output_name} is " + f"duplicated in these variables: {var_list_as_string}. " + "Disambiguate by using the add_timeseries_output " + "output_name option.") else: matches = [var] @@ -453,7 +453,8 @@ def configure_timeseries_outputs(self, phase): add_connection = timeseries_comp._add_output_configure(output_name, units, shape, desc='', - src=f'rhs_all.{v}') + src=f'rhs_all.{v}', + rate=is_rate) if add_connection: phase.connect(src_name=f'rhs_all.{v}', diff --git a/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py b/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py index 62afbc6af..3c5c99876 100644 --- a/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py +++ b/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py @@ -31,13 +31,7 @@ def setup(self): else: self.num_nodes = grid_data.num_segments * self.options['output_nodes_per_seg'] - for (name, kwargs) in self._timeseries_outputs: - units = kwargs['units'] - desc = kwargs['units'] - shape = kwargs['shape'] - self._add_output_configure(name, units, shape, desc) - - def _add_output_configure(self, name, units, shape, desc): + def _add_output_configure(self, name, units, shape, desc, rate=False): """ Add a single timeseries output. @@ -55,20 +49,20 @@ def _add_output_configure(self, name, units, shape, desc): Default is None, which means it has no units. desc : str description of the timeseries output variable. + rate : bool + If True, timeseries output is a rate. """ - num_nodes = self.num_nodes + if rate: + raise NotImplementedError("Timeseries output rates are not currently supported for " + "SolveIVP transcriptions.") + nodeshape = (self.num_nodes,) input_name = f'all_values:{name}' - self.add_input(input_name, - shape=(num_nodes,) + shape, - units=units, desc=desc) - output_name = name - self.add_output(output_name, - shape=(num_nodes,) + shape, - units=units, desc=desc) + self.add_input(input_name, shape=nodeshape + shape, units=units, desc=desc) + self.add_output(name, shape=nodeshape + shape, units=units, desc=desc) - self._vars[name] = (input_name, output_name, shape) + self._vars[name] = (input_name, name, shape, rate) def compute(self, inputs, outputs): """ @@ -81,5 +75,4 @@ def compute(self, inputs, outputs): outputs : `Vector` `Vector` containing outputs. """ - for (input_name, output_name, _) in self._vars.values(): - outputs[output_name] = inputs[input_name] + outputs.set_val(inputs.asarray()) diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py index a890cfc41..be10e8f6e 100644 --- a/dymos/transcriptions/solve_ivp/solve_ivp.py +++ b/dymos/transcriptions/solve_ivp/solve_ivp.py @@ -633,22 +633,18 @@ def configure_timeseries_outputs(self, phase): src_idxs_raw = np.zeros(num_seg * output_nodes_per_seg, dtype=int) src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']).ravel() - # tgt_name = f'all_values:parameters:{name}' - # phase.promotes('timeseries', inputs=[(tgt_name, prom_name)], - # src_indices=(src_idxs,), flat_src_indices=True) - # print(src_idxs) phase.connect(f'parameter_vals:{name}', f'timeseries.all_values:parameters:{name}', src_indices=src_idxs, flat_src_indices=True) - for ts_output in phase._timeseries['timeseries']['outputs']: + for output_name, ts_output in phase._timeseries['timeseries']['outputs'].items(): var = ts_output['name'] - output_name = ts_output['output_name'] units = ts_output['units'] wildcard_units = ts_output['wildcard_units'] shape = ts_output['shape'] if '*' in var: # match outputs from the ODE - matches = filter(list(ode_outputs.keys()), var) + # TODO: is filter still case INSENSTIVE on windows? If so, fix this + matches = filter(ode_outputs.keys(), var) else: matches = [var]