# Race car Lap Simulation¶

Things you'll learn through this example

• Optimizing trajectories for cyclic events (such as race car laps)
• Linking a phase with itself so it is ensured that start and end states are identical
• Using an independent variable other than time

The race track problem

The race track problem is different from other dymos example problems in that it is a cyclic problem. The solution must be such that the start and end states are identical. The problem for this example is given a race track shape, what is the optimum values for the controls such that the time elapsed to complete the track is a minimum.

## Background of the problem¶

This example is based on the work of Pieter de Buck The goal of his problem was to find an optimal trajectory for the race car with active aerodynamic and four wheel drive controls.

Pieter posted an issue to the Dymos repository asking how to use dymos for cyclic events such as this. His solution involved using an iterative solution, which worked but dymos includes a feature that allows for simple and more efficient solution. This example shows how to handle this situation using simple dymos calls without iteration.

If you are interested, Pieter and Adeline Shin have created an interactive visualization using his version of the code.

While Pieter's code includes several actual race tracks shapes, for simplicity and speed of solution, this dymos example uses an oval track. If you are interested in trying other race tracks, they are defined in the tracks.py file.

It is easy to use a different track. For example, in the test_doc_racecar.py file,

modify the lines:

from dymos.examples.racecar.tracks import ovaltrack
track = ovaltrack


For example, to use the Monaco track,

from dymos.examples.racecar.tracks import Monaco
track = Monaco


## Optimizing trajectories for cyclic events¶

For trajectory optimization of cyclic events (such as race car laps), there needs to be a way that ensure that the start and end states are identical.

In this race car example, the states and controls such as velocity, thrust, and steering angle need to be the same at the start and finish.

Linking the controls is not really needed since it generally follows from the linked states.

There is a simple way to link the states of cyclic events in dymos. You can use the Trajectory.link_phases functionality to "link" a phase with itself. In this example this is accomplished with the following code which links all 7 state variables.

traj.link_phases(phases=['phase0', 'phase0'], vars=['V','n','alpha','omega','lambda','ax','ay'], locs=('final', 'initial'))


This code causes the value of V at the end of phase0 ('final') to be linked via constraint to the value of V at the beginning ('initial') of phase0. This call does the same for the other 6 state variables as well.

See the Trajectory API page for more information about the link_phases method.

## Using an independent variable other than time¶

In all the other dymos examples, time is the independent variable. In this race car example, the independent variable is the arc length along the track. Integration is done from the start of the track to the end. The “typical” equations of motion are computed, so the rate of change of x wrt time ($\frac{dx}{dt}$).

By dividing the rates of change w.r.t. time by ($\frac{ds}{dt}$), the rates of change of the states w.r.t. track arclength traversed are computed.

## Definition of the problem¶

Here are some figures which help define the problem.

This table describes the controls, states, and constraints. The symbols match up with the variable names in the code shown at the bottom of this page.

###### Table 1. Optimization problem formulation¶
Symbol Lower Upper Units
Minimize Time t s
With respect to
Controls Front left thrust $T_{fl}$
Front right thrust $T_{fr}$
Rear left thrust $T_{rl}$
Rear right thrust $T_{rr}$
Rear wing flap angle $\gamma$ 0 50 deg
Steering angle $\delta$ rad
States Speed V $ms^{-1}$
Longitudinal acceleration $a_x$ $ms^{-2}$
Lateral acceleration $a_y$ $ms^{-2}$
Normal distance to centerline n -4 4 m
Angle relative to centerline $\alpha$ rad
Slip angle $\lambda$ rad
Yaw rate $\Omega$ $rad^{-1}$
Subject to Front left adherence $c_{fl}$ 0 1
Front right adherence $c_{fr}$ 0 1
Rear left adherence $c_{rl}$ 0 1
Rear right adherence $c_{rr}$ 0 1
Front left power $P_{fl}$ 75 kW
Front right power $P_{fr}$ 75 kW
Rear left power $P_{rl}$ 75 kW
Rear right power $P_{rr}$ 75 kW
Periodic state constraints
System dynamics
Track Layout

