Cyclic constraints

[pwmdebuck]

Summary of Issue

For trajectory optimization of cyclic events (such as race car laps), it would be helpful to be able to define constraints that ensure that the start and end states are identical

Issue Type

• [] Bug
• [x] Enhancement

Description

I have been working on trajectory optimization for race cars using Dymos. I would first like to thank you for the great tool that you have put together.

I was wondering if it would be possible to define a boundary constraint that links the initial and final states together. For a cyclic event like a race-car going around a lap I need the states and controls such as velocity, thrust, and steering angle to be the same at the start and finish (see barca_track_vel_annotated.pdf for an example result).

Currently I can do this by manually setting the initial and final values to be the same (fix_initial/fix_final=True) but this is an iterative process. Ideally I would be able to constrain the initial and final values to be the same, but not a fixed value. Is this currently possible in Dymos?

Additionally, I would be glad to provide the problem setup if interested. I am able to find an optimal trajectory for the race car with active aerodynamic and four wheel drive controls.

Environment

Operating System: OS X 10.15.5 Python environment: Python 3.8.5 Packages: dymos @ git+https://github.com/OpenMDAO/dymos.git@5d7e79ab8a89d0c47c3d75aabe299d0a36546fa8 mpi4py==3.0.3 numpy==1.19.1 openmdao==3.2.1 pyoptsparse==2.1.5 scipy==1.5.2

[JustinSGray]

Thanks for your kind words. I happen to like racing quite a lot ---I used to own a miata--- so Id love to see some of your results! Perhaps we could adapt your models to go-karting and use it to optimize our path on the local go-kart track!

There is not a built in way in Dymos to setup period boundary conditions. At least I don't think there is, but @robfalck can correct me if Im wrong. However, its pretty easy to create them using normal OpenMDAO APIs in combination with the timeseries outputs from the phases.

here is a sample of how to to do it in the brachistochrone problem. I modified it so that the final y state value must match the initial one. In this example fix_initial=True, but that does not fundamentally change anything. If your problem is well posed without fixed boundary conditions, the same technique will give you periodic bcs.

The 4 new lines are these:

p.model.add_subsystem('periodic_bc', om.ExecComp('bc_defect=final-initial'))
p.model.connect('traj0.phase0.timeseries.states:y', 'periodic_bc.initial', src_indices=0)
p.model.connect('traj0.phase0.timeseries.states:y', 'periodic_bc.final', src_indices=-1)
p.model.add_constraint('periodic_bc.bc_defect', equals=0)

Here is a full runable script:

import matplotlib
# matplotlib.use('Agg')
import matplotlib.pyplot as plt

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.pyOptSparseDriver()
p.driver.options['optimizer'] = 'SLSQP'
p.driver.declare_coloring(tol=1.0E-12)

traj = dm.Trajectory()
traj.add_phase('phase0', phase)
p.model.add_subsystem('traj0', traj)

t = dm.Radau(num_segments=8,order=3)

phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)

# p.model.add_subsystem('traj0', traj)

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

phase.add_state('x', rate_source=BrachistochroneODE.states['x']['rate_source'],
                units=BrachistochroneODE.states['x']['units'],
                fix_initial=True, fix_final=False)
phase.add_state('y', rate_source=BrachistochroneODE.states['y']['rate_source'],
                units=BrachistochroneODE.states['y']['units'],
                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', rate_source=BrachistochroneODE.states['v']['rate_source'],
                units=BrachistochroneODE.states['v']['units'],
                fix_initial=True, fix_final=False)

phase.add_control('theta', targets=BrachistochroneODE.parameters['theta']['targets'],
                  continuity=True, rate_continuity=True,
                  units='deg', lower=0.01, upper=179.9)

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

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

##############################################################################
# remove this and to use a periodic boundary condition instead
##############################################################################
# phase.add_boundary_constraint('y', loc='final', equals=5)

##############################################################################
# Use normal OpenMDAO APIs to create a periodic boundary condition
##############################################################################
p.model.add_subsystem('periodic_bc', om.ExecComp('bc_defect=final-initial'))
p.model.connect('traj0.phase0.timeseries.states:y', 'periodic_bc.initial', src_indices=0)
p.model.connect('traj0.phase0.timeseries.states:y', 'periodic_bc.final', src_indices=-1)
p.model.add_constraint('periodic_bc.bc_defect', equals=0)

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

p.setup(check=['unconnected_inputs'])

p['traj0.phase0.t_initial'] = 0.0
p['traj0.phase0.t_duration'] = 2.0

p['traj0.phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input')
p['traj0.phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input')
p['traj0.phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input')
p['traj0.phase0.controls:theta'] = phase.interpolate(ys=[5, 100], nodes='control_input')
p['traj0.phase0.parameters:g'] = 9.80665

