(1) is the [-1] index of the timeseries variables always the last time point in the transcription?
Yes, it is.
(2) Can I add one or more specific time points to the nodes? Is it possible to have nodes guaranteed to be at specific fixed time points that might not match up with collocation points?
Yes...
There is a "right" way to do this which you should most certainly try first. Then there is a "wrong" way which I'll mention for completeness, and might work for you but also likely has some numerical challenges that I won't go into detail on but are the reason its not the right way.
--- The Right Way ---
If you have some intermediate times where you need to measure/constrain things, then you should break your trajectory up into multiple phases, and you can set the breakpoints to be at the key times you need to consider. So if you have a trajectory where the duration is uncertain but will be at least 10 seconds long, and you know you will want to do some post-processing on the performance at t=10, then make it have two phases. One phase from t=0 to t=10, and a second phase from t=10 to t=t_final. Even if t_final happens to equal 10, you're still ok. Dymos allows for 0 duration phases.
--- The other way ---
If you set up your problem with a fixed duration, then this is definitely possible. Its possible that you get lucky and collocation points line up with the times you want. If you don't get lucky though, you can use another technique to get Dymos to interpolate the time history onto a secondary grid for you. That grid can have a different transcription and number of segments than the primary collocation one (but will have the same start and end times). You can choose a Radau transcription with an order of 1, which gives you a simple fixed time-step. Then the number of segments basically equals the number of time-steps, and you could then in theory ensure you have an interpolated value at any time that you can reach with an integer division of your total duration. If you needed two times, that didn't line up on the same integer division then you could add two additional timeseries with different numbers of segments.
I mention this method, only because there are some problems where you don't just want to take a measurement at one time, but you want it at some evenly spaced distribution of times (maybe you want to take a FFT of the data...). In these cases, the secondary time-series is a good option.
If your desired time doesn't line up on an integer division, then definitely don't use this approach. Also, if your phase has a non-fixed duration, then REALLY definitely don't use this approach. It simply won't work, since you can't hold the time points fixed when the duration is changing.
(3.1) Can you have multiple objectives?
No. Dymos uses gradient based optimization, and optimal control problems are not well posed for generic multi-objective problems. You need to compose a composite objective via some sort of explicit combinations of outputs.
One objective only. Call EITHER the dymos objective api, or the OpenMDAO one. Not both.
(3.2), How can you set up postprocessing with an integrated value?
Dymos is fundamentally an integrator. If you need a value integrated, you simply add it as a new state to the phase. You don't need to edit the ODE at all. Dymos understands what you are trying to do if you pass it an existing state as the rate source for a new state.
So for the brachistochrone problem if you wanted to track total integrated y motion you would add a new state (we'll call it y_tot, with the state_rate being y itself). Then the last time point in the timeseries for y_tot will have the integrated value.
Here is a script that shows the techniques I described for (2) and (3). Adding more phases with breakpoints at specific times is covered in this Dymos docs example.
import numpy as np
import openmdao.api as om
import dymos as dm
import matplotlib.pyplot as plt
# First define a system which computes the equations of motion
class BrachistochroneEOM(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_nodes', types=int)
def setup(self):
nn = self.options['num_nodes']
# Inputs
self.add_input('v', val=np.zeros(nn), units='m/s', desc='velocity')
self.add_input('theta', val=np.zeros(nn), units='rad', desc='angle of wire')
self.add_output('xdot', val=np.zeros(nn), units='m/s', desc='x rate of change')
self.add_output('ydot', val=np.zeros(nn), units='m/s', desc='y rate of change')
self.add_output('vdot', val=np.zeros(nn), units='m/s**2', desc='v rate of change')
# Ask OpenMDAO to compute the partial derivatives using complex-step
# with a partial coloring algorithm for improved performance
self.declare_partials(of='*', wrt='*', method='cs')
self.declare_coloring(wrt='*', method='cs', show_summary=True)
def compute(self, inputs, outputs):
v, theta = inputs.values()
outputs['vdot'] = 9.80665 * np.cos(theta)
outputs['xdot'] = v * np.sin(theta)
outputs['ydot'] = -v * np.cos(theta)
p = om.Problem()
# Define a Trajectory object
traj = p.model.add_subsystem('traj', dm.Trajectory())
# Define a Dymos Phase object with GaussLobatto Transcription
tx = dm.GaussLobatto(num_segments=10, order=3)
phase = dm.Phase(ode_class=BrachistochroneEOM, transcription=tx)
traj.add_phase(name='phase0', phase=phase)
# Set the time options
phase.set_time_options(fix_initial=True,
duration_bounds=(0.5, 10.0))
# Set the state options
phase.add_state('x', rate_source='xdot',
fix_initial=True, fix_final=True)
phase.add_state('y', rate_source='ydot',
fix_initial=True, fix_final=True)
phase.add_state('y_tot', rate_source='y',
fix_initial=True, fix_final=True)
phase.add_state('v', rate_source='vdot',
fix_initial=True, fix_final=False)
# Define theta as a control.
phase.add_control(name='theta', units='rad',
lower=0, upper=np.pi)
# compressed transcription so we don't get duplicate nodes between each segment
phase.add_timeseries('fixed_step', dm.Radau(num_segments=10, order=1))
phase.add_timeseries_output('v', timeseries="fixed_step")
######################################################
# Post trajectory calculations for composite objective
######################################################
p.model.add_subsystem('other_calcs', om.ExecComp("f=v[0]+v[1]", v={"shape":(2,)}))
p.model.connect('traj.phase0.timeseries.states:v', 'other_calcs.v', src_indices=[0,-1])
###################################################################
# Standard Dymos API for objective function - minimize final time.
###################################################################
phase.add_objective('time', loc='final')
###################################################################
# OpenMDAO API for objective functions that require post-dymos calculations
###################################################################
# p.model.add_objective('other_calcs.f')
# Set the driver.
p.driver = om.ScipyOptimizeDriver()
# Allow OpenMDAO to automatically determine total
# derivative sparsity pattern.
# This works in conjunction with partial derivative
# coloring to give a large speedup
p.driver.declare_coloring()
# Setup the problem
p.setup()
# Now that the OpenMDAO problem is setup, we can guess the
# values of time, states, and controls.
p.set_val('traj.phase0.t_duration', 2.0)
# States and controls here use a linearly interpolated
# initial guess along the trajectory.
p.set_val('traj.phase0.states:x',
phase.interp('x', ys=[0, 10]),
units='m')
p.set_val('traj.phase0.states:y',
phase.interp('y', ys=[10, 5]),
units='m')
p.set_val('traj.phase0.states:y_tot',
phase.interp('y', ys=[0, 0]),
units='m*s')
p.set_val('traj.phase0.states:v',
phase.interp('v', ys=[0, 5]),
units='m/s')
# constant initial guess for control
p.set_val('traj.phase0.controls:theta', 90, units='deg')
# Run the driver to solve the problem and generate default plots of
# state and control values vs time
dm.run_problem(p, make_plots=True, simulate=True)
print(p.get_val("traj.phase0.fixed_step.states:v"))
time = p.get_val("traj.phase0.fixed_step.time")[:,0]
dt = time[1:] - time[0:-1]
print("time", time)
# note that dymos keeps duplicates of values between each segment, so you see double values
# which gives 0 time steps between each segment
print("dt: ", dt)
print(p.get_val("traj.phase0.timeseries.states:y_tot"))