The next figure defines the vehicle parameters. Again, the symbols match up with the variable names in the code.

###### Table 2. Vehicle parameters¶
Parameter Value Units Description
$M$ 1184 $kg$ Vehicle mass
$a$ 1.404 $m$ CoG to front axle distance
$b$ 1.356 $m$ CoG to rear axle distance
$t_w$ 0.807 $m$ Half track width
$h$ 0.4 $m$ CoG height
$I_z$ 1775 $kg\ m^{2}$ Yaw inertia
$\beta$ 0.62 - Brake balance
$\chi$ 0.5 - Roll stiffness
$\rho$ 1.2 $kg\ m^{-3}$ Air density
$\mu_{0}^{x}$ 1.68 - Longitudinal base friction coefficient
$\mu_{0}^{y}$ 1.68 - Lateral base friction coefficient
$K_{\mu}$ -0.5 - Tire load sensitivity
$K_{\lambda}$ 44 - Tire lateral stiffness
$\tau_{a_{x}}$ 0.2 $s$ Longitudinal load transfer time constant
$\tau_{a_{y}}$ 0.2 $s$ Lateral load transfer time constant
$S_{w}$ 0.8 $m^{2}$ Wing planform area
$CoP$ 1.404 $m$ Center of pressure to front axle distance

There is also a figure which describes the states that govern the vehicle position relative to the track.

Here is the code to solve this problem:

import numpy as np
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal
import dymos as dm
import matplotlib.pyplot as plt
import matplotlib as mpl

from dymos.examples.racecar.combinedODE import CombinedODE
from dymos.examples.racecar.spline import (get_spline, get_track_points, get_gate_normals,
reverse_transform_gates, set_gate_displacements,
transform_gates, )
from dymos.examples.racecar.linewidthhelper import linewidth_from_data_units
from dymos.examples.racecar.tracks import ovaltrack  # track curvature imports

# change track here and in curvature.py. Tracks are defined in tracks.py
track = ovaltrack

# generate nodes along the centerline for curvature calculation (different
# than collocation nodes)
points = get_track_points(track)

# fit the centerline spline.
finespline, gates, gatesd, curv, slope = get_spline(points, s=0.0)
# by default 10000 points
s_final = track.get_total_length()

# Define the OpenMDAO problem
p = om.Problem(model=om.Group())

# Define a Trajectory object
traj = dm.Trajectory()

# Define a Dymos Phase object with GaussLobatto Transcription
phase = dm.Phase(ode_class=CombinedODE,
transcription=dm.GaussLobatto(num_segments=50, order=3, compressed=True))

# Set the time options, in this problem we perform a change of variables. So 'time' is
# actually 's' (distance along the track centerline)
# This is done to fix the collocation nodes in space, which saves us the calculation of
# the rate of change of curvature.
# The state equations are written with respect to time, the variable change occurs in
# timeODE.py
phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=s_final,
targets=['curv.s'], units='m', duration_ref=s_final,
duration_ref0=10)

# Define states
rate_source='dt_ds', ref=100)  # time
phase.add_state('n', fix_initial=False, fix_final=False, units='m', upper=4.0, lower=-4.0,
rate_source='dn_ds', targets=['n'],
ref=4.0)  # normal distance to centerline. The bounds on n define the
# width of the track
phase.add_state('V', fix_initial=False, fix_final=False, units='m/s', ref=40, ref0=5,
rate_source='dV_ds', targets=['V'])  # velocity
rate_source='dalpha_ds', targets=['alpha'],
ref=0.15)  # vehicle heading angle with respect to centerline
rate_source='dlambda_ds', targets=['lambda'],
ref=0.01)  # vehicle slip angle, or angle between the axis of the vehicle
# and velocity vector (all cars drift a little)
rate_source='domega_ds', targets=['omega'], ref=0.3)  # yaw rate
rate_source='dax_ds', targets=['ax'], ref=8)  # longitudinal acceleration
rate_source='day_ds', targets=['ay'], ref=8)  # lateral acceleration

