Calculation Expressions#
Sometimes a user may wish to calculate some additional output that is not inherently provided by the ODE system. In cases where they have access to the model, it might be trivial to add it to the ODE’s outputs. In other cases, they might not want to modify a model provided by someone else.
Calculation expressions provide the user with a way of computing additional outputs that is extremely general. Older version of dymos supported expressions for use as constraints, objectives, or timeseries outputs. Now dymos just computes these expressions along-side the ODE, making them available to basically all of dymos (constraints, objectives, timeseries outputs, state rates, etc.)
The most general way to add a calculated expression to an ODE output is to use the phase.add_calc_expr method.
- Phase.add_calc_expr(expr, add_timeseries=True, **kwargs)[source]
Adds an expression to be computed immediately after the user-given ODE.
Internally, dymos will wrap the user-given ODE in a group along with an ExecComp that evalutes the expressions given.
Unlike standard OpenMDAO ExecComp expressions, those specified here may use the ODE-relative path of variables (if they are not promoted to the top of the ODE).
This function may be called more than one time with the same expression. In that case, any new info from kwargs will be updated.
- Parameters:
- exprstr
The expression to be computed.
- add_timeseriesbool
If True, add the output of the expression to the timeseries.
- **kwargsdict
Any arguments to be forwarded to the ExecComp when the expression in added.
Providing the calculated expression as the name in add_boundary_constraint, add_path_constraint, add_objective, or add_timeseries_output will still work for backward compatibility.
In general, using add_calc_expr will just always work, and makes that output variable available across the phase.
Units#
In general, a given calculation that a user might apply assumes specific units for each of the input variables. In dymos, introspection for units happens after the ODE is completely defined. Because the ExecComp that computes the calcuation expressions is itself part of the ODE, we cannot wait for this introspection step. Units for inputs will default to None if not otherwise specified.
If you’re not certain of the units of an input variable upon which your calculation relies, you should explicitly provide it via keyword arguments to the add_calc_expr call (or one of the associated calls that accepts calculation expressions, mentioned above).
Example - Adding an additional calculation to the brachistochrone#
Let’s assume that we have a brachistochrone ODE that was provided to us, but for configuration control reasons, we might not want to edit the source code.
Here we’re adding a calculation expression for the additional output k, which is the ratio
In the brachistochrone solution, \(k\) is approximately constant, and its value is dependent upon the start and end point of the trajectory, as well as the gravitational parameter. In practice, it’s value is somewhat large at the start since we’re starting with nearly zero velocity.
import openmdao.api as om
import dymos as dm
from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
p = om.Problem(model=om.Group())
p.driver = om.ScipyOptimizeDriver()
p.driver.declare_coloring()
# Create the transcription so we can get the number of nodes for the downstream analysis
tx = dm.Radau(num_segments=20, order=3, compressed=False)
traj = dm.Trajectory()
phase = dm.Phase(transcription=tx, ode_class=BrachistochroneODE)
traj.add_phase('phase0', phase)
p.model.add_subsystem('traj', traj)
phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
phase.add_state('x', units='m', rate_source='xdot', fix_initial=True, fix_final=True)
phase.add_state('y', units='m', rate_source='ydot', fix_initial=True, fix_final=True)
phase.add_state('v', units='m/s', rate_source='vdot', fix_initial=True, fix_final=False)
phase.add_control('theta', units='deg', lower=0.01, upper=179.9,
continuity=True, rate_continuity=True)
phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
phase.add_calc_expr('k = sin(theta) / v', theta={'units': 'rad'}, v={'units': 'm/s'})
# Minimize time at the end of the phase
phase.add_objective('time', loc='final', scaler=1)
p.setup(check=True)
phase.set_time_val(initial=0.0, duration=2.0)
phase.set_state_val('x', [0, 10])
phase.set_state_val('y', [10, 5])
phase.set_state_val('v', [1.0E-1, 9.9])
phase.set_control_val('theta', [5, 100])
phase.set_parameter_val('g', 9.80665)
dm.run_problem(p, simulate=True, make_plots=True)
INFO: checking out_of_order...
INFO: out_of_order check complete (0.000314 sec).
INFO: checking system...
INFO: system check complete (0.000022 sec).
INFO: checking solvers...
INFO: solvers check complete (0.000261 sec).
INFO: checking dup_inputs...
INFO: dup_inputs check complete (0.000078 sec).
INFO: checking missing_recorders...
INFO: missing_recorders check complete (0.000002 sec).
INFO: checking unserializable_options...
INFO: unserializable_options check complete (0.001265 sec).
INFO: checking comp_has_no_outputs...
INFO: comp_has_no_outputs check complete (0.000047 sec).
INFO: checking auto_ivc_warnings...
INFO: auto_ivc_warnings check complete (0.000003 sec).
INFO: checking out_of_order...
INFO: out_of_order check complete (0.000236 sec).
INFO: checking system...
INFO: system check complete (0.000013 sec).
INFO: checking solvers...
INFO: solvers check complete (0.000127 sec).
INFO: checking dup_inputs...
INFO: dup_inputs check complete (0.000043 sec).
INFO: checking missing_recorders...
INFO: missing_recorders check complete (0.000002 sec).
INFO: checking unserializable_options...
INFO: unserializable_options check complete (0.000985 sec).
INFO: checking comp_has_no_outputs...
INFO: comp_has_no_outputs check complete (0.000032 sec).
INFO: checking auto_ivc_warnings...
INFO: auto_ivc_warnings check complete (0.000002 sec).
Jacobian shape: (219, 296) (2.15% nonzero)
FWD solves: 11 REV solves: 0
Total colors vs. total size: 11 vs 296 (96.28% improvement)
Sparsity computed using tolerance: 1e-25.
Dense total jacobian for Problem 'problem' was computed 3 times.
Time to compute sparsity: 0.3754 sec
Time to compute coloring: 0.2817 sec
Memory to compute coloring: 0.2500 MB
Coloring created on: 2025-03-13 20:13:07
Optimization terminated successfully (Exit mode 0)
Current function value: 1.7931378025881326
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
Optimization Complete
-----------------------------------
/usr/share/miniconda/envs/test/lib/python3.11/site-packages/openmdao/visualization/opt_report/opt_report.py:619: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
ax.set_ylim([ymin_plot, ymax_plot])
/usr/share/miniconda/envs/test/lib/python3.11/site-packages/openmdao/core/group.py:1166: DerivativesWarning:Constraints or objectives [ode_eval.control_interp.control_rates:theta_rate, ode_eval.control_interp.control_rates:theta_rate2, ode_eval.control_interp.control_values:theta] cannot be impacted by the design variables of the problem because no partials were defined for them in their parent component(s).
Simulating trajectory traj
Done simulating trajectory traj
Problem: problem
Driver: ScipyOptimizeDriver
success : True
iterations : 6
runtime : 3.2336E+00 s
model_evals : 6
model_time : 4.5761E-03 s
deriv_evals : 5
deriv_time : 4.2186E-02 s
exit_status : SUCCESS
Show code cell source
from IPython.display import HTML
# Define the path to the HTML file
html_file_path = p.get_reports_dir() / 'traj_results_report.html'
html_content = html_file_path.read_text()
# Inject CSS to control the output cell height and avoid scrollbars
html_with_custom_height = f"""
<div style="height: 800px; overflow: auto;">
{html_content}
</div>
"""
HTML(html_with_custom_height)