dm.run_problem(p, run_driver=True)

exp_out = traj.simulate()

fig, ax = plt.subplots()
fig.suptitle('Brachistochrone Solution')

x_imp = p.get_val('traj0.phase0.timeseries.states:x')
y_imp = p.get_val('traj0.phase0.timeseries.states:y')

x_exp = exp_out.get_val('traj0.phase0.timeseries.states:x')
y_exp = exp_out.get_val('traj0.phase0.timeseries.states:y')

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('traj0.phase0.timeseries.time_phase')
y_imp = p.get_val('traj0.phase0.timeseries.controls:theta')

x_exp = exp_out.get_val('traj0.phase0.timeseries.time_phase')
y_exp = exp_out.get_val('traj0.phase0.timeseries.controls:theta')

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 (rad)')
ax.grid(True)
ax.legend(loc='lower right')

plt.show()

It produces a trajectory like this:

[robfalck]

Justin's example is a good example of how to accomplish this using native OpenMDAO machinery.

It's not well documented (we'll work on that) but you can use the Trajectory.link_phases functionality to "link" a phase with itself. In the example Justin posted above, if you omit the four lines he inserted and replace them with the following:

traj.link_phases(phases=['phase0', 'phase0'], vars=['y'], locs=('++', '--'))

Then the value of y at the end of phase0 ('++') is linked via constraint to the value of y at the beginning ('--') of phase0. This will produce the same plot Justin shows above.

We're always looking for example cases that demonstrate the capabilities of Dymos. Would you be interested in including your problem setup as an example problem in Dymos?

[pwmdebuck]

Thanks for the quick answers! If you have a specific go-kart track in mind I can create a track that resembles it and approximate some go-kart parameters and see what happens!

I will try the phase linking option first, thanks for the suggestions!

I would be happy to help with an example case. I have uploaded my scripts + results here. The runscript is carODEsolver.py and all the subsystems are connected in combinedODE.py. I will go through some steps to make the scripts more presentable/intuitive.

[robfalck]

@pwmdebuck Thanks for the example, it's really impressive...especially converging with 1300 segments! Do you happen to have a description for the states and controls used? I'm going to write a documentation page for this case (using the oval track for more simplicity), but I need to know how to describe the variables, and which variables should be linked via cyclic constraints.

[pwmdebuck]

Hi Rob,

Thanks a lot! Graph coloring is obviously essential at 1300 segments, but once it is computed this problem is actually quite well behaved! It works especially well with IPOPT which makes sense to me because the process of lowering the barrier parameter, finding the new limit, and repeating until sufficiently converged resembles a similar process a race car driver might go through when trying to minimize their lap time. If it is of interest, I have some comparison data between SNOPT and IPOPT for numbers of segments and convergence criteria.

This problem has been submitted as a paper to Vehicle System Dynamics, which is under review. I'm not familiar with how much I am allowed to share while it's under review, so I won't link the pdf right now. However, I have two tables which describe the states, controls, and vehicle parameters. The symbols match up with the variable names in the code. There is also a figure which describes the states that govern the vehicle position relative to the track.

I linked all the states (7) using your traj.link_phases() suggestion above. Linking the controls is not really needed since it generally follows from the linked states.

I have made some changes since I uploaded my scripts here. If you want I would be happy to clean it such that only the essentials remain. The oval track could be substituted for a very simple kart-like track for more interesting results.

Finally (sorry for the long message), I used this problem to construct a design of experiments given varying key parameters such as the lift coefficient and maximum engine power. I fit a surrogate model to this data which results in an interactive lap simulation. Check it out here: http://www.formulae.one/lapsimulation. Let me know what you think!

Thanks! Pieter

In this figure the curved line represents the centerline of the track with curvature kappa. n is the normal distance to the centerline, and s is the distance along the centerline.

[robfalck]

This looks great! This is also the first example we would have that replaces time with some other variable, so it would be really valuable for our documentation. I would really appreciate it if you could clean up a case - I'd prefer the oval track for the sake of simplicity and having our CI tests run in reasonable time. If you have something reasonable already I don't mind cleaning it up a bit before we put it on the examples.

The website looks great, and I'm looking forward to reading your paper.

[pwmdebuck]

Hi Rob,

I've cleaned up all the files for an oval track case, leaving only what is necessary for the optimization and basic plots. It can be found here: https://github.com/pwmdebuck/dymos-1/tree/master/dymos/examples/racecar

My Dymos fork is quite a bit behind so I hope that it works on the current build!

[robfalck]

Thank you Pieter,

Even if it's a bit behind I'll be able to take it from here, thank you.