# Define Controls
fix_final=False, targets=['delta'], ref=0.04)  # steering angle
'thrust'])  # the thrust controls the longitudinal force of the rear tires and is
# positive while accelerating, negative while braking

# Performance Constraints
pmax = 960000  # W
phase.add_path_constraint('power', upper=pmax, ref=100000)  # engine power limit

# The following four constraints are the tire friction limits, with 'rr' designating the
# rear right wheel etc. This limit is computed in tireConstraintODE.py

# Some of the vehicle design parameters are available to set here. Other parameters can
# be found in their respective ODE files.
targets=['car.M', 'tire.M', 'tireconstraint.M', 'normal.M'],
dynamic=False)  # vehicle mass
dynamic=False)  # brake bias
dynamic=False)  # center of pressure location
dynamic=False)  # center of gravity height
dynamic=False)  # roll stiffness
dynamic=False)  # downforce coefficient*area
dynamic=False)  # drag coefficient*area

# Minimize final time.
# note that we use the 'state' time instead of Dymos 'time'

# Link the states at the start and end of the phase in order to ensure a continous lap
vars=['V', 'n', 'alpha', 'omega', 'lambda', 'ax', 'ay'],
locs=('final', 'initial'))

# Set the driver. IPOPT or SNOPT are recommended but SLSQP might work.
p.driver = om.pyOptSparseDriver(optimizer='IPOPT')

p.driver.opt_settings['mu_init'] = 1e-3
p.driver.opt_settings['max_iter'] = 500
p.driver.opt_settings['acceptable_tol'] = 1e-3
p.driver.opt_settings['constr_viol_tol'] = 1e-3
p.driver.opt_settings['compl_inf_tol'] = 1e-3
p.driver.opt_settings['acceptable_iter'] = 0
p.driver.opt_settings['tol'] = 1e-3
p.driver.opt_settings['nlp_scaling_method'] = 'none'
p.driver.opt_settings['print_level'] = 5
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'  # for faster convergence
p.driver.options['print_results'] = False

# Allow OpenMDAO to automatically determine our sparsity pattern.
# Doing so can significant speed up the execution of Dymos.
p.driver.declare_coloring()

# Setup the problem
p.setup(check=True)  # force_alloc_complex=True
# Now that the OpenMDAO problem is setup, we can set the values of the states.

