Dymos: How can I set object to minimize a function value of state variables at the end of simulation time

Question

I want to use Dymos to solve the optimal control problem:

subject to the dynamic system:

I have two questions:

(1) How to set the object function (object1) of V(X(T)), where the value of state variables at final time point are the inputs of the function V.

I think the setting

phase.add_objective('V(X)', loc='final', scale=1)

will minimize the integral from t=0 to t=T of a variable whose rate is V(X), is that right?

(2) I think I can setup an output objectval2 with time change rate L(X, u) then set

phase.add_objective('objectval2', loc='final', scale=1)

then how do I let the optimizer know that I want to minimize the sum of objectval1 + objectval2, instead of objectval1*objectval2 or any other function of objectval1 and objectval2.

I read the tutorial and I didn't find an example similar to my case.

#################

Comments: Thanks. I cannot add the following clarifying questions at the comments to the answer, it said these are too long. So I added them here.

I thinks the key code here is:

p.model.add_subsystem('other_calcs', om.ExecComp("f=a+b"))
p.model.connect('traj.phase0.timeseries.states:v', 'other_calcs.a', src_indices=[-1])
p.model.connect('traj.pahse0.timeseries.states:v', 'other_calcs.b', src_indices=[0])

I saw this in the tutorial and I have 3 questions:

(1) For the keyword input src_indices, I think src_indices=[-1] passes the last value of v(v[-1]) to a and src_indices=[0] passes the first value of v(v[0]) to b. I think the time point at v[-1] depends on the number of nodes, which decided by parameter num_segments and order. So does transcription function always set the last node position to the last time point? For example, if I solve a control problem from t=0 to t=2000, no matter how I set the num_segments and order, e.g. num_segments=2 and order=3 or num_segments=20 and order=5, v[-1] will always be the value of v at t=2000 and v[0] will always be the value of v at t=0?

(2) Can I add one or more specific time points to the nodes? For example, if I solve a control problem from t=0 to t=2000, I want to minimize the sum of the value of v at t=500, t=1000 and t=2000, can I required the nodes must includes the time points at t=500, t=1000 and t=2000, no matter how I set the num_segments and order?

(3) Last question is about the objective function setting I asked first time. In the example, you commented out # phase.add_objective('time', loc='final'). What will the optimizer do if I set the code:

phase.add_objective('x', loc='final')
p.model.add_objective('other_calcs.f')

Does the optimizer minimizes the sum of x integral from t=0 to t=T and v(T)? If not, how can I set the optimizer to minimizes the sum of x integral from t=0 to t=T and v(T)?

Justin Gray · Accepted Answer

(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"))

Dymos: How can I set object to minimize a function value of state variables at the end of simulation time

Answers (2)

Related Questions