From 133f106d0c2960afe919b816ec433df57206f7f7 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 26 Jul 2023 13:45:49 -0400 Subject: [PATCH 01/20] Timeseries inputs and outputs are now tagged with dymos variable type and constraint type --- ...doc_brachistochrone_polynomial_controls.py | 159 +++--------------- dymos/phase/options.py | 3 + dymos/phase/phase.py | 33 ++-- dymos/phase/test/test_timeseries.py | 16 ++ .../analytic_timeseries_output_comp.py | 8 +- .../explicit_shooting_timeseries_comp.py | 8 +- .../pseudospectral_timeseries_output_comp.py | 8 +- .../components/solve_ivp_timeseries_comp.py | 8 +- dymos/transcriptions/transcription_base.py | 71 ++++++-- 9 files changed, 130 insertions(+), 184 deletions(-) diff --git a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py index d5a87e489..11a880941 100644 --- a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py +++ b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py @@ -538,10 +538,6 @@ def test_brachistochrone_polynomial_control_radau(self): class TestBrachistochronePolynomialControlRatePathConstrained(unittest.TestCase): def test_brachistochrone_polynomial_control_gauss_lobatto(self): - import numpy as np - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -594,42 +590,7 @@ def test_brachistochrone_polynomial_control_gauss_lobatto(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') - - plt.show() - def test_brachistochrone_polynomial_control_radau(self): - import numpy as np - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -648,7 +609,7 @@ def test_brachistochrone_polynomial_control_radau(self): phase.add_state('x', fix_initial=True, fix_final=True) - phase.add_state('y', fix_initial=True, fix_final=True) + phase.add_state('y', fix_initial=True, fix_final=False) phase.add_state('v', fix_initial=True, fix_final=False) @@ -657,6 +618,8 @@ def test_brachistochrone_polynomial_control_radau(self): phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665) phase.add_path_constraint('theta_rate', lower=0, upper=120) + phase.add_path_constraint('y', lower=0, upper=10) + phase.add_boundary_constraint('y', loc='final', equals=5) # Minimize time at the end of the phase phase.add_objective('time', loc='final', scaler=10) @@ -682,46 +645,32 @@ def test_brachistochrone_polynomial_control_radau(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') + for prob in [p, exp_out]: + ts_group = prob.model._get_subsystem('phase0.timeseries') - fig, ax = plt.subplots() + map_type_to_promnames = {'dymos.type:time': {'time'}, + 'dymos.type:t_phase': set(), + 'dymos.type:control': set(), + 'dymos.type:polynomial_control': {'theta'}, + 'dymos.type:polynomial_control_rate': set(), + 'dymos.type:polynomial_control_rate2': set(), + 'dymos.type:parameter': set(), + 'dymos.type:state': {'x', 'y', 'v'}, + 'dymos.initial_boundary_constraint': set(), + 'dymos.final_boundary_constraint': {'y'}, + 'dymos.path_constraint': {'theta_rate', 'y'}} - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') - - plt.show() + for dymos_type, prom_names in map_type_to_promnames.items(): + prom_outputs = {meta['prom_name'] for abs_path, meta in + ts_group.list_outputs(tags=[dymos_type], out_stream=None)} + self.assertSetEqual(prom_outputs, prom_names, + msg=f'\n{dymos_type}\nin outputs: {prom_outputs}\nexpected: {prom_names}') @use_tempdirs class TestBrachistochronePolynomialControlRate2PathConstrained(unittest.TestCase): def test_brachistochrone_polynomial_control_gauss_lobatto(self): - import numpy as np - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -774,42 +723,7 @@ def test_brachistochrone_polynomial_control_gauss_lobatto(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') - - plt.show() - def test_brachistochrone_polynomial_control_radau(self): - import numpy as np - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -862,37 +776,6 @@ def test_brachistochrone_polynomial_control_radau(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') - - plt.show() - @use_tempdirs class TestBrachistochronePolynomialControlSimulation(unittest.TestCase): diff --git a/dymos/phase/options.py b/dymos/phase/options.py index 48fe7cead..d032d9eca 100644 --- a/dymos/phase/options.py +++ b/dymos/phase/options.py @@ -745,6 +745,9 @@ def __init__(self, read_only=False): self.declare(name='expr_kwargs', default={}, allow_none=False, desc='Options to be passed to the timeseries expression comp when adding the expression.') + self.declare(name='tags', default=None, allow_none=True, + desc='Tags to be applied to the timeseries inputs and outputs.') + class PhaseTimeseriesOptionsDictionary(om.OptionsDictionary): """ diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index bb32cb829..c6a96f0fd 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -1182,12 +1182,6 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None, bc['flat_indices'] = flat_indices bc['is_expr'] = is_expr - # 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': - 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, lower=None, upper=None, equals=None, scaler=None, adder=None, ref=None, ref0=None, linear=False, flat_indices=False): @@ -1289,14 +1283,8 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None pc['flat_indices'] = flat_indices pc['is_expr'] = is_expr - # 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': - 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, - timeseries='timeseries', **kwargs): + timeseries='timeseries', tags=None, **kwargs): r""" Add a variable to the timeseries outputs of the phase. @@ -1324,6 +1312,8 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap 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. + tags : str or list of str or None + Any tags to be applied to the timeseries output. **kwargs Additional arguments passed to the exec comp. """ @@ -1342,7 +1332,8 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap shape=shape, timeseries=timeseries, rate=False, - expr=expr) + expr=expr, + tags=tags) # Handle specific units for wildcard names. if oname is not None and '*' in name_i: @@ -1356,7 +1347,8 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap timeseries=timeseries, rate=False, expr=expr, - expr_kwargs=kwargs) + expr_kwargs=kwargs, + tags=tags) def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, timeseries='timeseries'): @@ -1414,7 +1406,7 @@ def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, rate=True) def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, - timeseries='timeseries', rate=False, expr=False, expr_kwargs=None): + timeseries='timeseries', rate=False, expr=False, expr_kwargs=None, tags=None): r""" Add a single variable or rate to the timeseries outputs of the phase. @@ -1444,7 +1436,13 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha 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. - expr : + expr : bool + If True, the given name is an equation to be evaluated by an ExecComp. + expr_kwargs : dict or None + A dictionary of keyword arguments to be passed to the ExecComp when expr = True. + tags : str or list of str or None + A string or list of strings used as tags for the timeseries inputs and outputs. Dymos uses + this to tag the timeseries outputs with the variable type (state, time, control, parameter). Returns ------- @@ -1475,6 +1473,7 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha ts_output['is_rate'] = rate ts_output['is_expr'] = expr ts_output['expr_kwargs'] = expr_kwargs + ts_output['tags'] = [tags] if isinstance(tags, str) else tags self._timeseries[timeseries]['outputs'][output_name] = ts_output diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py index 235f7ee7b..21124c2bd 100644 --- a/dymos/phase/test/test_timeseries.py +++ b/dymos/phase/test/test_timeseries.py @@ -31,6 +31,8 @@ def test_timeseries_gl(self, test_smaller_timeseries=False): phase = dm.Phase(ode_class=BrachistochroneODE, transcription=dm.GaussLobatto(num_segments=8, order=3, compressed=True)) + phase.timeseries_options['include_t_phase'] = True + p.model.add_subsystem('phase0', phase) phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) @@ -53,6 +55,7 @@ def test_timeseries_gl(self, test_smaller_timeseries=False): phase.add_objective('time_phase', loc='final', scaler=10) p.model.linear_solver = om.DirectSolver() + p.setup(check=True) p['phase0.t_initial'] = 0.0 @@ -105,6 +108,19 @@ def test_timeseries_gl(self, test_smaller_timeseries=False): else: # no error accessing timseries.parameter exp_out.get_val(f'phase0.timeseries.{dp}') + ts_group = exp_out.model._get_subsystem('phase0.timeseries') + + map_type_to_promnames = {'dymos.type:time': {'time'}, + 'dymos.type:t_phase': {'time_phase'}, + 'dymos.type:control': {'theta'}, + 'dymos.type:parameter': {'g'}, + 'dymos.type:state': {'x', 'y', 'v'}} + + for dymos_type, prom_names in map_type_to_promnames.items(): + prom_outputs = {meta['prom_name'] for abs_path, meta in + ts_group.list_outputs(tags=[dymos_type], out_stream=None)} + self.assertSetEqual(prom_outputs, prom_names) + def test_timeseries_gl_smaller_timeseries(self): self.test_timeseries_gl(test_smaller_timeseries=True) diff --git a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py index 23dab54d6..9fbfc8f00 100644 --- a/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py +++ b/dymos/transcriptions/analytic/analytic_timeseries_output_comp.py @@ -95,7 +95,7 @@ def setup(self): 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, rate=False): + def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False, tags=None): """ Add a single timeseries output. @@ -117,6 +117,8 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals The src path of the variables input, used to prevent redundant inputs. rate : bool If True, timeseries output is a rate. + tags : list of str or None + The tags to be applied to the inputs and outputs of the timeseries. Returns ------- @@ -139,13 +141,13 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals input_name = f'input_values:{name}' self.add_input(input_name, shape=(input_num_nodes,) + shape, - units=units, desc=desc) + units=units, desc=desc, tags=tags) self._sources[src] = input_name input_units = self._units[input_name] = units 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, tags=tags) self._vars[name] = (input_name, output_name, shape, rate) diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting_timeseries_comp.py b/dymos/transcriptions/explicit_shooting/explicit_shooting_timeseries_comp.py index 786adceb2..568b30852 100644 --- a/dymos/transcriptions/explicit_shooting/explicit_shooting_timeseries_comp.py +++ b/dymos/transcriptions/explicit_shooting/explicit_shooting_timeseries_comp.py @@ -49,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, rate=False): + def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False, tags=None): """ Add a single timeseries output. @@ -71,6 +71,8 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals The src path of the variables input, used to prevent redundant inputs. rate : bool If True, timeseries output is a rate. + tags : list of str or None + The tags to be applied to the inputs and outputs of the timeseries. Returns ------- @@ -97,12 +99,12 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals input_name = f'input_values:{name}' self.add_input(input_name, shape=(input_num_nodes,) + shape, - units=units, desc=desc) + units=units, desc=desc, tags=tags) self._sources[src] = input_name input_units = self._units[input_name] = units added_source = True - self.add_output(name, shape=(output_num_nodes,) + shape, units=units, desc=desc) + self.add_output(name, shape=(output_num_nodes,) + shape, units=units, desc=desc, tags=tags) self._vars[name] = (input_name, name, shape, rate) diff --git a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py index c56a7ed0f..0ee7f42c3 100644 --- a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py +++ b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py @@ -95,7 +95,7 @@ def setup(self): 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, rate=False): + def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False, tags=None): """ Add a single timeseries output. @@ -117,6 +117,8 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals The src path of the variables input, used to prevent redundant inputs. rate : bool If True, timeseries output is a rate. + tags : list of str or None + The tags to be applied to the inputs and outputs of the timeseries. Returns ------- @@ -139,13 +141,13 @@ def _add_output_configure(self, name, units, shape, desc='', src=None, rate=Fals input_name = f'input_values:{name}' self.add_input(input_name, shape=(input_num_nodes,) + shape, - units=units, desc=desc) + units=units, desc=desc, tags=tags) self._sources[src] = input_name input_units = self._units[input_name] = units 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, tags=tags) self._vars[name] = (input_name, output_name, shape, rate) 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 a3b499c25..721889eef 100644 --- a/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py +++ b/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py @@ -93,7 +93,7 @@ def setup(self): self.add_input('dt_dstau', shape=(self.output_num_nodes,), units=self.options['time_units']) - def _add_output_configure(self, name, units, shape, desc, src=None, rate=False): + def _add_output_configure(self, name, units, shape, desc, src=None, rate=False, tags=None): """ Add a single timeseries output. @@ -115,6 +115,8 @@ def _add_output_configure(self, name, units, shape, desc, src=None, rate=False): The source of the timeseries output. rate : bool If True, timeseries output is a rate. + tags : list of str or None + The tags to be applied to the inputs and outputs of the timeseries. """ if rate: om.issue_warning(f'Timeseries rate outputs not currently supported in simulate: {name} being skipped.') @@ -135,13 +137,13 @@ def _add_output_configure(self, name, units, shape, desc, src=None, rate=False): input_name = f'input_values:{name}' self.add_input(input_name, shape=(input_num_nodes,) + shape, - units=units, desc=desc) + units=units, desc=desc, tags=tags) self._sources[src] = input_name input_units = self._units[input_name] = units 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, tags=tags) self._vars[name] = (input_name, output_name, shape, rate) diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py index bda172659..a23d8237c 100644 --- a/dymos/transcriptions/transcription_base.py +++ b/dymos/transcriptions/transcription_base.py @@ -88,10 +88,13 @@ def setup_time(self, phase): for ts_name, ts_options in phase._timeseries.items(): if t_name not in ts_options['outputs']: - phase.add_timeseries_output(t_name, timeseries=ts_name) + phase.add_timeseries_output(t_name, timeseries=ts_name, + tags=['dymos.type:time']) if t_phase_name not in ts_options['outputs'] and \ - (phase.timeseries_options['include_t_phase'] or time_options['time_phase_targets']): - phase.add_timeseries_output(t_phase_name, timeseries=ts_name) + (phase.timeseries_options['include_t_phase'] or + time_options['time_phase_targets'] is not _unspecified): + phase.add_timeseries_output(t_phase_name, timeseries=ts_name, + tags=['dymos.type:t_phase']) def configure_time(self, phase): """ @@ -166,15 +169,17 @@ def setup_controls(self, phase): for ts_name, ts_options in phase._timeseries.items(): if f'{control_prefix}{name}' not in ts_options['outputs']: phase.add_timeseries_output(name, output_name=f'{control_prefix}{name}', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:control']) if f'{control_rate_prefix}{name}_rate' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or options['rate_targets']): + (phase.timeseries_options['include_control_rates'] or + options['rate_targets'] is not _unspecified): phase.add_timeseries_output(f'{name}_rate', output_name=f'{control_rate_prefix}{name}_rate', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:control_rate']) if f'{control_rate_prefix}{name}_rate2' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or options['rate2_targets']): + (phase.timeseries_options['include_control_rates'] or + options['rate2_targets'] is not _unspecified): phase.add_timeseries_output(f'{name}_rate2', output_name=f'{control_rate_prefix}{name}_rate2', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:control_rate2']) def configure_controls(self, phase): """ @@ -214,15 +219,17 @@ def setup_polynomial_controls(self, phase): for ts_name, ts_options in phase._timeseries.items(): if f'{prefix}{name}' not in ts_options['outputs']: phase.add_timeseries_output(name, output_name=f'{prefix}{name}', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:polynomial_control']) if f'{rate_prefix}{name}_rate' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or options['rate_targets']): + (phase.timeseries_options['include_control_rates'] or + options['rate_targets'] is not _unspecified): phase.add_timeseries_output(f'{name}_rate', output_name=f'{rate_prefix}{name}_rate', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:polynomial_control_rate']) if f'{rate_prefix}{name}_rate2' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or options['rate2_targets']): + (phase.timeseries_options['include_control_rates'] or + options['rate2_targets'] is not _unspecified): phase.add_timeseries_output(f'{name}_rate2', output_name=f'{rate_prefix}{name}_rate2', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:polynomial_control_rate2']) def configure_polynomial_controls(self, phase): """ @@ -254,7 +261,7 @@ def setup_parameters(self, phase): for ts_name, ts_options in phase._timeseries.items(): if f'{param_prefix}{name}' not in ts_options['outputs']: phase.add_timeseries_output(name, output_name=f'{param_prefix}{name}', - timeseries=ts_name) + timeseries=ts_name, tags=['dymos.type:parameter']) def configure_parameters(self, phase): """ @@ -350,13 +357,15 @@ def configure_states(self, phase): for ts_name, ts_options in phase._timeseries.items(): if f'{state_prefix}{name}' not in ts_options['outputs']: phase.add_timeseries_output(name, output_name=f'{state_prefix}{name}', - timeseries=ts_name) + timeseries=ts_name, + tags=['dymos.type:state']) if options['rate_source'] and phase.timeseries_options['include_state_rates']: output_name = f'{state_rate_prefix}{name}' if state_rate_prefix else options['rate_source'] if output_name not in ts_options['outputs']: phase.add_timeseries_output(name=options['rate_source'], output_name=output_name, - timeseries=ts_name) + timeseries=ts_name, + tags=['dymos.type:state_rate']) def setup_ode(self, phase): """ @@ -447,13 +456,15 @@ def configure_timeseries_outputs(self, phase): shape = ts_output['shape'] src = ts_output['src'] is_rate = ts_output['is_rate'] + tags = ts_output['tags'] added_src = timeseries_comp._add_output_configure(name, shape=shape, units=units, desc='', src=src, - rate=is_rate) + rate=is_rate, + tags=tags) if added_src: phase.connect(src_name=src, tgt_name=f'{timeseries_name}.input_values:{name}', @@ -475,10 +486,28 @@ def configure_boundary_constraints(self, phase): for ibc in phase._initial_boundary_constraints: con_output, constraint_kwargs = self._get_constraint_kwargs('initial', ibc, phase) phase.add_constraint(con_output, **constraint_kwargs) + con_name = ibc['constraint_name'] + # Automatically add the requested variable to the timeseries outputs if it's an ODE output. + if con_name not in phase._timeseries['timeseries']['outputs']: + phase.add_timeseries_output(con_name, output_name=con_name, + units=ibc['units'], shape=ibc['shape'], + tags=[f'dymos.initial_boundary_constraint']) + else: + phase._timeseries['timeseries']['outputs'][con_name]['tags'] += \ + [f'dymos.initial_boundary_constraint'] for fbc in phase._final_boundary_constraints: con_output, constraint_kwargs = self._get_constraint_kwargs('final', fbc, phase) phase.add_constraint(con_output, **constraint_kwargs) + # Automatically add the requested variable to the timeseries outputs if it's not already. + con_name = fbc['constraint_name'] + if con_name not in phase._timeseries['timeseries']['outputs']: + phase.add_timeseries_output(con_name, output_name=con_name, + units=fbc['units'], shape=fbc['shape'], + tags=[f'dymos.final_boundary_constraint']) + else: + phase._timeseries['timeseries']['outputs'][con_name]['tags'] += \ + [f'dymos.final_boundary_constraint'] def _get_constraint_kwargs(self, constraint_type, options, phase): """ @@ -577,6 +606,14 @@ def configure_path_constraints(self, phase): for pc in phase._path_constraints: con_output, constraint_kwargs = self._get_constraint_kwargs('path', pc, phase) phase.add_constraint(con_output, **constraint_kwargs) + con_name = pc['constraint_name'] + if con_name not in phase._timeseries['timeseries']['outputs']: + phase.add_timeseries_output(con_name, output_name=con_name, + units=pc['units'], shape=pc['shape'], + tags=[f'dymos.path_constraint']) + else: + phase._timeseries['timeseries']['outputs'][con_name]['tags'] += \ + [f'dymos.path_constraint'] def configure_objective(self, phase): """ From 00efc2fa660df76d99558cc5d67a4c2f11178a40 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 26 Jul 2023 14:03:15 -0400 Subject: [PATCH 02/20] reverted addition of constraint type to timeseries tags --- dymos/phase/phase.py | 7 +++++ dymos/transcriptions/transcription_base.py | 26 ------------------- .../timeseries/bokeh_timeseries_report.py | 2 -- 3 files changed, 7 insertions(+), 28 deletions(-) diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index c6a96f0fd..f9f457d4d 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -1182,6 +1182,10 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None, bc['flat_indices'] = flat_indices bc['is_expr'] = is_expr + 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, lower=None, upper=None, equals=None, scaler=None, adder=None, ref=None, ref0=None, linear=False, flat_indices=False): @@ -1283,6 +1287,9 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None pc['flat_indices'] = flat_indices pc['is_expr'] = is_expr + 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, timeseries='timeseries', tags=None, **kwargs): r""" diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py index a23d8237c..5f59b0797 100644 --- a/dymos/transcriptions/transcription_base.py +++ b/dymos/transcriptions/transcription_base.py @@ -486,28 +486,10 @@ def configure_boundary_constraints(self, phase): for ibc in phase._initial_boundary_constraints: con_output, constraint_kwargs = self._get_constraint_kwargs('initial', ibc, phase) phase.add_constraint(con_output, **constraint_kwargs) - con_name = ibc['constraint_name'] - # Automatically add the requested variable to the timeseries outputs if it's an ODE output. - if con_name not in phase._timeseries['timeseries']['outputs']: - phase.add_timeseries_output(con_name, output_name=con_name, - units=ibc['units'], shape=ibc['shape'], - tags=[f'dymos.initial_boundary_constraint']) - else: - phase._timeseries['timeseries']['outputs'][con_name]['tags'] += \ - [f'dymos.initial_boundary_constraint'] for fbc in phase._final_boundary_constraints: con_output, constraint_kwargs = self._get_constraint_kwargs('final', fbc, phase) phase.add_constraint(con_output, **constraint_kwargs) - # Automatically add the requested variable to the timeseries outputs if it's not already. - con_name = fbc['constraint_name'] - if con_name not in phase._timeseries['timeseries']['outputs']: - phase.add_timeseries_output(con_name, output_name=con_name, - units=fbc['units'], shape=fbc['shape'], - tags=[f'dymos.final_boundary_constraint']) - else: - phase._timeseries['timeseries']['outputs'][con_name]['tags'] += \ - [f'dymos.final_boundary_constraint'] def _get_constraint_kwargs(self, constraint_type, options, phase): """ @@ -606,14 +588,6 @@ def configure_path_constraints(self, phase): for pc in phase._path_constraints: con_output, constraint_kwargs = self._get_constraint_kwargs('path', pc, phase) phase.add_constraint(con_output, **constraint_kwargs) - con_name = pc['constraint_name'] - if con_name not in phase._timeseries['timeseries']['outputs']: - phase.add_timeseries_output(con_name, output_name=con_name, - units=pc['units'], shape=pc['shape'], - tags=[f'dymos.path_constraint']) - else: - phase._timeseries['timeseries']['outputs'][con_name]['tags'] += \ - [f'dymos.path_constraint'] def configure_objective(self, phase): """ diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index 175f2b6ca..1bdc68bcc 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -154,8 +154,6 @@ def _load_data_sources(prob, solution_record_file=None, simulation_record_file=N phase_name = phase.pathname.split('.')[-1] data_dict[traj_name]['param_data_by_phase'][phase_name] = {'param': [], 'val': [], 'units': []} - phase_sol_data = data_dict[traj_name]['sol_data_by_phase'][phase_name] = {} - phase_sim_data = data_dict[traj_name]['sim_data_by_phase'][phase_name] = {} ts_units_dict = data_dict[traj_name]['timeseries_units'] param_outputs = {op: meta for op, meta in outputs.items() From 33786891e1332f655304e66bed7e38c49baad6aa Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 26 Jul 2023 16:06:46 -0400 Subject: [PATCH 03/20] added dashed lines representing bounds and circle_crosses representing fixed design variables in the timeseries plot report --- .../test/test_ex_min_time_climb.py | 4 +- .../timeseries/bokeh_timeseries_report.py | 65 ++++++++++++++++++- 2 files changed, 65 insertions(+), 4 deletions(-) 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 1a1401dca..d7b77b58f 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 @@ -117,12 +117,12 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', p['traj.phase0.states:m'] = phase.interp('m', [19030.468, 16841.431]) p['traj.phase0.controls:alpha'] = phase.interp('alpha', [0.0, 0.0]) - dm.run_problem(p, simulate=True) + dm.run_problem(p, simulate=True, make_plots=True) return p -@use_tempdirs +# @use_tempdirs class TestMinTimeClimb(unittest.TestCase): def _test_results(self, p, time_name='time'): diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index 1bdc68bcc..7ab848de2 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -14,6 +14,7 @@ except ImportError: _NO_BOKEH = True +import numpy as np import openmdao.api as om from openmdao.utils.units import conversion_to_base_units import dymos as dm @@ -148,6 +149,7 @@ def _load_data_sources(prob, solution_record_file=None, simulation_record_file=N data_dict[traj_name] = {'param_data_by_phase': {}, 'sol_data_by_phase': {}, 'sim_data_by_phase': {}, + 'bounds_by_phase': {}, 'timeseries_units': {}} for phase in traj.system_iter(include_self=True, recurse=True, typ=dm.Phase): @@ -199,7 +201,6 @@ def _load_data_sources(prob, solution_record_file=None, simulation_record_file=N ts_outputs = {op: meta for op, meta in outputs.items() if op.startswith(f'{phase.pathname}.timeseries')} for output_name in sorted(ts_outputs.keys(), key=str.casefold): - meta = ts_outputs[output_name] prom_name = abs2prom_map['output'][output_name] var_name = prom_name.split('.')[-1] @@ -279,6 +280,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi figures = [] x_range = None + size = 6 for var_name in sorted(ts_units_dict.keys(), key=str.casefold): fig_kwargs = {'x_range': x_range} if x_range is not None else {} @@ -309,9 +311,68 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi if x_name in sol_data and var_name in sol_data: legend_items = [] if sol_data: - sol_plot = fig.circle(x='time', y=var_name, source=sol_source, color=color) + lower = upper = None + fix_initial = False + fix_final = False + opt = True + if var_name in phase.state_options: + lower = phase.state_options[var_name]['lower'] + upper = phase.state_options[var_name]['upper'] + fix_initial = phase.state_options[var_name]['fix_initial'] + fix_final = phase.state_options[var_name]['fix_final'] + elif var_name in phase.control_options: + lower = phase.control_options[var_name]['lower'] + upper = phase.control_options[var_name]['upper'] + fix_initial = phase.control_options[var_name]['fix_initial'] + fix_final = phase.control_options[var_name]['fix_final'] + opt = phase.control_options[var_name]['opt'] + elif var_name in phase.polynomial_control_options: + lower = phase.polynomial_control_options[var_name]['lower'] + upper = phase.polynomial_control_options[var_name]['upper'] + fix_initial = phase.polynomial_control_options[var_name]['fix_initial'] + fix_final = phase.polynomial_control_options[var_name]['fix_final'] + opt = phase.polynomial_control_options[var_name]['opt'] + elif var_name in phase.parameter_options: + lower = phase.parameter_options[var_name]['lower'] + upper = phase.parameter_options[var_name]['upper'] + opt = phase.parameter_options[var_name]['opt'] + + if opt: + sol_plot = fig.circle(x=x_name, y=var_name, source=sol_source, + color=color, size=size) + else: + sol_plot = fig.circle_cross(x=x_data[0], y=sol_data[var_name][0, ...], + color=color, fill_color='white', + size=size+2, line_width=1) + sol_plot.tags.extend(['sol', f'phase:{phase_name}']) legend_items.append(sol_plot) + + # Plot the bounds if available + x_data = sol_data[x_name].ravel() + if lower: + lb_plot = fig.line(x=x_data, y=lower * np.ones_like(x_data), + line_dash="dashed") + lb_plot.tags.extend(['bounds']) + lb_plot.tags.extend(['bounds', f'phase:{phase_name}']) + + if upper: + ub_plot = fig.line(x=x_data, y=upper * np.ones_like(x_data), + line_dash="dashed") + ub_plot.tags.extend(['bounds', f'phase:{phase_name}']) + + if fix_initial: + fix_initial_plot = fig.circle_cross(x=x_data[0], y=sol_data[var_name][0, ...], + color=color, fill_color='white', + size=size+2, line_width=2) + fix_initial_plot.tags.extend(['sol', f'phase:{phase_name}']) + + if fix_final: + fix_final_plot = fig.circle_cross(x=x_data[-1], y=sol_data[var_name][-1, ...], + color=color, fill_color='white', + size=size+2, line_width=2) + fix_final_plot.tags.extend(['sol', f'phase:{phase_name}']) + if sim_data: sim_plot = fig.line(x='time', y=var_name, source=sim_source, color=color) sim_plot.tags.extend(['sim', f'phase:{phase_name}']) From 14a28af12f81482ce3075c9c09eb17f0f04d580c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 26 Jul 2023 16:34:21 -0400 Subject: [PATCH 04/20] cleanup 1 --- dymos/phase/test/test_timeseries.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py index 21124c2bd..541a6b06e 100644 --- a/dymos/phase/test/test_timeseries.py +++ b/dymos/phase/test/test_timeseries.py @@ -113,7 +113,7 @@ def test_timeseries_gl(self, test_smaller_timeseries=False): map_type_to_promnames = {'dymos.type:time': {'time'}, 'dymos.type:t_phase': {'time_phase'}, 'dymos.type:control': {'theta'}, - 'dymos.type:parameter': {'g'}, + 'dymos.type:parameter': set() if test_smaller_timeseries else {'g'}, 'dymos.type:state': {'x', 'y', 'v'}} for dymos_type, prom_names in map_type_to_promnames.items(): @@ -158,6 +158,7 @@ def test_timeseries_radau(self, test_smaller_timeseries=False): phase.add_objective('time_phase', loc='final', scaler=10) phase.timeseries_options['include_state_rates'] = True + phase.timeseries_options['include_t_phase'] = True p.model.options['assembled_jac_type'] = 'csc' p.model.linear_solver = om.DirectSolver() From 8480d650e894924f9178d92b721c72ac3f33a983 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 26 Jul 2023 16:38:10 -0400 Subject: [PATCH 05/20] tests for some control rate issues --- .../test/test_brachistochrone_control_rate_targets.py | 2 ++ .../test/test_brachistochrone_vector_path_constraints.py | 5 +++++ 2 files changed, 7 insertions(+) diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py index be026f048..5e4d13146 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py @@ -126,6 +126,7 @@ def test_brachistochrone_control_rate_targets_gauss_lobatto(self): p.model.add_subsystem('phase0', phase) + phase.timeseries_options['include_control_rates'] = True phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) phase.add_state('x', rate_source='xdot', @@ -411,6 +412,7 @@ def test_brachistochrone_control_rate_targets_radau(self): phase = dm.Phase(ode_class=BrachistochroneRateTargetODE, transcription=dm.Radau(num_segments=10)) + phase.timeseries_options['include_control_rates'] = True p.model.add_subsystem('phase0', phase) phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_vector_path_constraints.py b/dymos/examples/brachistochrone/test/test_brachistochrone_vector_path_constraints.py index 70e77ce5c..e9a84e6ff 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_vector_path_constraints.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_vector_path_constraints.py @@ -32,6 +32,7 @@ def test_brachistochrone_vector_state_path_constraints_radau_partial_indices(sel p.model.add_subsystem('phase0', phase) + phase.timeseries_options['include_control_rates'] = True phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) phase.add_state('pos', fix_initial=True, fix_final=True) @@ -144,6 +145,7 @@ def test_brachistochrone_vector_ode_path_constraints_radau_partial_indices(self) phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE, transcription=dm.Radau(num_segments=20, order=3)) + phase.timeseries_options['include_control_rates'] = True p.model.add_subsystem('phase0', phase) phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) @@ -257,6 +259,7 @@ def test_brachistochrone_vector_ode_path_constraints_radau_no_indices(self): transcription=dm.Radau(num_segments=20, order=3)) p.model.add_subsystem('phase0', phase) + phase.timeseries_options['include_control_rates'] = True phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) @@ -369,6 +372,7 @@ def test_brachistochrone_vector_state_path_constraints_gl_partial_indices(self): phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE, transcription=dm.GaussLobatto(num_segments=20, order=3)) + phase.timeseries_options['include_control_rates'] = True p.model.add_subsystem('phase0', phase) phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) @@ -598,6 +602,7 @@ def test_brachistochrone_vector_ode_path_constraints_gl_no_indices(self): phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE, transcription=dm.GaussLobatto(num_segments=20, order=3)) + phase.timeseries_options['include_control_rates'] = True p.model.add_subsystem('phase0', phase) phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) From cca2e04085add30e132b6fc0f71f8c1ca83833f9 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 7 Aug 2023 13:14:49 -0400 Subject: [PATCH 06/20] add_timeseries_output now updates tags of the given output name, if provided. Changed some tests to look for control rates in the control group since they are no longer added to the timeseries outputs by default. --- ...st_brachistochrone_control_rate_targets.py | 8 +- ...histochrone_vector_boundary_constraints.py | 158 +++++------------- ...doc_brachistochrone_polynomial_controls.py | 4 +- dymos/phase/options.py | 2 +- dymos/phase/phase.py | 27 ++- dymos/phase/test/test_phase.py | 15 -- dymos/trajectory/test/test_trajectory.py | 23 ++- dymos/trajectory/trajectory.py | 8 +- dymos/transcriptions/transcription_base.py | 60 +++---- 9 files changed, 111 insertions(+), 194 deletions(-) diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py index 5e4d13146..f0d6bd77f 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py @@ -575,9 +575,9 @@ def test_brachistochrone_polynomial_control_rate_targets_gauss_lobatto(self): fig, ax = plt.subplots() t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') + theta_imp = p.get_val('phase0.polynomial_control_rates:theta_rate') t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') + theta_exp = exp_out.get_val('phase0.polynomial_control_rates:theta_rate') ax.plot(t_imp, theta_imp, 'ro', label='solution') ax.plot(t_exp, theta_exp, 'b-', label='simulated') @@ -669,9 +669,9 @@ def test_brachistochrone_polynomial_control_rate_targets_radau(self): fig, ax = plt.subplots() t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') + theta_imp = p.get_val('phase0.polynomial_control_rates:theta_rate') t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') + theta_exp = exp_out.get_val('phase0.polynomial_control_rates:theta_rate') ax.plot(t_imp, theta_imp, 'ro', label='solution') ax.plot(t_exp, theta_exp, 'b-', label='simulated') diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_vector_boundary_constraints.py b/dymos/examples/brachistochrone/test/test_brachistochrone_vector_boundary_constraints.py index 37db29d46..9d1c27117 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_vector_boundary_constraints.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_vector_boundary_constraints.py @@ -1,8 +1,5 @@ import unittest -import matplotlib -import matplotlib.pyplot as plt - import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal from openmdao.utils.testing_utils import use_tempdirs @@ -10,10 +7,10 @@ import dymos as dm from dymos.examples.brachistochrone.brachistochrone_vector_states_ode \ import BrachistochroneVectorStatesODE +from dymos.utils.testing_utils import assert_timeseries_near_equal SHOW_PLOTS = True -matplotlib.use('Agg') @use_tempdirs @@ -65,46 +62,25 @@ def test_brachistochrone_vector_boundary_constraints_radau_no_indices(self): assert_near_equal(p.get_val('phase0.t')[-1], 1.8016, tolerance=1.0E-3) - # Plot results - if SHOW_PLOTS: - p.run_driver() - exp_out = phase.simulate(times_per_seg=10) - - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.pos')[:, 0] - y_imp = p.get_val('phase0.timeseries.pos')[:, 1] - - x_exp = exp_out.get_val('phase0.timeseries.pos')[:, 0] - y_exp = exp_out.get_val('phase0.timeseries.pos')[:, 1] - - ax.plot(x_imp, y_imp, 'ro', label='implicit') - ax.plot(x_exp, y_exp, 'b-', label='explicit') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') + exp_out = phase.simulate(times_per_seg=10) - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') + x_imp = p.get_val('phase0.timeseries.time') + y_imp = p.get_val('phase0.control_rates:theta_rate2') - x_imp = p.get_val('phase0.timeseries.time') - y_imp = p.get_val('phase0.timeseries.theta_rate2') + x_exp = exp_out.get_val('phase0.timeseries.time') + y_exp = exp_out.get_val('phase0.control_rates:theta_rate2') - x_exp = exp_out.get_val('phase0.timeseries.time') - y_exp = exp_out.get_val('phase0.timeseries.theta_rate2') + igd = phase.options['transcription'].grid_data + egd = exp_out.model._get_subsystem('phase0').options['transcription'].options['output_grid'] - ax.plot(x_imp, y_imp, 'ro', label='implicit') - ax.plot(x_exp, y_exp, 'b-', label='explicit') - - ax.set_xlabel('time (s)') - ax.set_ylabel('theta rate2 (rad/s**2)') - ax.grid(True) - ax.legend(loc='lower right') - - plt.show() + for row in range(igd.num_segments): + istart, iend = igd.segment_indices[row, :] + estart, eend = egd.segment_indices[row, :] + assert_timeseries_near_equal(t_ref=x_imp[istart:iend, :], + x_ref=y_imp[istart:iend, :], + t_check=x_exp[estart:eend, :], + x_check=y_exp[estart:eend, :], + rel_tolerance=1.0E-9) return p @@ -153,46 +129,25 @@ def test_brachistochrone_vector_boundary_constraints_radau_full_indices(self): assert_near_equal(p.get_val('phase0.t')[-1], 1.8016, tolerance=1.0E-3) - # Plot results - if SHOW_PLOTS: - p.run_driver() - exp_out = phase.simulate(times_per_seg=10) - - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.pos')[:, 0] - y_imp = p.get_val('phase0.timeseries.pos')[:, 1] - - x_exp = exp_out.get_val('phase0.timeseries.pos')[:, 0] - y_exp = exp_out.get_val('phase0.timeseries.pos')[:, 1] + exp_out = phase.simulate(times_per_seg=10) - ax.plot(x_imp, y_imp, 'ro', label='implicit') - ax.plot(x_exp, y_exp, 'b-', label='explicit') + x_imp = p.get_val('phase0.timeseries.time') + y_imp = p.get_val('phase0.control_rates:theta_rate2') - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') + x_exp = exp_out.get_val('phase0.timeseries.time') + y_exp = exp_out.get_val('phase0.control_rates:theta_rate2') - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') + igd = phase.options['transcription'].grid_data + egd = exp_out.model._get_subsystem('phase0').options['transcription'].options['output_grid'] - x_imp = p.get_val('phase0.timeseries.time') - y_imp = p.get_val('phase0.timeseries.theta_rate2') - - x_exp = exp_out.get_val('phase0.timeseries.time') - y_exp = exp_out.get_val('phase0.timeseries.theta_rate2') - - ax.plot(x_imp, y_imp, 'ro', label='implicit') - ax.plot(x_exp, y_exp, 'b-', label='explicit') - - ax.set_xlabel('time (s)') - ax.set_ylabel('theta rate2 (rad/s**2)') - ax.grid(True) - ax.legend(loc='lower right') - - plt.show() + for row in range(igd.num_segments): + istart, iend = igd.segment_indices[row, :] + estart, eend = egd.segment_indices[row, :] + assert_timeseries_near_equal(t_ref=x_imp[istart:iend, :], + x_ref=y_imp[istart:iend, :], + t_check=x_exp[estart:eend, :], + x_check=y_exp[estart:eend, :], + rel_tolerance=1.0E-9) return p @@ -239,46 +194,25 @@ def test_brachistochrone_vector_boundary_constraints_radau_partial_indices(self) assert_near_equal(p.get_val('phase0.t')[-1], 1.8016, tolerance=1.0E-3) - # Plot results - if SHOW_PLOTS: - p.run_driver() - exp_out = phase.simulate(times_per_seg=20) - - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.pos')[:, 0] - y_imp = p.get_val('phase0.timeseries.pos')[:, 1] - - x_exp = exp_out.get_val('phase0.timeseries.pos')[:, 0] - y_exp = exp_out.get_val('phase0.timeseries.pos')[:, 1] - - ax.plot(x_imp, y_imp, 'ro', label='implicit') - ax.plot(x_exp, y_exp, 'b-', label='explicit') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.time') - y_imp = p.get_val('phase0.timeseries.theta_rate2') + exp_out = phase.simulate(times_per_seg=10) - x_exp = exp_out.get_val('phase0.timeseries.time') - y_exp = exp_out.get_val('phase0.timeseries.theta_rate2') + x_imp = p.get_val('phase0.timeseries.time') + y_imp = p.get_val('phase0.control_rates:theta_rate2') - ax.plot(x_imp, y_imp, 'ro', label='implicit') - ax.plot(x_exp, y_exp, 'b-', label='explicit') + x_exp = exp_out.get_val('phase0.timeseries.time') + y_exp = exp_out.get_val('phase0.control_rates:theta_rate2') - ax.set_xlabel('time (s)') - ax.set_ylabel('theta rate2 (rad/s**2)') - ax.grid(True) - ax.legend(loc='lower right') + igd = phase.options['transcription'].grid_data + egd = exp_out.model._get_subsystem('phase0').options['transcription'].options['output_grid'] - plt.show() + for row in range(igd.num_segments): + istart, iend = igd.segment_indices[row, :] + estart, eend = egd.segment_indices[row, :] + assert_timeseries_near_equal(t_ref=x_imp[istart:iend, :], + x_ref=y_imp[istart:iend, :], + t_check=x_exp[estart:eend, :], + x_check=y_exp[estart:eend, :], + rel_tolerance=1.0E-9) return p diff --git a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py index 11a880941..419601d86 100644 --- a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py +++ b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py @@ -648,6 +648,8 @@ def test_brachistochrone_polynomial_control_radau(self): for prob in [p, exp_out]: ts_group = prob.model._get_subsystem('phase0.timeseries') + ts_group.list_outputs() + map_type_to_promnames = {'dymos.type:time': {'time'}, 'dymos.type:t_phase': set(), 'dymos.type:control': set(), @@ -657,7 +659,7 @@ def test_brachistochrone_polynomial_control_radau(self): 'dymos.type:parameter': set(), 'dymos.type:state': {'x', 'y', 'v'}, 'dymos.initial_boundary_constraint': set(), - 'dymos.final_boundary_constraint': {'y'}, + 'dymos.final_boundary_constraint': set('y'), 'dymos.path_constraint': {'theta_rate', 'y'}} for dymos_type, prom_names in map_type_to_promnames.items(): diff --git a/dymos/phase/options.py b/dymos/phase/options.py index d032d9eca..3bef8b21f 100644 --- a/dymos/phase/options.py +++ b/dymos/phase/options.py @@ -745,7 +745,7 @@ def __init__(self, read_only=False): self.declare(name='expr_kwargs', default={}, allow_none=False, desc='Options to be passed to the timeseries expression comp when adding the expression.') - self.declare(name='tags', default=None, allow_none=True, + self.declare(name='tags', default=[], allow_none=False, desc='Tags to be applied to the timeseries inputs and outputs.') diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index f9f457d4d..d96a480a7 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -1182,8 +1182,9 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None, bc['flat_indices'] = flat_indices bc['is_expr'] = is_expr - if constraint_name not in self._timeseries['timeseries']['outputs']: - self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape) + print('adding timeseries output', name, constraint_name) + self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape, + tags=[f'dymos.{loc}_boundary_constraint']) def add_path_constraint(self, name, constraint_name=None, units=None, shape=None, indices=None, @@ -1287,8 +1288,8 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None pc['flat_indices'] = flat_indices pc['is_expr'] = is_expr - if constraint_name not in self._timeseries['timeseries']['outputs']: - self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape) + self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape, + tags=['dymos.path_constraint']) def add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, timeseries='timeseries', tags=None, **kwargs): @@ -1324,6 +1325,13 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap **kwargs Additional arguments passed to the exec comp. """ + if tags is None: + _tags = [] + elif isinstance(tags, str): + _tags = [tags] + else: + _tags = tags + if type(name) is list: for i, name_i in enumerate(name): expr = True if '=' in name_i else False @@ -1340,7 +1348,7 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap timeseries=timeseries, rate=False, expr=expr, - tags=tags) + tags=_tags) # Handle specific units for wildcard names. if oname is not None and '*' in name_i: @@ -1355,7 +1363,7 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap rate=False, expr=expr, expr_kwargs=kwargs, - tags=tags) + tags=_tags) def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified, timeseries='timeseries'): @@ -1480,11 +1488,14 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha ts_output['is_rate'] = rate ts_output['is_expr'] = expr ts_output['expr_kwargs'] = expr_kwargs - ts_output['tags'] = [tags] if isinstance(tags, str) else tags + ts_output['tags'] = tags self._timeseries[timeseries]['outputs'][output_name] = ts_output - return output_name + elif tags is not None: + # Output already present, update tags + self._timeseries[timeseries]['outputs'][output_name]['tags'].extend(tags) + def add_timeseries(self, name, transcription, subset='all'): r""" diff --git a/dymos/phase/test/test_phase.py b/dymos/phase/test/test_phase.py index a9cd30bc4..d079a3354 100644 --- a/dymos/phase/test/test_phase.py +++ b/dymos/phase/test/test_phase.py @@ -472,21 +472,6 @@ def test_control_rate2_boundary_constraint_gl(self): p.run_driver() - plt.plot(p.get_val('phase0.timeseries.x'), - p.get_val('phase0.timeseries.y'), 'ko') - - plt.figure() - - plt.plot(p.get_val('phase0.timeseries.time'), - p.get_val('phase0.timeseries.theta'), 'ro') - - plt.plot(p.get_val('phase0.timeseries.time'), - p.get_val('phase0.timeseries.theta_rate'), 'bo') - - plt.plot(p.get_val('phase0.timeseries.time'), - p.get_val('phase0.timeseries.theta_rate2'), 'go') - plt.show() - assert_near_equal(p.get_val('phase0.timeseries.theta_rate2')[-1], 0, tolerance=1.0E-6) diff --git a/dymos/trajectory/test/test_trajectory.py b/dymos/trajectory/test/test_trajectory.py index 6a55e98f0..5e5875352 100644 --- a/dymos/trajectory/test/test_trajectory.py +++ b/dymos/trajectory/test/test_trajectory.py @@ -627,9 +627,6 @@ def test_linked_control_rate(self): self.traj.link_phases(phases=['burn1', 'burn2'], vars=['accel']) self.traj.link_phases(phases=['burn1', 'burn2'], vars=['u1_rate']) - for phs in burn1, coast, burn2: - phs.timeseries_options['include_control_rates'] = True - # Finish Problem Setup p.model.linear_solver = om.DirectSolver() @@ -674,8 +671,10 @@ def test_linked_control_rate(self): p.run_model() - burn1_u1_final = p.get_val('burn1.timeseries.u1_rate')[-1, ...] - burn2_u1_initial = p.get_val('burn2.timeseries.u1_rate')[0, ...] + p.model.list_outputs() + + burn1_u1_final = p.get_val('burn1.control_rates:u1_rate')[-1, ...] + burn2_u1_initial = p.get_val('burn2.control_rates:u1_rate')[0, ...] u1_linkage_error = p.get_val('linkages.burn1:u1_rate_final|burn2:u1_rate_initial') assert_near_equal(u1_linkage_error, burn1_u1_final - burn2_u1_initial) @@ -804,8 +803,8 @@ def test_linked_control_rate2(self): p.run_model() - burn1_u1_final = p.get_val('burn1.timeseries.u1_rate2')[-1, ...] - burn2_u1_initial = p.get_val('burn2.timeseries.u1_rate2')[0, ...] + burn1_u1_final = p.get_val('burn1.control_rates:u1_rate2')[-1, ...] + burn2_u1_initial = p.get_val('burn2.control_rates:u1_rate2')[0, ...] u1_linkage_error = p.get_val('linkages.burn1:u1_rate2_final|burn2:u1_rate2_initial') assert_near_equal(u1_linkage_error, burn1_u1_final - burn2_u1_initial) @@ -934,8 +933,8 @@ def test_linked_polynomial_control_rate(self): p.run_model() - burn1_u1_final = p.get_val('burn1.timeseries.u1_rate')[-1, ...] - burn2_u1_initial = p.get_val('burn2.timeseries.u1_rate')[0, ...] + burn1_u1_final = p.get_val('burn1.control_rates:u1_rate')[-1, ...] + burn2_u1_initial = p.get_val('burn2.polynomial_control_rates:u1_rate')[0, ...] u1_linkage_error = p.get_val('linkages.burn1:u1_rate_final|burn2:u1_rate_initial') assert_near_equal(u1_linkage_error, burn1_u1_final - burn2_u1_initial) @@ -1063,8 +1062,8 @@ def test_linked_polynomial_control_rate2(self): p.run_model() - burn1_u1_final = p.get_val('burn1.timeseries.u1_rate2')[-1, ...] - burn2_u1_initial = p.get_val('burn2.timeseries.u1_rate2')[0, ...] + burn1_u1_final = p.get_val('burn1.control_rates:u1_rate2')[-1, ...] + burn2_u1_initial = p.get_val('burn2.polynomial_control_rates:u1_rate2')[0, ...] u1_linkage_error = p.get_val('linkages.burn1:u1_rate2_final|burn2:u1_rate2_initial') assert_near_equal(u1_linkage_error, burn1_u1_final - burn2_u1_initial) @@ -1884,8 +1883,6 @@ def compute(self, inputs, outputs): ##################################################### dm.run_problem(p, run_driver=False, simulate=False) - p.model.list_inputs() - ascent_h_m = p.get_val('traj.ascent.timeseries.h', units='m')[[0, -1], ...] descent_h_m = p.get_val('traj.descent.timeseries.h', units='m')[[0, -1], ...] diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py index 62de5b655..ffab7c13f 100644 --- a/dymos/trajectory/trajectory.py +++ b/dymos/trajectory/trajectory.py @@ -573,8 +573,8 @@ def _update_linkage_options_configure(self, linkage_options): units[i] = phases[i].control_options[vars[i]]['units'] shapes[i] = phases[i].control_options[vars[i]]['shape'] elif classes[i] in {'control_rate', 'control_rate2'}: - prefix = 'control_rates:' if use_prefix[i] else '' - sources[i] = f'timeseries.{prefix}{vars[i]}' + prefix = 'control_rates:' + sources[i] = f'{prefix}{vars[i]}' control_name = vars[i][:-5] if classes[i] == 'control_rate' else vars[i][:-6] units[i] = phases[i].control_options[control_name]['units'] deriv = 1 if classes[i].endswith('rate') else 2 @@ -586,8 +586,8 @@ def _update_linkage_options_configure(self, linkage_options): units[i] = phases[i].polynomial_control_options[vars[i]]['units'] shapes[i] = phases[i].polynomial_control_options[vars[i]]['shape'] elif classes[i] in {'polynomial_control_rate', 'polynomial_control_rate2'}: - prefix = 'polynomial_control_rates:' if use_prefix[i] else '' - sources[i] = f'timeseries.{prefix}{vars[i]}' + prefix = 'polynomial_control_rates:' + sources[i] = f'{prefix}{vars[i]}' control_name = vars[i][:-5] if classes[i] == 'polynomial_control_rate' else vars[i][:-6] control_units = phases[i].polynomial_control_options[control_name]['units'] time_units = phases[i].time_options['units'] diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py index 5f59b0797..68198a04e 100644 --- a/dymos/transcriptions/transcription_base.py +++ b/dymos/transcriptions/transcription_base.py @@ -87,12 +87,9 @@ def setup_time(self, phase): promotes_inputs=['*'], promotes_outputs=['*']) for ts_name, ts_options in phase._timeseries.items(): - if t_name not in ts_options['outputs']: - phase.add_timeseries_output(t_name, timeseries=ts_name, - tags=['dymos.type:time']) - if t_phase_name not in ts_options['outputs'] and \ - (phase.timeseries_options['include_t_phase'] or - time_options['time_phase_targets'] is not _unspecified): + phase.add_timeseries_output(t_name, timeseries=ts_name, tags=['dymos.type:time']) + if (phase.timeseries_options['include_t_phase'] or + time_options['time_phase_targets'] is not _unspecified): phase.add_timeseries_output(t_phase_name, timeseries=ts_name, tags=['dymos.type:t_phase']) @@ -167,17 +164,14 @@ def setup_controls(self, phase): for name, options in phase.control_options.items(): for ts_name, ts_options in phase._timeseries.items(): - if f'{control_prefix}{name}' not in ts_options['outputs']: - phase.add_timeseries_output(name, output_name=f'{control_prefix}{name}', - timeseries=ts_name, tags=['dymos.type:control']) - if f'{control_rate_prefix}{name}_rate' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or - options['rate_targets'] is not _unspecified): + phase.add_timeseries_output(name, output_name=f'{control_prefix}{name}', + timeseries=ts_name, tags=['dymos.type:control']) + if phase.timeseries_options['include_control_rates'] \ + or options['rate_targets'] is not _unspecified: phase.add_timeseries_output(f'{name}_rate', output_name=f'{control_rate_prefix}{name}_rate', timeseries=ts_name, tags=['dymos.type:control_rate']) - if f'{control_rate_prefix}{name}_rate2' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or - options['rate2_targets'] is not _unspecified): + if phase.timeseries_options['include_control_rates'] \ + or options['rate2_targets'] is not _unspecified: phase.add_timeseries_output(f'{name}_rate2', output_name=f'{control_rate_prefix}{name}_rate2', timeseries=ts_name, tags=['dymos.type:control_rate2']) @@ -217,17 +211,14 @@ def setup_polynomial_controls(self, phase): for name, options in phase.polynomial_control_options.items(): for ts_name, ts_options in phase._timeseries.items(): - if f'{prefix}{name}' not in ts_options['outputs']: - phase.add_timeseries_output(name, output_name=f'{prefix}{name}', - timeseries=ts_name, tags=['dymos.type:polynomial_control']) - if f'{rate_prefix}{name}_rate' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or - options['rate_targets'] is not _unspecified): + phase.add_timeseries_output(name, output_name=f'{prefix}{name}', + timeseries=ts_name, tags=['dymos.type:polynomial_control']) + if phase.timeseries_options['include_control_rates'] or \ + options['rate_targets'] is not _unspecified: phase.add_timeseries_output(f'{name}_rate', output_name=f'{rate_prefix}{name}_rate', timeseries=ts_name, tags=['dymos.type:polynomial_control_rate']) - if f'{rate_prefix}{name}_rate2' not in ts_options['outputs'] and \ - (phase.timeseries_options['include_control_rates'] or - options['rate2_targets'] is not _unspecified): + if phase.timeseries_options['include_control_rates'] or \ + options['rate2_targets'] is not _unspecified: phase.add_timeseries_output(f'{name}_rate2', output_name=f'{rate_prefix}{name}_rate2', timeseries=ts_name, tags=['dymos.type:polynomial_control_rate2']) @@ -259,9 +250,8 @@ def setup_parameters(self, phase): for name, options in phase.parameter_options.items(): if (options['include_timeseries'] is None and include_params) or options['include_timeseries']: for ts_name, ts_options in phase._timeseries.items(): - if f'{param_prefix}{name}' not in ts_options['outputs']: - phase.add_timeseries_output(name, output_name=f'{param_prefix}{name}', - timeseries=ts_name, tags=['dymos.type:parameter']) + phase.add_timeseries_output(name, output_name=f'{param_prefix}{name}', + timeseries=ts_name, tags=['dymos.type:parameter']) def configure_parameters(self, phase): """ @@ -355,17 +345,15 @@ def configure_states(self, phase): for name, options in phase.state_options.items(): for ts_name, ts_options in phase._timeseries.items(): - if f'{state_prefix}{name}' not in ts_options['outputs']: - phase.add_timeseries_output(name, output_name=f'{state_prefix}{name}', - timeseries=ts_name, - tags=['dymos.type:state']) + phase.add_timeseries_output(name, output_name=f'{state_prefix}{name}', + timeseries=ts_name, + tags=['dymos.type:state']) if options['rate_source'] and phase.timeseries_options['include_state_rates']: output_name = f'{state_rate_prefix}{name}' if state_rate_prefix else options['rate_source'] - if output_name not in ts_options['outputs']: - phase.add_timeseries_output(name=options['rate_source'], - output_name=output_name, - timeseries=ts_name, - tags=['dymos.type:state_rate']) + phase.add_timeseries_output(name=options['rate_source'], + output_name=output_name, + timeseries=ts_name, + tags=['dymos.type:state_rate']) def setup_ode(self, phase): """ From 8c2d0e669288267e90fec91865eaaa3db46c0cda Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 7 Aug 2023 13:18:22 -0400 Subject: [PATCH 07/20] pep8 --- dymos/phase/phase.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index d96a480a7..6774710fb 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -1186,7 +1186,6 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None, self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape, tags=[f'dymos.{loc}_boundary_constraint']) - def add_path_constraint(self, name, constraint_name=None, units=None, shape=None, indices=None, lower=None, upper=None, equals=None, scaler=None, adder=None, ref=None, ref0=None, linear=False, flat_indices=False): @@ -1496,7 +1495,6 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha # Output already present, update tags self._timeseries[timeseries]['outputs'][output_name]['tags'].extend(tags) - def add_timeseries(self, name, transcription, subset='all'): r""" Adds a new timeseries output upon which outputs can be provided. From 6e953d80e44f7cc19d22919db114156d88f7dd93 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 7 Aug 2023 14:59:28 -0400 Subject: [PATCH 08/20] more test cleanup --- ...st_brachistochrone_control_rate_targets.py | 193 +++--------------- ...doc_brachistochrone_polynomial_controls.py | 31 ++- .../test/test_state_rate_introspection.py | 3 +- dymos/phase/phase.py | 2 +- dymos/phase/test/test_phase.py | 17 -- 5 files changed, 66 insertions(+), 180 deletions(-) diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py index f0d6bd77f..fabc56457 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py @@ -4,6 +4,10 @@ import openmdao.api as om from openmdao.utils.testing_utils import use_tempdirs +from openmdao.utils.assert_utils import assert_near_equal + +import dymos as dm +from dymos.utils.testing_utils import assert_timeseries_near_equal class BrachistochroneRateTargetODE(om.ExplicitComponent): @@ -286,9 +290,9 @@ def test_brachistochrone_control_rate_targets_radau(self): fig, ax = plt.subplots() t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') + theta_imp = p.get_val('phase0.control_rates:theta_rate') t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') + theta_exp = exp_out.get_val('phase0.control_rates:theta_rate') ax.plot(t_imp, theta_imp, 'ro', label='solution') ax.plot(t_exp, theta_exp, 'b-', label='simulated') @@ -305,12 +309,6 @@ def test_brachistochrone_control_rate_targets_radau(self): class TestBrachistochroneExplicitControlRateTargets(unittest.TestCase): def test_brachistochrone_control_rate_targets_gauss_lobatto(self): - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - p = om.Problem(model=om.Group()) p.driver = om.ScipyOptimizeDriver() p.driver.declare_coloring() @@ -363,44 +361,25 @@ def test_brachistochrone_control_rate_targets_gauss_lobatto(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') - t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') + x_imp = p.get_val('phase0.timeseries.time') + y_imp = p.get_val('phase0.control_rates:theta_rate2') - ax.plot(t_imp, theta_imp, 'ro', label='solution') - ax.plot(t_exp, theta_exp, 'b-', label='simulated') + x_exp = exp_out.get_val('phase0.timeseries.time') + y_exp = exp_out.get_val('phase0.control_rates:theta_rate2') - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') + igd = phase.options['transcription'].grid_data + egd = exp_out.model._get_subsystem('phase0').options['transcription'].options['output_grid'] - plt.show() + for row in range(igd.num_segments): + istart, iend = igd.segment_indices[row, :] + estart, eend = egd.segment_indices[row, :] + assert_timeseries_near_equal(t_ref=x_imp[istart:iend, :], + x_ref=y_imp[istart:iend, :], + t_check=x_exp[estart:eend, :], + x_check=y_exp[estart:eend, :], + rel_tolerance=1.0E-9) def test_brachistochrone_control_rate_targets_radau(self): - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - import matplotlib.pyplot as plt import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -458,40 +437,23 @@ def test_brachistochrone_control_rate_targets_radau(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') - t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - ax.plot(t_exp, theta_exp, 'b-', label='simulated') + x_imp = p.get_val('phase0.timeseries.time') + y_imp = p.get_val('phase0.control_rates:theta_rate2') - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') + x_exp = exp_out.get_val('phase0.timeseries.time') + y_exp = exp_out.get_val('phase0.control_rates:theta_rate2') - plt.show() + igd = phase.options['transcription'].grid_data + egd = exp_out.model._get_subsystem('phase0').options['transcription'].options['output_grid'] + for row in range(igd.num_segments): + istart, iend = igd.segment_indices[row, :] + estart, eend = egd.segment_indices[row, :] + assert_timeseries_near_equal(t_ref=x_imp[istart:iend, :], + x_ref=y_imp[istart:iend, :], + t_check=x_exp[estart:eend, :], + x_check=y_exp[estart:eend, :], + rel_tolerance=1.0E-9) @use_tempdirs class TestBrachistochronePolynomialControlRateTargets(unittest.TestCase): @@ -688,12 +650,6 @@ def test_brachistochrone_polynomial_control_rate_targets_radau(self): class TestBrachistochronePolynomialControlExplicitRateTargets(unittest.TestCase): def test_brachistochrone_polynomial_control_rate_targets_gauss_lobatto(self): - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - p = om.Problem(model=om.Group()) p.driver = om.ScipyOptimizeDriver() p.driver.declare_coloring() @@ -743,51 +699,7 @@ def test_brachistochrone_polynomial_control_rate_targets_gauss_lobatto(self): # Test the results assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3) - # Generate the explicitly simulated trajectory - exp_out = phase.simulate() - - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') - t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - ax.plot(t_exp, theta_exp, 'b-', label='simulated') - - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') - - plt.show() - def test_brachistochrone_polynomial_control_rate_targets_radau(self): - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - import matplotlib.pyplot as plt - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - p = om.Problem(model=om.Group()) p.driver = om.ScipyOptimizeDriver() p.driver.declare_coloring() @@ -837,43 +749,6 @@ def test_brachistochrone_polynomial_control_rate_targets_radau(self): # Test the results assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3) - # Generate the explicitly simulated trajectory - exp_out = phase.simulate() - - fig, ax = plt.subplots() - fig.suptitle('Brachistochrone Solution') - - x_imp = p.get_val('phase0.timeseries.x') - y_imp = p.get_val('phase0.timeseries.y') - - x_exp = exp_out.get_val('phase0.timeseries.x') - y_exp = exp_out.get_val('phase0.timeseries.y') - - ax.plot(x_imp, y_imp, 'ro', label='solution') - ax.plot(x_exp, y_exp, 'b-', label='simulated') - - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.grid(True) - ax.legend(loc='upper right') - - fig, ax = plt.subplots() - - t_imp = p.get_val('phase0.timeseries.time') - theta_imp = p.get_val('phase0.timeseries.theta_rate') - t_exp = exp_out.get_val('phase0.timeseries.time') - theta_exp = exp_out.get_val('phase0.timeseries.theta_rate') - - ax.plot(t_imp, theta_imp, 'ro', label='solution') - ax.plot(t_exp, theta_exp, 'b-', label='simulated') - - ax.set_xlabel('time (s)') - ax.set_ylabel(r'$\theta$ (deg)') - ax.grid(True) - ax.legend(loc='upper right') - - plt.show() - @use_tempdirs class TestBrachistochronePolynomialControlExplicitRate2Targets(unittest.TestCase): diff --git a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py index 419601d86..5c3f2f48c 100644 --- a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py +++ b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py @@ -590,6 +590,31 @@ def test_brachistochrone_polynomial_control_gauss_lobatto(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() + io = p.model.get_io_metadata(metadata_keys=['tags']) + from pprint import pprint + pprint(io) + + for prob in [p, exp_out]: + ts_group = prob.model._get_subsystem('phase0.timeseries') + + map_type_to_promnames = {'dymos.type:time': {'time'}, + 'dymos.type:t_phase': set(), + 'dymos.type:control': set(), + 'dymos.type:polynomial_control': {'theta'}, + 'dymos.type:polynomial_control_rate': set(), + 'dymos.type:polynomial_control_rate2': set(), + 'dymos.type:parameter': set(), + 'dymos.type:state': {'x', 'y', 'v'}, + 'dymos.initial_boundary_constraint': set(), + 'dymos.final_boundary_constraint': set(), + 'dymos.path_constraint': {'theta_rate'}} + + for dymos_type, prom_names in map_type_to_promnames.items(): + prom_outputs = {meta['prom_name'] for abs_path, meta in + ts_group.list_outputs(tags=[dymos_type], out_stream=None)} + self.assertSetEqual(prom_outputs, prom_names, + msg=f'\n{dymos_type}\nin outputs: {prom_outputs}\nexpected: {prom_names}') + def test_brachistochrone_polynomial_control_radau(self): import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal @@ -645,11 +670,13 @@ def test_brachistochrone_polynomial_control_radau(self): # Generate the explicitly simulated trajectory exp_out = phase.simulate() + io = p.model.get_io_metadata(metadata_keys=['tags']) + from pprint import pprint + pprint(io) + for prob in [p, exp_out]: ts_group = prob.model._get_subsystem('phase0.timeseries') - ts_group.list_outputs() - map_type_to_promnames = {'dymos.type:time': {'time'}, 'dymos.type:t_phase': set(), 'dymos.type:control': set(), diff --git a/dymos/examples/brachistochrone/test/test_state_rate_introspection.py b/dymos/examples/brachistochrone/test/test_state_rate_introspection.py index 4c335feb1..8e9320a1c 100644 --- a/dymos/examples/brachistochrone/test/test_state_rate_introspection.py +++ b/dymos/examples/brachistochrone/test/test_state_rate_introspection.py @@ -722,10 +722,11 @@ def _test_transcription(self, transcription=dm.GaussLobatto, time_name='time'): # # Set the time options # Time has no targets in our ODE. - # We fix the initial time so that the it is not a design variable in the optimization. + # We fix the initial time so that it is not a design variable in the optimization. # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. # phase.set_time_options(fix_initial=True, fix_duration=True, units='s', name=time_name) + phase.timeseries_options['include_t_phase'] = True # # Set the time options diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index 6774710fb..152f26f3c 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -92,6 +92,7 @@ def __init__(self, from_phase=None, **kwargs): self.refine_options.update(from_phase.refine_options) self.simulate_options.update(from_phase.simulate_options) + self.timeseries_options.update(from_phase.timeseries_options) self._initial_boundary_constraints = from_phase._initial_boundary_constraints.copy() self._final_boundary_constraints = from_phase._final_boundary_constraints.copy() @@ -1182,7 +1183,6 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None, bc['flat_indices'] = flat_indices bc['is_expr'] = is_expr - print('adding timeseries output', name, constraint_name) self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape, tags=[f'dymos.{loc}_boundary_constraint']) diff --git a/dymos/phase/test/test_phase.py b/dymos/phase/test/test_phase.py index d079a3354..d67e612fa 100644 --- a/dymos/phase/test/test_phase.py +++ b/dymos/phase/test/test_phase.py @@ -405,23 +405,6 @@ def test_control_rate_boundary_constraint_gl(self): p.run_driver() - import matplotlib.pyplot as plt - - plt.plot(p.get_val('phase0.timeseries.x'), - p.get_val('phase0.timeseries.y'), 'ko') - - plt.figure() - - plt.plot(p.get_val('phase0.timeseries.time'), - p.get_val('phase0.timeseries.theta'), 'ro') - - plt.plot(p.get_val('phase0.timeseries.time'), - p.get_val('phase0.timeseries.theta_rate'), 'bo') - - plt.plot(p.get_val('phase0.timeseries.time'), - p.get_val('phase0.timeseries.theta_rate2'), 'go') - plt.show() - assert_near_equal(p.get_val('phase0.timeseries.theta_rate')[-1], 0, tolerance=1.0E-6) From c03adf59dc45b56a27d2ca1fb750d6c18aae4a18 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 06:42:55 -0400 Subject: [PATCH 09/20] fix for logic on when to include control rates in the timeseries outputs. --- dymos/transcriptions/transcription_base.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py index 68198a04e..a7c16acb9 100644 --- a/dymos/transcriptions/transcription_base.py +++ b/dymos/transcriptions/transcription_base.py @@ -88,8 +88,7 @@ def setup_time(self, phase): for ts_name, ts_options in phase._timeseries.items(): phase.add_timeseries_output(t_name, timeseries=ts_name, tags=['dymos.type:time']) - if (phase.timeseries_options['include_t_phase'] or - time_options['time_phase_targets'] is not _unspecified): + if phase.timeseries_options['include_t_phase']: phase.add_timeseries_output(t_phase_name, timeseries=ts_name, tags=['dymos.type:t_phase']) @@ -166,12 +165,10 @@ def setup_controls(self, phase): for ts_name, ts_options in phase._timeseries.items(): phase.add_timeseries_output(name, output_name=f'{control_prefix}{name}', timeseries=ts_name, tags=['dymos.type:control']) - if phase.timeseries_options['include_control_rates'] \ - or options['rate_targets'] is not _unspecified: + if phase.timeseries_options['include_control_rates']: phase.add_timeseries_output(f'{name}_rate', output_name=f'{control_rate_prefix}{name}_rate', timeseries=ts_name, tags=['dymos.type:control_rate']) - if phase.timeseries_options['include_control_rates'] \ - or options['rate2_targets'] is not _unspecified: + if phase.timeseries_options['include_control_rates']: phase.add_timeseries_output(f'{name}_rate2', output_name=f'{control_rate_prefix}{name}_rate2', timeseries=ts_name, tags=['dymos.type:control_rate2']) @@ -213,12 +210,10 @@ def setup_polynomial_controls(self, phase): for ts_name, ts_options in phase._timeseries.items(): phase.add_timeseries_output(name, output_name=f'{prefix}{name}', timeseries=ts_name, tags=['dymos.type:polynomial_control']) - if phase.timeseries_options['include_control_rates'] or \ - options['rate_targets'] is not _unspecified: + if phase.timeseries_options['include_control_rates']: phase.add_timeseries_output(f'{name}_rate', output_name=f'{rate_prefix}{name}_rate', timeseries=ts_name, tags=['dymos.type:polynomial_control_rate']) - if phase.timeseries_options['include_control_rates'] or \ - options['rate2_targets'] is not _unspecified: + if phase.timeseries_options['include_control_rates']: phase.add_timeseries_output(f'{name}_rate2', output_name=f'{rate_prefix}{name}_rate2', timeseries=ts_name, tags=['dymos.type:polynomial_control_rate2']) From 690b5a756b080e331d9413357eba7a48fdef7575 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 06:43:47 -0400 Subject: [PATCH 10/20] pep8 --- .../test/test_brachistochrone_control_rate_targets.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py index fabc56457..15cc7c7a4 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py @@ -455,6 +455,7 @@ def test_brachistochrone_control_rate_targets_radau(self): x_check=y_exp[estart:eend, :], rel_tolerance=1.0E-9) + @use_tempdirs class TestBrachistochronePolynomialControlRateTargets(unittest.TestCase): From c64fcc7243bcf51d82efaec2907172741db1c069 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 13:57:59 -0400 Subject: [PATCH 11/20] Added ability to toggle bounds to timeseries report --- .../doc/test_doc_balanced_field_length.py | 2 +- .../vanderpol/doc/test_doc_vanderpol.py | 10 +- dymos/examples/vanderpol/vanderpol_dymos.py | 7 +- .../timeseries/bokeh_timeseries_report.py | 104 ++++++++++-------- 4 files changed, 71 insertions(+), 52 deletions(-) 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 3b07d2f1c..b20ce591a 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 @@ -7,7 +7,7 @@ SHOW_PLOTS = True -@use_tempdirs +# @use_tempdirs class TestBalancedFieldLengthForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='IPOPT') diff --git a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py index 18d8d364e..82d5681cc 100644 --- a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py +++ b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py @@ -5,7 +5,7 @@ from openmdao.utils.mpi import MPI -@use_tempdirs +# @use_tempdirs class TestVanderpolForDocs(unittest.TestCase): def tearDown(self): for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']: @@ -45,10 +45,10 @@ 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) + # 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.") diff --git a/dymos/examples/vanderpol/vanderpol_dymos.py b/dymos/examples/vanderpol/vanderpol_dymos.py index 34a5e2572..c7a21b130 100644 --- a/dymos/examples/vanderpol/vanderpol_dymos.py +++ b/dymos/examples/vanderpol/vanderpol_dymos.py @@ -57,7 +57,7 @@ def vanderpol(transcription='gauss-lobatto', num_segments=40, transcription_orde units='V', targets='x1') # target required because x1 is an input phase.add_state('J', fix_initial=False, fix_final=False, - rate_source='Jdot', + rate_source='Jdot', lower=-10, units=None) # define the control @@ -69,9 +69,12 @@ def vanderpol(transcription='gauss-lobatto', num_segments=40, transcription_orde 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('x0', loc='final', lower=0, upper=0.1) phase.add_boundary_constraint('x1', loc='final', equals=0.0) + phase.add_path_constraint('x1', lower=-0.5, upper=1.5) + phase.add_path_constraint('u', lower=-1, upper=1.5) + # define objective to minimize phase.add_objective('J', loc='final') diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index 5d639e7d2..2b35aab64 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -5,8 +5,8 @@ try: from bokeh.io import output_notebook, output_file, save, show from bokeh.layouts import column, grid, row - from bokeh.models import Legend, DataTable, Div, ColumnDataSource, TableColumn, \ - TabPanel, Tabs, CheckboxButtonGroup, CustomJS, MultiChoice + from bokeh.models import Arrow, Legend, DataTable, Div, ColumnDataSource, TableColumn, \ + TabPanel, Tabs, CheckboxButtonGroup, CustomJS, MultiChoice, Whisker, VeeHead from bokeh.plotting import figure, curdoc import bokeh.palettes as bp import bokeh.resources as bokeh_resources @@ -20,6 +20,9 @@ import dymos as dm +BOUNDS_ALPHA = 0.08 +MARKER_SIZE = 8 + _js_show_renderer = """ function show_renderer(renderer, phases_to_show, kinds_to_show) { var tags = renderer.tags; @@ -36,29 +39,6 @@ """ -# Javascript Callback when the solution/simulation checkbox buttons are toggled -# args: (figures) -_SOL_SIM_TOGGLE_JS = """ -// Loop through figures and toggle the visibility of the renderers -const active = cb_obj.active; -var figures = figures; -var renderer; - -for (var i = 0; i < figures.length; i++) { - if (figures[i]) { - for (var j =0; j < figures[i].renderers.length; j++) { - renderer = figures[i].renderers[j] - if (renderer.tags.includes('sol')) { - renderer.visible = active.includes(0); - } - else if (renderer.tags.includes('sim')) { - renderer.visible = active.includes(1); - } - } - } -} -""" - # Javascript Callback when the solution/simulation checkbox buttons are toggled # args: (figures) _js_show_figures = """ @@ -77,8 +57,10 @@ } } return ((tags.includes('sol') && kinds_to_show.includes(0)) || - (tags.includes('sim') && kinds_to_show.includes(1))) && - phases_to_show.includes(renderer_phase); + (tags.includes('sim') && kinds_to_show.includes(1)) || + (tags.includes('bounds') && kinds_to_show.includes(2)) || + (tags.includes('constraints') && kinds_to_show.includes(3))) && + phases_to_show.includes(renderer_phase); } for (var i = 0; i < figures.length; i++) { @@ -149,7 +131,6 @@ def _load_data_sources(prob, solution_record_file=None, simulation_record_file=N data_dict[traj_name] = {'param_data_by_phase': {}, 'sol_data_by_phase': {}, 'sim_data_by_phase': {}, - 'bounds_by_phase': {}, 'timeseries_units': {}} for phase_name, phase in traj._phases.items(): @@ -278,7 +259,6 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi figures = [] x_range = None - size = 6 for var_name in sorted(ts_units_dict.keys(), key=str.casefold): fig_kwargs = {'x_range': x_range} if x_range is not None else {} @@ -298,9 +278,14 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi fig.yaxis.axis_label_text_font_size = '10pt' fig.toolbar.autohide = True legend_data = [] + renderers = [] if x_range is None: x_range = fig.x_range for i, phase_name in enumerate(phase_names): + phase = traj._phases[phase_name] + bcis = {con['constraint_name'] for con in phase._initial_boundary_constraints} + bcfs = {con['constraint_name'] for con in phase._final_boundary_constraints} + paths = {con['constraint_name'] for con in phase._path_constraints} color = colors[i % 20] sol_data = source_data[traj_name]['sol_data_by_phase'][phase_name] sim_data = source_data[traj_name]['sim_data_by_phase'][phase_name] @@ -310,6 +295,8 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi legend_items = [] if sol_data: lower = upper = None + bci_lower = bci_upper = bci_equals = bcf_lower = bcf_upper = \ + bcf_equals = path_lower = path_upper = path_equals = None fix_initial = False fix_final = False opt = True @@ -335,48 +322,76 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi upper = phase.parameter_options[var_name]['upper'] opt = phase.parameter_options[var_name]['opt'] + if var_name in bcis: + bci = [c for c in phase._initial_boundary_constraints + if c['constraint_name'] == var_name][0] + bci_lower = bci['equals'] or bci['lower'] + bci_upper = bci['equals'] or bci['upper'] + if var_name in bcfs: + bcf = [c for c in phase._final_boundary_constraints + if c['constraint_name'] == var_name][0] + bcf_lower = bcf['equals'] or bcf['lower'] + bcf_upper = bcf['equals'] or bcf['upper'] + if var_name in paths: + path = [c for c in phase._path_constraints + if c['constraint_name'] == var_name][0] + path_lower = path['equals'] or path['lower'] + path_upper = path['equals'] or path['upper'] + + x_data = sol_data[x_name].ravel() + if opt: sol_plot = fig.circle(x=x_name, y=var_name, source=sol_source, - color=color, size=size) + color=color, size=MARKER_SIZE) else: sol_plot = fig.circle_cross(x=x_data[0], y=sol_data[var_name][0, ...], color=color, fill_color='white', - size=size+2, line_width=1) + size=MARKER_SIZE+2, line_width=1) + renderers.append(sol_plot) sol_plot.tags.extend(['sol', f'phase:{phase_name}']) legend_items.append(sol_plot) # Plot the bounds if available - x_data = sol_data[x_name].ravel() - if lower: - lb_plot = fig.line(x=x_data, y=lower * np.ones_like(x_data), - line_dash="dashed") - lb_plot.tags.extend(['bounds']) - lb_plot.tags.extend(['bounds', f'phase:{phase_name}']) + if lower is not None: + sol_source.data[f'{var_name}:lower'] = lower * np.ones_like(x_data) + else: + sol_source.data[f'{var_name}:lower'] = -1.E8 * np.ones_like(x_data) - if upper: - ub_plot = fig.line(x=x_data, y=upper * np.ones_like(x_data), - line_dash="dashed") - ub_plot.tags.extend(['bounds', f'phase:{phase_name}']) + if upper is not None: + sol_source.data[f'{var_name}:upper'] = upper * np.ones_like(x_data) + else: + sol_source.data[f'{var_name}:upper'] = 1.E8 * np.ones_like(x_data) + if lower is not None or upper is not None: + band = fig.varea(x=x_name, y1=f'{var_name}:lower', y2=f'{var_name}:upper', + source=sol_source, fill_alpha=BOUNDS_ALPHA, + fill_color=color) + band.tags.extend([f'phase:{phase_name}', 'bounds']) + + # Plot fixed endpoints if appropriate if fix_initial: fix_initial_plot = fig.circle_cross(x=x_data[0], y=sol_data[var_name][0, ...], color=color, fill_color='white', - size=size+2, line_width=2) + size=MARKER_SIZE+2, line_width=2) fix_initial_plot.tags.extend(['sol', f'phase:{phase_name}']) if fix_final: fix_final_plot = fig.circle_cross(x=x_data[-1], y=sol_data[var_name][-1, ...], color=color, fill_color='white', - size=size+2, line_width=2) + size=MARKER_SIZE+2, line_width=2) fix_final_plot.tags.extend(['sol', f'phase:{phase_name}']) if sim_data: sim_plot = fig.line(x='time', y=var_name, source=sim_source, color=color) sim_plot.tags.extend(['sim', f'phase:{phase_name}']) legend_items.append(sim_plot) + renderers.append(sim_plot) legend_data.append((phase_name, legend_items)) + # Only scale the y axis based on the sol and sim data plots, not any bounds plots + fig.y_range.renderers = renderers + legend = Legend(items=legend_data, location='center', label_text_font_size='8pt') fig.add_layout(legend, 'right') figures.append(fig) @@ -388,7 +403,8 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi param_panels = [TabPanel(child=table, title=f'{phase_names[i]} parameters') for i, table in enumerate(param_tables)] - sol_sim_toggle = CheckboxButtonGroup(labels=['Solution', 'Simulation'], active=[0, 1]) + sol_sim_toggle = CheckboxButtonGroup(labels=['Solution', 'Simulation', 'Bounds'], + active=[0, 1, 2]) sol_sim_row = row(children=[Div(text='Display data:', sizing_mode='stretch_height'), sol_sim_toggle], From 36117c84cbea9378dea01142702739b581f0085d Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 14:16:48 -0400 Subject: [PATCH 12/20] reimplement use_tempdirs --- .../doc/test_doc_balanced_field_length.py | 2 +- .../test/test_double_integrator_timeseries_units.py | 2 +- .../min_time_climb/test/test_ex_min_time_climb.py | 2 +- dymos/examples/vanderpol/doc/test_doc_vanderpol.py | 2 +- .../pseudospectral/pseudospectral_base.py | 10 +++++----- 5 files changed, 9 insertions(+), 9 deletions(-) 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 b20ce591a..3b07d2f1c 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 @@ -7,7 +7,7 @@ SHOW_PLOTS = True -# @use_tempdirs +@use_tempdirs class TestBalancedFieldLengthForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='IPOPT') diff --git a/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py b/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py index 488a898f5..376991054 100644 --- a/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py +++ b/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py @@ -55,7 +55,7 @@ def double_integrator_direct_collocation(transcription=dm.GaussLobatto, compress return p -# @use_tempdirs +@use_tempdirs class TestDoubleIntegratorExample(unittest.TestCase): @classmethod 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 d7b77b58f..47a33cdfb 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 @@ -122,7 +122,7 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', return p -# @use_tempdirs +@use_tempdirs class TestMinTimeClimb(unittest.TestCase): def _test_results(self, p, time_name='time'): diff --git a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py index 82d5681cc..d06a61574 100644 --- a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py +++ b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py @@ -5,7 +5,7 @@ from openmdao.utils.mpi import MPI -# @use_tempdirs +@use_tempdirs class TestVanderpolForDocs(unittest.TestCase): def tearDown(self): for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']: diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py index cd744e9c2..5ee161106 100644 --- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py +++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py @@ -186,12 +186,12 @@ def configure_states(self, phase): f'initial_bounds for state {name} in phase {phase.name}') if options['input_initial']: raise ValueError('Cannot specify \'fix_initial=True\' and specify ' - f'\'connected_initial=True\' for state {name} ' + f'\'input_initial=True\' for state {name} ' f'in phase {phase.name}') idx_mask[0, ...] = np.asarray(np.logical_not(options['fix_initial']), dtype=int) elif options['input_initial']: if options['initial_bounds'] is not None: - raise ValueError('Cannot specify \'connected_initial=True\' and specify ' + raise ValueError('Cannot specify \'input_initial=True\' and specify ' f'initial_bounds for state {name} in phase {phase.name}') idx_mask[0, ...] = np.asarray(np.logical_not(options['input_initial']), dtype=int) @@ -596,10 +596,10 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None): shape = phase.state_options[var]['shape'] units = phase.state_options[var]['units'] solve_segments = phase.state_options[var]['solve_segments'] - connected_initial = phase.state_options[var]['input_initial'] - if not solve_segments and not connected_initial: + input_initial = phase.state_options[var]['input_initial'] + if not solve_segments and not input_initial: linear = True - elif solve_segments in {'forward'} and not connected_initial and loc == 'initial': + elif solve_segments in {'forward'} and not input_initial and loc == 'initial': linear = True elif solve_segments == 'backward' and loc == 'final': linear = True From c077ba80e2c395a0f65ccf8ca9fb4ec3a022e203 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 15:18:43 -0400 Subject: [PATCH 13/20] vanderpol problem reset --- dymos/examples/vanderpol/vanderpol_dymos.py | 2 +- dymos/trajectory/trajectory.py | 13 +++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/dymos/examples/vanderpol/vanderpol_dymos.py b/dymos/examples/vanderpol/vanderpol_dymos.py index c7a21b130..bd7dfe129 100644 --- a/dymos/examples/vanderpol/vanderpol_dymos.py +++ b/dymos/examples/vanderpol/vanderpol_dymos.py @@ -69,7 +69,7 @@ def vanderpol(transcription='gauss-lobatto', num_segments=40, transcription_orde 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', lower=0, upper=0.1) + phase.add_boundary_constraint('x0', loc='final', equals=0.0) phase.add_boundary_constraint('x1', loc='final', equals=0.0) phase.add_path_constraint('x1', lower=-0.5, upper=1.5) diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py index ffab7c13f..3a73ce655 100644 --- a/dymos/trajectory/trajectory.py +++ b/dymos/trajectory/trajectory.py @@ -51,7 +51,16 @@ class Trajectory(om.Group): _phase_graph : nx.DiGraph A graph of linked phases. """ - def __init__(self, **kwargs): + def __init__(self, parallel_phases=True, **kwargs): + """ + Initialize a Trajectory instance. + + Parameters + ---------- + parallel_phases : bool + If True, make the phases subsystem a ParallelGroup otherwise + it will be a standard OpenMDAO Group. + """ super(Trajectory, self).__init__(**kwargs) self._linkages = {} @@ -59,7 +68,7 @@ def __init__(self, **kwargs): self._phase_graph = nx.DiGraph() self.parameter_options = {} - self.phases = om.ParallelGroup() + self.phases = om.ParallelGroup() if parallel_phases else om.Group() def initialize(self): """ From 52b36a7030796006471ca0c2a5c467386c46d465 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 15:43:59 -0400 Subject: [PATCH 14/20] cleanup --- dymos/trajectory/trajectory.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py index 3a73ce655..ffab7c13f 100644 --- a/dymos/trajectory/trajectory.py +++ b/dymos/trajectory/trajectory.py @@ -51,16 +51,7 @@ class Trajectory(om.Group): _phase_graph : nx.DiGraph A graph of linked phases. """ - def __init__(self, parallel_phases=True, **kwargs): - """ - Initialize a Trajectory instance. - - Parameters - ---------- - parallel_phases : bool - If True, make the phases subsystem a ParallelGroup otherwise - it will be a standard OpenMDAO Group. - """ + def __init__(self, **kwargs): super(Trajectory, self).__init__(**kwargs) self._linkages = {} @@ -68,7 +59,7 @@ def __init__(self, parallel_phases=True, **kwargs): self._phase_graph = nx.DiGraph() self.parameter_options = {} - self.phases = om.ParallelGroup() if parallel_phases else om.Group() + self.phases = om.ParallelGroup() def initialize(self): """ From 50ed0ba2c8608077a7542b77df174b49efc14d02 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Tue, 8 Aug 2023 16:42:40 -0400 Subject: [PATCH 15/20] doc cleanup --- .../ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb b/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb index 70dc7310e..3884e3a62 100644 --- a/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb +++ b/docs/dymos_book/examples/ssto_moon_linear_tangent/ssto_moon_linear_tangent.ipynb @@ -305,6 +305,8 @@ "phase = dm.Phase(ode_class=LaunchVehicleLinearTangentODE,\n", " transcription=dm.GaussLobatto(num_segments=10, order=5, compressed=True))\n", "\n", + "phase.timeseries_options['include_t_phase'] = True\n", + "\n", "traj.add_phase('phase0', phase)\n", "\n", "phase.set_time_options(fix_initial=True, duration_bounds=(10, 1000),\n", @@ -482,7 +484,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.11.4" } }, "nbformat": 4, From a37c8be31a8b3f8ea1ca82b9077886d7ad3ae88c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 9 Aug 2023 08:52:18 -0400 Subject: [PATCH 16/20] Fixed issues with timeseries plots when the name of time in a phase is not "time". --- dymos/visualization/timeseries/bokeh_timeseries_report.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index 2b35aab64..f7b8e1705 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -192,7 +192,7 @@ def _load_data_sources(prob, solution_record_file=None, simulation_record_file=N def make_timeseries_report(prob, solution_record_file=None, simulation_record_file=None, - x_name='time', ncols=2, margin=10, theme='light_minimal'): + ncols=2, margin=10, theme='light_minimal'): """ Create the bokeh-based timeseries results report. @@ -260,6 +260,9 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi figures = [] x_range = None + # The x_axis label will use the name for time in the first phase. + x_name = list(traj._phases.values())[0].time_options['name'] + for var_name in sorted(ts_units_dict.keys(), key=str.casefold): fig_kwargs = {'x_range': x_range} if x_range is not None else {} From e16fd24e007c63b0a9dc0c845f4682256bad6435 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 9 Aug 2023 10:12:37 -0400 Subject: [PATCH 17/20] Cleaned up representation of bound and path constraints, still need bci/bcf --- .../test/test_ex_min_time_climb.py | 2 +- .../timeseries/bokeh_timeseries_report.py | 51 +++++++++++++++---- 2 files changed, 42 insertions(+), 11 deletions(-) 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 47a33cdfb..d7b77b58f 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 @@ -122,7 +122,7 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', return p -@use_tempdirs +# @use_tempdirs class TestMinTimeClimb(unittest.TestCase): def _test_results(self, p, time_name='time'): diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index f7b8e1705..a29d5104f 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -5,8 +5,8 @@ try: from bokeh.io import output_notebook, output_file, save, show from bokeh.layouts import column, grid, row - from bokeh.models import Arrow, Legend, DataTable, Div, ColumnDataSource, TableColumn, \ - TabPanel, Tabs, CheckboxButtonGroup, CustomJS, MultiChoice, Whisker, VeeHead + from bokeh.models import BoxAnnotation, Legend, DataTable, Div, ColumnDataSource, TableColumn, \ + TabPanel, Tabs, CheckboxButtonGroup, CustomJS, MultiChoice from bokeh.plotting import figure, curdoc import bokeh.palettes as bp import bokeh.resources as bokeh_resources @@ -20,8 +20,15 @@ import dymos as dm -BOUNDS_ALPHA = 0.08 +BOUNDS_ALPHA = 0.1 +BOUNDS_HATCH_PATTERN = '/' +BOUNDS_HATCH_ALPHA = 0.5 +PATH_ALPHA = 0.1 +PATH_HATCH_PATTERN = 'x' +PATH_HATCH_ALPHA = 0.5 MARKER_SIZE = 8 +MIN_Y = -1.0E8 +MAX_Y = 1.0E8 _js_show_renderer = """ function show_renderer(renderer, phases_to_show, kinds_to_show) { @@ -63,11 +70,13 @@ phases_to_show.includes(renderer_phase); } +console.log(figures[0]) for (var i = 0; i < figures.length; i++) { if (figures[i]) { for (var j=0; j < figures[i].renderers.length; j++) { renderer = figures[i].renderers[j]; // Get the phase with which this renderer is associated + console.log(renderer.tags) renderer.visible = show_renderer(renderer, phases_to_show, kinds_to_show); } } @@ -286,6 +295,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi x_range = fig.x_range for i, phase_name in enumerate(phase_names): phase = traj._phases[phase_name] + x_name = phase.time_options['name'] bcis = {con['constraint_name'] for con in phase._initial_boundary_constraints} bcfs = {con['constraint_name'] for con in phase._final_boundary_constraints} paths = {con['constraint_name'] for con in phase._path_constraints} @@ -340,6 +350,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi if c['constraint_name'] == var_name][0] path_lower = path['equals'] or path['lower'] path_upper = path['equals'] or path['upper'] + print(var_name, path_lower, path_upper, path['equals']) x_data = sol_data[x_name].ravel() @@ -366,11 +377,31 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi else: sol_source.data[f'{var_name}:upper'] = 1.E8 * np.ones_like(x_data) - if lower is not None or upper is not None: - band = fig.varea(x=x_name, y1=f'{var_name}:lower', y2=f'{var_name}:upper', - source=sol_source, fill_alpha=BOUNDS_ALPHA, - fill_color=color) - band.tags.extend([f'phase:{phase_name}', 'bounds']) + if lower is not None: + fig.varea(x=[x_data[0], x_data[-1]], y1=[lower, lower], y2=[MIN_Y, MIN_Y], + fill_alpha=BOUNDS_ALPHA, fill_color=color, + hatch_pattern=BOUNDS_HATCH_PATTERN, hatch_color=color, + hatch_alpha=BOUNDS_HATCH_ALPHA, + hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'bounds']) + + if upper is not None: + fig.varea(x=[x_data[0], x_data[-1]], y1=[upper, upper], y2=[MAX_Y, MAX_Y], + fill_alpha=BOUNDS_ALPHA, fill_color=color, + hatch_pattern=BOUNDS_HATCH_PATTERN, hatch_color=color, + hatch_alpha=BOUNDS_HATCH_ALPHA, + hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'bounds']) + + if path_lower is not None: + fig.varea(x=[x_data[0], x_data[-1]], y1=[path_lower, path_lower], y2=[MIN_Y, MIN_Y], + fill_alpha=PATH_ALPHA, fill_color=color, + hatch_pattern=PATH_HATCH_PATTERN, hatch_color=color, hatch_alpha=PATH_HATCH_ALPHA, + hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'constraints']) + + if path_upper is not None: + fig.varea(x=[x_data[0], x_data[-1]], y1=[path_upper, path_upper], y2=[MAX_Y, MAX_Y], + fill_alpha=PATH_ALPHA, fill_color=color, + hatch_pattern=PATH_HATCH_PATTERN, hatch_color=color, hatch_alpha=PATH_HATCH_ALPHA, + hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'constraints']) # Plot fixed endpoints if appropriate if fix_initial: @@ -386,7 +417,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi fix_final_plot.tags.extend(['sol', f'phase:{phase_name}']) if sim_data: - sim_plot = fig.line(x='time', y=var_name, source=sim_source, color=color) + sim_plot = fig.line(x=x_name, y=var_name, source=sim_source, color=color) sim_plot.tags.extend(['sim', f'phase:{phase_name}']) legend_items.append(sim_plot) renderers.append(sim_plot) @@ -406,7 +437,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi param_panels = [TabPanel(child=table, title=f'{phase_names[i]} parameters') for i, table in enumerate(param_tables)] - sol_sim_toggle = CheckboxButtonGroup(labels=['Solution', 'Simulation', 'Bounds'], + sol_sim_toggle = CheckboxButtonGroup(labels=['Solution', 'Simulation', 'Bounds', 'Constraints'], active=[0, 1, 2]) sol_sim_row = row(children=[Div(text='Display data:', sizing_mode='stretch_height'), From 1fc0103976fa0ee8d4f93b32fae92e57a2c01c5c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 9 Aug 2023 12:56:09 -0400 Subject: [PATCH 18/20] Added path constraints and some cleanup --- .../test/test_ex_min_time_climb.py | 2 +- dymos/test/test_run_problem.py | 1 + .../test/test_timeseries_plots.py | 5 +- .../timeseries/bokeh_timeseries_report.py | 62 ++++++++++--------- 4 files changed, 38 insertions(+), 32 deletions(-) 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 d7b77b58f..47a33cdfb 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 @@ -122,7 +122,7 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', return p -# @use_tempdirs +@use_tempdirs class TestMinTimeClimb(unittest.TestCase): def _test_results(self, p, time_name='time'): diff --git a/dymos/test/test_run_problem.py b/dymos/test/test_run_problem.py index 95fd0a950..f351450cf 100755 --- a/dymos/test/test_run_problem.py +++ b/dymos/test/test_run_problem.py @@ -660,6 +660,7 @@ def setUp(self): phase0.set_refine_options(refine=True) phase0.timeseries_options['include_state_rates'] = True + phase0.timeseries_options['include_control_rates'] = True phase0.timeseries_options['include_t_phase'] = True p.model.linear_solver = om.DirectSolver() diff --git a/dymos/visualization/test/test_timeseries_plots.py b/dymos/visualization/test/test_timeseries_plots.py index fdd997c93..6145cd833 100644 --- a/dymos/visualization/test/test_timeseries_plots.py +++ b/dymos/visualization/test/test_timeseries_plots.py @@ -58,6 +58,8 @@ def setUp(self): compressed=compressed), subset='control_input') + phase.timeseries_options['include_control_rates'] = 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 @@ -83,7 +85,7 @@ def test_brachistochrone_timeseries_plots(self): temp = dm.options['plots'] dm.options['plots'] = 'matplotlib' - dm.run_problem(self.p, make_plots=False) + dm.run_problem(self.p, make_plots=True) timeseries_plots('dymos_solution.db', problem=self.p) plot_dir = pathlib.Path(_get_reports_dir(self.p)).joinpath('plots') @@ -256,6 +258,7 @@ def test_trajectory_linked_phases_make_plot(self): for phase_name, phase in self.traj._phases.items(): phase.timeseries_options['use_prefix'] = True phase.timeseries_options['include_state_rates'] = True + phase.timeseries_options['include_control_rates'] = True phase.timeseries_options['include_t_phase'] = True p.setup(check=True) diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index a29d5104f..7d060813c 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -27,8 +27,8 @@ PATH_HATCH_PATTERN = 'x' PATH_HATCH_ALPHA = 0.5 MARKER_SIZE = 8 -MIN_Y = -1.0E8 -MAX_Y = 1.0E8 +MIN_Y = -1.0E6 +MAX_Y = 1.0E6 _js_show_renderer = """ function show_renderer(renderer, phases_to_show, kinds_to_show) { @@ -213,8 +213,6 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi The path to the solution record file, if available. simulation_record_file : str The path to the simulation record file, if available. - x_name : str - Name of the horizontal axis variable in the timeseries. ncols : int The number of columns of timeseries output plots. margin : int @@ -338,19 +336,18 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi if var_name in bcis: bci = [c for c in phase._initial_boundary_constraints if c['constraint_name'] == var_name][0] - bci_lower = bci['equals'] or bci['lower'] - bci_upper = bci['equals'] or bci['upper'] + bci_lower = bci['equals'] if bci['equals'] is not None else bci['lower'] + bci_upper = bci['equals'] if bci['equals'] is not None else bci['upper'] if var_name in bcfs: bcf = [c for c in phase._final_boundary_constraints if c['constraint_name'] == var_name][0] - bcf_lower = bcf['equals'] or bcf['lower'] - bcf_upper = bcf['equals'] or bcf['upper'] + bcf_lower = bcf['equals'] if bcf['equals'] is not None else bcf['lower'] + bcf_upper = bcf['equals'] if bcf['equals'] is not None else bcf['upper'] if var_name in paths: path = [c for c in phase._path_constraints if c['constraint_name'] == var_name][0] - path_lower = path['equals'] or path['lower'] - path_upper = path['equals'] or path['upper'] - print(var_name, path_lower, path_upper, path['equals']) + path_lower = path['equals'] if path['equals'] is not None else path['lower'] + path_upper = path['equals'] if path['equals'] is not None else path['upper'] x_data = sol_data[x_name].ravel() @@ -378,30 +375,34 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi sol_source.data[f'{var_name}:upper'] = 1.E8 * np.ones_like(x_data) if lower is not None: - fig.varea(x=[x_data[0], x_data[-1]], y1=[lower, lower], y2=[MIN_Y, MIN_Y], - fill_alpha=BOUNDS_ALPHA, fill_color=color, - hatch_pattern=BOUNDS_HATCH_PATTERN, hatch_color=color, - hatch_alpha=BOUNDS_HATCH_ALPHA, - hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'bounds']) + area = fig.varea(x=[x_data[0], x_data[-1]], y1=[lower, lower], y2=[MIN_Y, MIN_Y], + fill_alpha=BOUNDS_ALPHA, fill_color=color, + hatch_pattern=BOUNDS_HATCH_PATTERN, hatch_color=color, + hatch_alpha=BOUNDS_HATCH_ALPHA, + hatch_weight=0.2) + area.tags.extend([f'phase:{phase_name}', 'bounds'],) if upper is not None: - fig.varea(x=[x_data[0], x_data[-1]], y1=[upper, upper], y2=[MAX_Y, MAX_Y], - fill_alpha=BOUNDS_ALPHA, fill_color=color, - hatch_pattern=BOUNDS_HATCH_PATTERN, hatch_color=color, - hatch_alpha=BOUNDS_HATCH_ALPHA, - hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'bounds']) + area = fig.varea(x=[x_data[0], x_data[-1]], y1=[upper, upper], y2=[MAX_Y, MAX_Y], + fill_alpha=BOUNDS_ALPHA, fill_color=color, + hatch_pattern=BOUNDS_HATCH_PATTERN, hatch_color=color, + hatch_alpha=BOUNDS_HATCH_ALPHA, + hatch_weight=0.2) + area.tags.extend([f'phase:{phase_name}', 'bounds']) if path_lower is not None: - fig.varea(x=[x_data[0], x_data[-1]], y1=[path_lower, path_lower], y2=[MIN_Y, MIN_Y], - fill_alpha=PATH_ALPHA, fill_color=color, - hatch_pattern=PATH_HATCH_PATTERN, hatch_color=color, hatch_alpha=PATH_HATCH_ALPHA, - hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'constraints']) + area = fig.varea(x=[x_data[0], x_data[-1]], y1=[path_lower, path_lower], y2=[MIN_Y, MIN_Y], + fill_alpha=PATH_ALPHA, fill_color=color, + hatch_pattern=PATH_HATCH_PATTERN, hatch_color=color, + hatch_alpha=PATH_HATCH_ALPHA, hatch_weight=0.2) + area.tags.extend([f'phase:{phase_name}', 'constraints']) if path_upper is not None: - fig.varea(x=[x_data[0], x_data[-1]], y1=[path_upper, path_upper], y2=[MAX_Y, MAX_Y], - fill_alpha=PATH_ALPHA, fill_color=color, - hatch_pattern=PATH_HATCH_PATTERN, hatch_color=color, hatch_alpha=PATH_HATCH_ALPHA, - hatch_weight=0.2).tags.extend([f'phase:{phase_name}', 'constraints']) + area = fig.varea(x=[x_data[0], x_data[-1]], y1=[path_upper, path_upper], y2=[MAX_Y, MAX_Y], + fill_alpha=PATH_ALPHA, fill_color=color, + hatch_pattern=PATH_HATCH_PATTERN, hatch_color=color, + hatch_alpha=PATH_HATCH_ALPHA, hatch_weight=0.2) + area.tags.extend([f'phase:{phase_name}', 'constraints']) # Plot fixed endpoints if appropriate if fix_initial: @@ -426,6 +427,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi # Only scale the y axis based on the sol and sim data plots, not any bounds plots fig.y_range.renderers = renderers + # Add a phase legend outside of the axes legend = Legend(items=legend_data, location='center', label_text_font_size='8pt') fig.add_layout(legend, 'right') figures.append(fig) @@ -438,7 +440,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi for i, table in enumerate(param_tables)] sol_sim_toggle = CheckboxButtonGroup(labels=['Solution', 'Simulation', 'Bounds', 'Constraints'], - active=[0, 1, 2]) + active=[0, 1, 2, 3]) sol_sim_row = row(children=[Div(text='Display data:', sizing_mode='stretch_height'), sol_sim_toggle], From 60a3912c12271eede0aa2b354dce66d593ed31e6 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Thu, 10 Aug 2023 10:01:58 -0400 Subject: [PATCH 19/20] Fix for older versions of OpenMDAO. --- dymos/phase/test/test_timeseries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py index 541a6b06e..7904b61c3 100644 --- a/dymos/phase/test/test_timeseries.py +++ b/dymos/phase/test/test_timeseries.py @@ -118,7 +118,7 @@ def test_timeseries_gl(self, test_smaller_timeseries=False): for dymos_type, prom_names in map_type_to_promnames.items(): prom_outputs = {meta['prom_name'] for abs_path, meta in - ts_group.list_outputs(tags=[dymos_type], out_stream=None)} + ts_group.list_outputs(tags=[dymos_type], out_stream=None, prom_name=True)} self.assertSetEqual(prom_outputs, prom_names) def test_timeseries_gl_smaller_timeseries(self): From 04086c5254c3eff8e01c677ca9fc133885cfa29c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Thu, 10 Aug 2023 11:10:18 -0400 Subject: [PATCH 20/20] more fixes for compatibility with older openmdao versions. --- .../test/test_doc_brachistochrone_polynomial_controls.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py index 5c3f2f48c..11fd4881e 100644 --- a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py +++ b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py @@ -611,7 +611,7 @@ def test_brachistochrone_polynomial_control_gauss_lobatto(self): for dymos_type, prom_names in map_type_to_promnames.items(): prom_outputs = {meta['prom_name'] for abs_path, meta in - ts_group.list_outputs(tags=[dymos_type], out_stream=None)} + ts_group.list_outputs(tags=[dymos_type], prom_name=True, out_stream=None)} self.assertSetEqual(prom_outputs, prom_names, msg=f'\n{dymos_type}\nin outputs: {prom_outputs}\nexpected: {prom_names}') @@ -691,7 +691,7 @@ def test_brachistochrone_polynomial_control_radau(self): for dymos_type, prom_names in map_type_to_promnames.items(): prom_outputs = {meta['prom_name'] for abs_path, meta in - ts_group.list_outputs(tags=[dymos_type], out_stream=None)} + ts_group.list_outputs(tags=[dymos_type], prom_name=True, out_stream=None)} self.assertSetEqual(prom_outputs, prom_names, msg=f'\n{dymos_type}\nin outputs: {prom_outputs}\nexpected: {prom_names}')