# States
p.set_val('traj.phase0.states:V', phase.interpolate(ys=[20, 20], nodes='state_input'),
units='m/s')  # non-zero velocity in order to protect against 1/0 errors.
p.set_val('traj.phase0.states:lambda',
phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
units='rad')  # all other states start at 0
p.set_val('traj.phase0.states:omega', phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
p.set_val('traj.phase0.states:alpha', phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
p.set_val('traj.phase0.states:ax', phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
units='m/s**2')
p.set_val('traj.phase0.states:ay', phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
units='m/s**2')
p.set_val('traj.phase0.states:n', phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
units='m')
p.set_val('traj.phase0.states:t', phase.interpolate(ys=[0.0, 100.0], nodes='state_input'),
units='s')  # initial guess for what the final time should be

# Controls
p.set_val('traj.phase0.controls:delta',
p.set_val('traj.phase0.controls:thrust',
phase.interpolate(ys=[0.1, 0.1], nodes='control_input'),
units=None)  # a small amount of thrust can speed up convergence

dm.run_problem(p, run_driver=True)
print('Optimization finished')

# Test this example in Dymos' continuous integration process

# Get optimized time series
n = p.get_val('traj.phase0.timeseries.states:n')
s = p.get_val('traj.phase0.timeseries.time')
V = p.get_val('traj.phase0.timeseries.states:V')
thrust = p.get_val('traj.phase0.timeseries.controls:thrust')
delta = p.get_val('traj.phase0.timeseries.controls:delta')
power = p.get_val('traj.phase0.timeseries.power', units='W')

print("Plotting")

# We know the optimal distance from the centerline (n). To transform this into the racing
# line we fit a spline to the displaced points. This will let us plot the racing line in
# x/y coordinates
normals = get_gate_normals(finespline, slope)
newgates = []
newnormals = []
newn = []
for i in range(len(n)):
index = ((s[i] / s_final) * np.array(finespline).shape[1]).astype(
int)  # interpolation to find the appropriate index
if index[0] == np.array(finespline).shape[1]:
index[0] = np.array(finespline).shape[1] - 1
if i > 0 and s[i] == s[i - 1]:
continue
else:
newgates.append([finespline[0][index[0]], finespline[1][index[0]]])
newnormals.append(normals[index[0]])
newn.append(n[i][0])

newgates = reverse_transform_gates(newgates)
displaced_gates = set_gate_displacements(newn, newgates, newnormals)
displaced_gates = np.array((transform_gates(displaced_gates)))

npoints = 1000
# fit the racing line spline to npoints
displaced_spline, gates, gatesd, curv, slope = get_spline(displaced_gates, 1 / npoints, 0)

plt.rcParams.update({'font.size': 12})

def plot_track_with_data(state, s):
# this function plots the track
state = np.array(state)[:, 0]
s = np.array(s)[:, 0]
s_new = np.linspace(0, s_final, npoints)

# Colormap and norm of the track plot
cmap = mpl.cm.get_cmap('viridis')
norm = mpl.colors.Normalize(vmin=np.amin(state), vmax=np.amax(state))

fig, ax = plt.subplots(figsize=(15, 6))
# establishes the figure axis limits needed for plotting the track below
plt.plot(displaced_spline[0], displaced_spline[1], linewidth=0.1, solid_capstyle="butt")

plt.axis('equal')
# the linewidth is set in order to match the width of the track
plt.plot(finespline[0], finespline[1], 'k',
linewidth=linewidth_from_data_units(8.5, ax),
solid_capstyle="butt")
plt.plot(finespline[0], finespline[1], 'w', linewidth=linewidth_from_data_units(8, ax),
solid_capstyle="butt")  # 8 is the width, and the 8.5 wide line draws 'kerbs'
plt.xlabel('x (m)')
plt.ylabel('y (m)')

# plot spline with color
for i in range(1, len(displaced_spline[0])):
s_spline = s_new[i]
index_greater = np.argwhere(s >= s_spline)[0][0]
index_less = np.argwhere(s < s_spline)[-1][0]

x = s_spline
xp = np.array([s[index_less], s[index_greater]])
fp = np.array([state[index_less], state[index_greater]])
interp_state = np.interp(x, xp,
fp)  # interpolate the given state to calculate the color

# calculate the appropriate color
state_color = norm(interp_state)
color = cmap(state_color)
color = mpl.colors.to_hex(color)

# the track plot consists of thousands of tiny lines:
point = [displaced_spline[0][i], displaced_spline[1][i]]
prevpoint = [displaced_spline[0][i - 1], displaced_spline[1][i - 1]]
if i <= 5 or i == len(displaced_spline[0]) - 1:
plt.plot([point[0], prevpoint[0]], [point[1], prevpoint[1]], color,
linewidth=linewidth_from_data_units(1.5, ax), solid_capstyle="butt",
antialiased=True)
else:
plt.plot([point[0], prevpoint[0]], [point[1], prevpoint[1]], color,
linewidth=linewidth_from_data_units(1.5, ax),
solid_capstyle="projecting", antialiased=True)

clb = plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), fraction=0.02,

if np.array_equal(state, V[:, 0]):
clb.set_label('Velocity (m/s)')
elif np.array_equal(state, thrust[:, 0]):
clb.set_label('Thrust')
elif np.array_equal(state, delta[:, 0]):
clb.set_label('Delta')

plt.tight_layout()
plt.grid()

# Create the plots
plot_track_with_data(V, s)
plot_track_with_data(thrust, s)
plot_track_with_data(delta, s)

# Plot the main vehicle telemetry
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(15, 8))

# Velocity vs s
axes[0].plot(s,
p.get_val('traj.phase0.timeseries.states:V'), label='solution')

axes[0].set_xlabel('s (m)')
axes[0].set_ylabel('V (m/s)')
axes[0].grid()
axes[0].set_xlim(0, s_final)

# n vs s
axes[1].plot(s,
p.get_val('traj.phase0.timeseries.states:n', units='m'), label='solution')

axes[1].set_xlabel('s (m)')
axes[1].set_ylabel('n (m)')
axes[1].grid()
axes[1].set_xlim(0, s_final)

# throttle vs s
axes[2].plot(s, thrust)

axes[2].set_xlabel('s (m)')
axes[2].set_ylabel('thrust')
axes[2].grid()
axes[2].set_xlim(0, s_final)

# delta vs s
axes[3].plot(s,
p.get_val('traj.phase0.timeseries.controls:delta', units=None),
label='solution')

axes[3].set_xlabel('s (m)')
axes[3].set_ylabel('delta')
axes[3].grid()
axes[3].set_xlim(0, s_final)

plt.tight_layout()

# Performance constraint plot. Tire friction and power constraints
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(15, 4))

