Skip to content

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.

Vehicle Position States

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()
p.model.add_subsystem('traj', subsys=traj)

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

traj.add_phase(name='phase0', phase=phase)

# 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
phase.add_state('t', fix_initial=True, fix_final=False, units='s', lower=0,
                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
phase.add_state('alpha', fix_initial=False, fix_final=False, units='rad',
                rate_source='dalpha_ds', targets=['alpha'],
                ref=0.15)  # vehicle heading angle with respect to centerline
phase.add_state('lambda', fix_initial=False, fix_final=False, units='rad',
                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)
phase.add_state('omega', fix_initial=False, fix_final=False, units='rad/s',
                rate_source='domega_ds', targets=['omega'], ref=0.3)  # yaw rate
phase.add_state('ax', fix_initial=False, fix_final=False, units='m/s**2',
                rate_source='dax_ds', targets=['ax'], ref=8)  # longitudinal acceleration
phase.add_state('ay', fix_initial=False, fix_final=False, units='m/s**2',
                rate_source='day_ds', targets=['ay'], ref=8)  # lateral acceleration

# Define Controls
phase.add_control(name='delta', units='rad', lower=None, upper=None, fix_initial=False,
                  fix_final=False, targets=['delta'], ref=0.04)  # steering angle
phase.add_control(name='thrust', units=None, fix_initial=False, fix_final=False, targets=[
    '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
phase.add_path_constraint('c_rr', upper=1)
phase.add_path_constraint('c_rl', upper=1)
phase.add_path_constraint('c_fr', upper=1)
phase.add_path_constraint('c_fl', upper=1)

# Some of the vehicle design parameters are available to set here. Other parameters can
# be found in their respective ODE files.
phase.add_parameter('M', val=800.0, units='kg', opt=False,
                    targets=['car.M', 'tire.M', 'tireconstraint.M', 'normal.M'],
                    dynamic=False)  # vehicle mass
phase.add_parameter('beta', val=0.62, units=None, opt=False, targets=['tire.beta'],
                    dynamic=False)  # brake bias
phase.add_parameter('CoP', val=1.6, units='m', opt=False, targets=['normal.CoP'],
                    dynamic=False)  # center of pressure location
phase.add_parameter('h', val=0.3, units='m', opt=False, targets=['normal.h'],
                    dynamic=False)  # center of gravity height
phase.add_parameter('chi', val=0.5, units=None, opt=False, targets=['normal.chi'],
                    dynamic=False)  # roll stiffness
phase.add_parameter('ClA', val=4.0, units='m**2', opt=False, targets=['normal.ClA'],
                    dynamic=False)  # downforce coefficient*area
phase.add_parameter('CdA', val=2.0, units='m**2', opt=False, targets=['car.CdA'],
                    dynamic=False)  # drag coefficient*area

# Minimize final time.
# note that we use the 'state' time instead of Dymos 'time'
phase.add_objective('t', loc='final')

# Add output timeseries
phase.add_timeseries_output('*')

# Link the states at the start and end of the phase in order to ensure a continous lap
traj.link_phases(phases=['phase0', 'phase0'],
                 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'),
          units='rad/s')
p.set_val('traj.phase0.states:alpha', phase.interpolate(ys=[0.0, 0.0], nodes='state_input'),
          units='rad')
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',
          phase.interpolate(ys=[0.0, 0.0], nodes='control_input'), units='rad')
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,
                       pad=0.04)  # add colorbar

    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))
plt.subplots_adjust(right=0.82, bottom=0.14, top=0.97, left=0.07)

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 '
--- Linkage Report [traj] ---
    --- 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