Change parameters before running simulation with Dymos for monte-carlo analysis

Question

when you run a Dymos problem and get a solution, you can simulate it and see that actual trajectories matches the collocation solution given by IPOPT.

But let's say I want to take things a bit further, and check out what my trajectory looks like if the initial conditions are slightly different, or some parameter like thrust or aero coefficients are a bit different from what is assigned.

Ideally, you'd get you optimal trajectory, and then be able to run a number of simulations after with the perturbed / randomized parameters and get out a "sheath" of trajectories.

For example, I want to change the initial fuel load a bit and run that simulation, instead of the original one.

How would I be able to change what gets actually run by the Dymos traj.simulate() function?

Answer 1

The dymos.run_problem function will automatically run the main problem followed by a simulation.

You can use the simulate method of Trajectory to run a simulation and return the resulting problem.

Once this problem is returned, you can modify inputs (like parameters) and rerun the simulation.

Consider the brachistochrone problem with different gravitational accelerations. First, solve the nominal optimal control problem.

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(tol=1.0E-12)

t = dm.Radau(num_segments=15, order=3, compressed=True)

traj = dm.Trajectory()
phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
p.model.add_subsystem('traj0', traj)
traj.add_phase('phase0', phase)

phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

phase.add_state('x', fix_initial=True, fix_final=False)
phase.add_state('y', fix_initial=True, fix_final=False)

# Note that by omitting the targets here Dymos will automatically attempt to connect
# to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
phase.add_state('v', fix_initial=True, fix_final=False)

phase.add_control('theta',
                  continuity=True, rate_continuity=True,
                  units='deg', lower=0.01, upper=179.9)

phase.add_parameter('g', targets=['g'], units='m/s**2')

phase.add_boundary_constraint('x', loc='final', equals=10)
phase.add_boundary_constraint('y', loc='final', equals=5)

# Minimize time at the end of the phase
phase.add_objective('time_phase', loc='final', scaler=10)

p.setup()

phase.set_simulate_options(method='RK45')

p.set_val('traj0.phase0.t_initial', 0.0)
p.set_val('traj0.phase0.t_duration', 2.0)

p.set_val('traj0.phase0.states:x', phase.interp('x', [0, 10]))
p.set_val('traj0.phase0.states:y', phase.interp('y', [10, 5]))
p.set_val('traj0.phase0.states:v', phase.interp('v', [0, 9.9]))

p.set_val('traj0.phase0.controls:theta', phase.interp('theta', [5, 100]))
p.set_val('traj0.phase0.parameters:g', 9.80665)

p.run_driver()

t_sol = p.get_val('traj0.phase0.timeseries.time').copy()
x_sol = p.get_val('traj0.phase0.timeseries.states:x').copy()
y_sol = p.get_val('traj0.phase0.timeseries.states:y').copy()

Note we used run_driver here rather than the dymos.run_problem interface. Next, we want to simulate the trajectory at the nominal solution.

sim_prob = traj.simulate()

Now extract the solution from the timeseries. We use the copy() method, otherwise we'll be getting a numpy view of these values that will be change with later calls to simulate.

t_nom = sim_prob.get_val('traj0.phase0.timeseries.time').copy()
x_nom = sim_prob.get_val('traj0.phase0.timeseries.states:x').copy()
y_nom = sim_prob.get_val('traj0.phase0.timeseries.states:y').copy()

Now change the parameter g in the simulation problem, and rerun the simulation. Note that there's no closed loop control modeled in this example, so the controls will respond as if they are optimized for the nominal value of gravity.

sim_prob.set_val('traj0.phase0.parameters:g', 9.0)
sim_prob.run_model()
t_g9 = sim_prob.get_val('traj0.phase0.timeseries.time').copy()
x_g9 = sim_prob.get_val('traj0.phase0.timeseries.states:x').copy()
y_g9 = sim_prob.get_val('traj0.phase0.timeseries.states:y').copy()

And again for another value of g:

sim_prob.set_val('traj0.phase0.parameters:g', 11.0)
sim_prob.run_model()
t_g11 = sim_prob.get_val('traj0.phase0.timeseries.time').copy()
x_g11 = sim_prob.get_val('traj0.phase0.timeseries.states:x').copy()
y_g11 = sim_prob.get_val('traj0.phase0.timeseries.states:y').copy()

Finally, plotting everything:

import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 1, figsize=(12, 6));

axes[0].set_xlabel('time (sec)')
axes[1].set_xlabel('time (sec)')
axes[0].set_ylabel('x (m)')
axes[1].set_ylabel('y (m)')

axes[0].plot(t_sol, x_sol, 'o', label='solution')
axes[1].plot(t_sol, y_sol, 'o', label='solution')

axes[0].plot(t_nom, x_nom, '-', label='sim: nominal')
axes[1].plot(t_nom, y_nom, '-', label='sim: nominal')

axes[0].plot(t_g9, x_g9, '-', label='sim: g=9.0')
axes[1].plot(t_g9, y_g9, '-', label='sim: g=9.0')

axes[0].plot(t_g11, x_g11, '-', label='sim: g=11.0')
axes[1].plot(t_g11, y_g11, '-', label='sim: g=11.0')

axes[0].grid()
axes[1].legend()

axes[1].grid()
axes[0].legend()

There are probably a few things we can do as developers to help you out for this sort of workflow. First, we might want to change the return signature of run_problem to allow this behavior (although changing the name of the recorder might be necessary then).

We should also allow the simulate method to return an instantiated, setup Problem instance without actually executing run_model, in case the user plans on changing things.