axes.plot(s,
p.get_val('traj.phase0.timeseries.c_fl', units=None), label='c_fl')
axes.plot(s,
p.get_val('traj.phase0.timeseries.c_fr', units=None), label='c_fr')
axes.plot(s,
p.get_val('traj.phase0.timeseries.c_rl', units=None), label='c_rl')
axes.plot(s,
p.get_val('traj.phase0.timeseries.c_rr', units=None), label='c_rr')

axes.plot(s, power / pmax, label='Power')

axes.legend(bbox_to_anchor=(1.04, 0.5), loc="center left")
axes.set_xlabel('s (m)')
axes.set_ylabel('Performance constraints')
axes.grid()
axes.set_xlim(0, s_final)

plt.show()

/home/travis/build/OpenMDAO/dymos/dymos/phase/phase.py:1555: UserWarning: Phase time options have no effect because fix_duration=True or input_duration=True for phase 'phase0': duration_ref, duration_ref0
warnings.warn('Phase time options have no effect because fix_duration=True or '
/home/travis/build/OpenMDAO/dymos/dymos/phase/phase.py:1555: UserWarning: Phase time options have no effect because fix_duration=True or input_duration=True for phase 'phase0': duration_ref, duration_ref0
warnings.warn('Phase time options have no effect because fix_duration=True or '
--- phase0 - phase0 ---
V      [final]  ==  V      [initial]
n      [final]  ==  n      [initial]
alpha  [final]  ==  alpha  [initial]
omega  [final]  ==  omega  [initial]
lambda [final]  ==  lambda [initial]
ax     [final]  ==  ax     [initial]
ay     [final]  ==  ay     [initial]
----------------------------
INFO: checking out_of_order
INFO: checking system
INFO: checking solvers
INFO: checking dup_inputs
INFO: checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
INFO: checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO: checking out_of_order
INFO: checking system
INFO: checking solvers
INFO: checking dup_inputs
INFO: checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
INFO: checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO: checking out_of_order
INFO: checking system
INFO: checking solvers
INFO: checking dup_inputs
INFO: checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
INFO: checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
Full total jacobian was computed 3 times, taking 8.269288 seconds.
Total jacobian shape: (1256, 609)

Jacobian shape: (1256, 609)  ( 5.55% nonzero)
FWD solves: 67   REV solves: 0
Total colors vs. total size: 67 vs 609  (89.0% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 8.269288 sec.
Time to compute coloring: 2.569503 sec.
/home/travis/build/OpenMDAO/dymos/OpenMDAO/openmdao/core/total_jac.py:1718: UserWarning:Constraints or objectives ['traj.phases.phase0.path_constraints.path:c_fr', 'traj.phases.phase0.path_constraints.path:c_fl'] cannot be impacted by the design variables of the problem.
/home/travis/build/OpenMDAO/dymos/OpenMDAO/openmdao/core/total_jac.py:1718: UserWarning:Constraints or objectives ['traj.phases.phase0.path_constraints.path:c_fr', 'traj.phases.phase0.path_constraints.path:c_fl'] cannot be impacted by the design variables of the problem.
Optimization finished
Plotting