Comparison of Numerical Integration Schemes (Adaptive Runge-Kutta)
Contents
Comparison of Numerical Integration Schemes (Adaptive Runge-Kutta)#
Thomas Bieniewicz#
Problem 1 of Homework 1 utilizes two ODE Numerical Integration Schemes are used: 1st Order Forward Euler and 4th Order Runge-Kutta.
These schemes appeared to work well enough for basic applications at a first glance but begged the question: What are the advantages and disadvantages of using other ODE Numerical Integrators with advanced features, such as an adaptie time step for calculations?
Importance of ODE Numerical Solvers/Integeators:
In practice we often plan in discrete time
When applying control methods in real world settings, digital controllers on computers are almost always used, which require discrete time and thus continuous time dynamics must be discretized for implementation
There are very few cases where discretizing a continuous time system can be done analytically, even fewer when strictly considering “real world settings” where it can be difficult to make assumptions about linearity and other simplifications to a problem.
Adaptive Runge-Kutta Method:
Concept
Before the integration an initial time step value, min value for change in x, max allowed relative difference in half step and full step, and min allowed relative difference in full step and double step
Every time the next value \(x_k+1\) is calculated a full step (same as regular Runge-Kutta), a half step, and double step are calculated. This is like moving in the same “direction” as with non-adaptive Runge-Kutta but at only half or double the “distance”
Then, we compare the difference in calcualted \(x_k+1\) between the full step and half step and full step and double step to determine if we should take a larger or shorter time step proceeding from the current time step.
Advantages
Potential savings in computation by specifying limits on the relative value of a full, half, and double step and defining the time step based upon those specifications as a variable. This allows the algorithm to take much larger steps when the dynamics are smooth and don’t require as high of a degree of resolution to capture its shape (function changing slowly with time); and take much smaller steps when a high degree of resiultion is required (function is changing a lot with very small changes in time)
In effect, the algorithm records more data points when they are needed to determine the behavior of the function and takes less data points to save on computation not needed to determine the behavior of the function. This becomes extrmely valuable when requiring a high degree of accuracy for simulations of long periods of time.
No new math is needed to apply this method; only new code to generate an actively changing time step (dt). The ODE Numerical Integrator should work with any system linear, nonlinear, or time varying by changing the continuous dynamics function.
Disadvantages
Because the dt (time step) is variable, the number of values in a contol sequence and their corresponding times are also variable, which makes defining an input sequence for u is impossible to do before the integration. This is not an issue if the control input u is constant for interval the numerical integral is evaluated over.
Instead, if an analytic function for the input u exists, a solution would be to slightly change the code so that the analytic function for u can be evaluated at each time t and then its output array could be fed as an argument to the dynamics functions.
import abc
from typing import Callable
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import functools
import cvxpy as cvx
For this example, the same code from HW1 Problem 1 (dynamic_unicycle_ode) is used to define the continuous time dynamics for easy comparison, and a slightly modified version of the runge_kutta_integrator function still utilizing the same 4th order Runge-Kutta scheme from that same problem is also used.
the control input u is held constant
the initial state \(x(0)\) = 0
system is simulated over 5 second duration (t = 0 to t = 5)
results are plotted with the individual data points marked
# Step 1 - define Continuous Dynamics (same as HW1-P1)
def dynamic_unicycle_ode(x, u, t):
x, y, theta, v = x # decompose state variables and control inputs
w, a = u
# still need to pass in t, its just not used because time invariant
xdot = v * jnp.cos(theta) # nonlinear state equations
ydot = v * jnp.sin(theta)
thetadot = w
vdot = a
state_dot = jnp.array([xdot,ydot,thetadot,vdot]) # derivative of state variables --> xdot = f(x,u,t)
return state_dot
# Step 2 - define RK4 integrator (converts continuous dyunamics to discrete dynamics)
def adaptive_runge_kutta_integrator(continuous_dynamics, dt, x, u, t):
f1 = continuous_dynamics(x,u,t) # still need to pass in t, not used in time invariant
f2 = continuous_dynamics(x + (dt/2)*f1, u, t + (dt/2))
f3 = continuous_dynamics(x + (dt/2)*f2, u, t + (dt/2))
f4 = continuous_dynamics(x + dt*f3, u, t + dt)
x_full = x + (dt/6) * (f1 + 2*f2 + 2*f3 + f4) # full step
x_double = x + (dt/3) * (f1 + 2*f2 + 2*f3 + f4) # double step
x_half = x + (dt/12) * (f1 + 2*f2 + 2*f3 + f4) # half step
return x_full, x_half, x_double # returns result of all step sizes
# Step 3 - define function to check differences in step magnitudes, and update dt and t
def step_compare(dt, x_steps, t):
x_full, x_half, x_double = x_steps # unpack step data
x_nxt = 0
# define tolerance specifications
x_min = 1e-3
dx_min = 0.008
dx_max = 0.01
dt_min = 1e-3
# prevent dt from becoming too small if x is small value
if abs(x_full[0]) < x_min or abs(x_full[1]) < x_min or abs(x_full[2]) < x_min or abs(x_full[3]) < x_min:
dt = dt_min
x_nxt = x_full
else:
# if full step is more than dx_max percent larger than half step, reduce dt (as long as dt > dt_min)
if jnp.linalg.norm(x_full) > x_min and (jnp.linalg.norm(x_full - x_half)/jnp.linalg.norm(x_full)) > dx_max and dt > dt_min:
dt = dt/2
x_nxt = x_half
# if double step is less that dx_min percent larger than full step, increase dt
elif jnp.linalg.norm(x_full) > x_min and (jnp.linalg.norm(x_full - x_double)/jnp.linalg.norm(x_full)) < dx_min:
dt = 2*dt
x_nxt = x_double
# if full step within defined percent ranges of x_half and x_double, keep dt constant
else:
x_nxt = x_full
t = t + dt # update t with whichever size of time step was taken
return x_nxt, dt, t
# Step 5 - create function to simulate dynamics on specified time interval for variable time step
def simulate(discrete_dynamics, continuous_dynamics, step_compare, initial_state, controls, initial_dt, tmax):
x = initial_state # set initial state
dt = initial_dt # set initial time step
t_data = np.zeros(1)
t = 0
while t <= tmax:
print(t,tmax)
if t == 0:
x_data = x # set initial value as first data point in data
t_data[0] = t
x_steps = discrete_dynamics(continuous_dynamics,dt,x,controls,t_data[-1]) # call dynamics to give full, half, and double steps
x_nxt, dt, t = step_compare(dt,x_steps,t_data[-1]) # calcualte x_nxt, dt for next step, current t
x_data = jnp.vstack((x_data,x_nxt)) # add new data point to data
t_data = jnp.append(t_data,t) # add new time value to data
x = x_nxt # set current state for next calculation # set current state for next calculation
return x_data
Now that the continuous dynamics, discrete dynamics (numerical integrator), function to determine whether to increase or decrease the time step size (dt) and, function to simulate the system have been written its time to specify needed initial values and simulate the system so we can plot the results.
Note: the simulate function prints the values of t and tmax repeatedly so you can watch the current time t appraoch tmax and stop the simulation
initial_state = jnp.array([0.0, 0.0, 0.0, 0.0]) # set initial state value
control = jnp.array([2.0, 1.0]) # constant control over the 5 second duration. (u)
tmax = 5 # specify duration for integration
controls = control # constant control over the 5 second duration. (u)
initial_dt = 0.1 # initial time step value for starting integration
# Call adaptive RK
xs_adaptRK = simulate(adaptive_runge_kutta_integrator,dynamic_unicycle_ode, step_compare, initial_state, controls, initial_dt,tmax)
# plot the trajectories
plt.plot(xs_adaptRK[:,0], xs_adaptRK[:,1], label="variable dt - Adaptive RK", color='r')
plt.plot(xs_adaptRK[:,0], xs_adaptRK[:,1], marker='+', color='b', alpha=0.45)
plt.legend()
plt.grid(alpha=0.4)
plt.axis("equal")
plt.xlabel("x [m]")
plt.ylabel("y [m]")
No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
0 5
0.001 5
0.002 5
0.003 5
0.004 5
0.0050000004 5
0.0060000005 5
0.0070000007 5
0.008 5
0.009000001 5
0.010000001 5
0.011000001 5
0.012000001 5
0.013000001 5
0.014000001 5
0.0150000015 5
0.016 5
0.017 5
0.018000001 5
0.019000001 5
0.020000001 5
0.021000002 5
0.022000002 5
0.023000002 5
0.024000002 5
0.025000002 5
0.026000002 5
0.028000003 5
0.030000003 5
0.032 5
0.034 5
0.036000002 5
0.038000003 5
0.040000003 5
0.042000003 5
0.044000003 5
0.046000004 5
0.048000004 5
0.050000004 5
0.052000005 5
0.054000005 5
0.056000005 5
0.058000006 5
0.060000006 5
0.062000006 5
0.064 5
0.066 5
0.067999996 5
0.06999999 5
0.07199999 5
0.073999986 5
0.07599998 5
0.07799998 5
0.079999976 5
0.08199997 5
0.08399997 5
0.085999966 5
0.08799996 5
0.08999996 5
0.091999955 5
0.09399995 5
0.09599995 5
0.097999945 5
0.09999994 5
0.10199994 5
0.103999935 5
0.10599993 5
0.10799993 5
0.109999925 5
0.11199992 5
0.11399992 5
0.115999915 5
0.11799991 5
0.11999991 5
0.121999905 5
0.1239999 5
0.1259999 5
0.1279999 5
0.1299999 5
0.13199991 5
0.13399991 5
0.13599992 5
0.13799992 5
0.13999993 5
0.14199993 5
0.14399993 5
0.14599994 5
0.14799994 5
0.14999995 5
0.15199995 5
0.15599994 5
0.15999994 5
0.16399993 5
0.16799992 5
0.17199992 5
0.17599991 5
0.1799999 5
0.1839999 5
0.18799989 5
0.19199988 5
0.19599988 5
0.19999987 5
0.20399986 5
0.20799986 5
0.21199985 5
0.21599984 5
0.21999983 5
0.22399983 5
0.22799982 5
0.23199981 5
0.23599981 5
0.2399998 5
0.2439998 5
0.24799979 5
0.2519998 5
0.2559998 5
0.2599998 5
0.26399982 5
0.26799983 5
0.27199984 5
0.27599984 5
0.27999985 5
0.28399986 5
0.28799987 5
0.29199988 5
0.29599988 5
0.2999999 5
0.3039999 5
0.3079999 5
0.31199992 5
0.31599993 5
0.31999993 5
0.32399994 5
0.32799995 5
0.33199996 5
0.33599997 5
0.33999997 5
0.34399998 5
0.348 5
0.352 5
0.356 5
0.36 5
0.36400002 5
0.36800003 5
0.37200004 5
0.37600005 5
0.38000005 5
0.38400006 5
0.38800007 5
0.39200008 5
0.3960001 5
0.4000001 5
0.4040001 5
0.4080001 5
0.4160001 5
0.42400008 5
0.43200007 5
0.44000006 5
0.44800004 5
0.45600003 5
0.46400002 5
0.472 5
0.48 5
0.48799998 5
0.49599996 5
0.50399995 5
0.51199996 5
0.52 5
0.528 5
0.536 5
0.544 5
0.55200005 5
0.56000006 5
0.5680001 5
0.5760001 5
0.5840001 5
0.5920001 5
0.60000014 5
0.60800016 5
0.6160002 5
0.6240002 5
0.6320002 5
0.6400002 5
0.64800024 5
0.65600026 5
0.6640003 5
0.6720003 5
0.6800003 5
0.6880003 5
0.69600034 5
0.70400035 5
0.71200037 5
0.7200004 5
0.7280004 5
0.7360004 5
0.74400043 5
0.75200045 5
0.76000047 5
0.7680005 5
0.7760005 5
0.7840005 5
0.79200053 5
0.80000055 5
0.80800056 5
0.8160006 5
0.8240006 5
0.8320006 5
0.8400006 5
0.84800065 5
0.85600066 5
0.8640007 5
0.8720007 5
0.8800007 5
0.8880007 5
0.89600074 5
0.90400076 5
0.9120008 5
0.9200008 5
0.9280008 5
0.9360008 5
0.94400084 5
0.95200086 5
0.9600009 5
0.9680009 5
0.9760009 5
0.9840009 5
1.000001 5
1.016001 5
1.032001 5
1.048001 5
1.0640011 5
1.0800011 5
1.0960011 5
1.1120012 5
1.1280012 5
1.1440012 5
1.1600013 5
1.1760013 5
1.1920013 5
1.2080014 5
1.2240014 5
1.2400014 5
1.2560015 5
1.2720015 5
1.2880015 5
1.3040016 5
1.3200016 5
1.3360016 5
1.3520017 5
1.3680017 5
1.3840017 5
1.4000018 5
1.4160018 5
1.4320018 5
1.4480019 5
1.4640019 5
1.4800019 5
1.496002 5
1.512002 5
1.528002 5
1.544002 5
1.5600021 5
1.5760021 5
1.5920022 5
1.6080022 5
1.6240022 5
1.6400023 5
1.6560023 5
1.6720023 5
1.6880023 5
1.7040024 5
1.7200024 5
1.7360024 5
1.7520025 5
1.7680025 5
1.7840025 5
1.8000026 5
1.8160026 5
1.8320026 5
1.8480027 5
1.8640027 5
1.8800027 5
1.8960028 5
1.9120028 5
1.9280028 5
1.9440029 5
1.9600029 5
1.9760029 5
1.992003 5
2.008003 5
2.024003 5
2.040003 5
2.056003 5
2.0720031 5
2.0880032 5
2.1040032 5
2.1200032 5
2.1360033 5
2.1520033 5
2.1680033 5
2.1840034 5
2.2000034 5
2.2160034 5
2.2320035 5
2.2480035 5
2.2640035 5
2.2800035 5
2.2960036 5
2.3120036 5
2.3280036 5
2.3440037 5
2.3600037 5
2.3760037 5
2.3920038 5
2.4080038 5
2.4240038 5
2.4400039 5
2.456004 5
2.472004 5
2.488004 5
2.504004 5
2.520004 5
2.536004 5
2.552004 5
2.5680041 5
2.5840042 5
2.6000042 5
2.6160042 5
2.6320043 5
2.6480043 5
2.6640043 5
2.6800044 5
2.6960044 5
2.7120044 5
2.7280045 5
2.7440045 5
2.7600045 5
2.7760046 5
2.7920046 5
2.8080046 5
2.8240047 5
2.8400047 5
2.8560047 5
2.8720047 5
2.8880048 5
2.9040048 5
2.9200048 5
2.9360049 5
2.952005 5
2.968005 5
2.984005 5
3.000005 5
3.016005 5
3.032005 5
3.048005 5
3.0640051 5
3.0800052 5
3.0960052 5
3.1120052 5
3.1280053 5
3.1440053 5
3.1600053 5
3.1760054 5
3.1920054 5
3.2080054 5
3.2240055 5
3.2400055 5
3.2560055 5
3.2720056 5
3.2880056 5
3.3040056 5
3.3200057 5
3.3360057 5
3.3520057 5
3.3680058 5
3.3840058 5
3.4000058 5
3.4160058 5
3.432006 5
3.448006 5
3.464006 5
3.480006 5
3.496006 5
3.512006 5
3.528006 5
3.544006 5
3.5600061 5
3.5760062 5
3.5920062 5
3.6080062 5
3.6240063 5
3.6400063 5
3.6560063 5
3.6720064 5
3.6880064 5
3.7040064 5
3.7200065 5
3.7360065 5
3.7520065 5
3.7680066 5
3.7840066 5
3.8000066 5
3.8160067 5
3.8320067 5
3.8480067 5
3.8640068 5
3.8800068 5
3.8960068 5
3.9120069 5
3.928007 5
3.944007 5
3.976007 5
4.008007 5
4.040007 5
4.072007 5
4.1040072 5
4.1360073 5
4.1680074 5
4.2000074 5
4.2320075 5
4.2640076 5
4.2960076 5
4.3280077 5
4.360008 5
4.392008 5
4.424008 5
4.456008 5
4.488008 5
4.520008 5
4.552008 5
4.584008 5
4.6160083 5
4.6480083 5
4.6800084 5
4.7120085 5
4.7440085 5
4.7760086 5
4.8080087 5
4.8400087 5
4.872009 5
4.904009 5
4.936009 5
4.968009 5
Text(0, 0.5, 'y [m]')
Plot: this plot is of the same variables, x and y position, as Problem 1 Part C of HW 1.
the red line is the graph of the spiral created by the discrete time dynamics
the blue cross markers are the individual plotted points generated by the discrete time dynamics to show their change in frequency
Observations: It is easy to see the difference in density of plotted points increases as the curvature of the spiral becomes more intense.
As the curvature of the spiral becomes more greater the time step (dt) decreases resulting in the ploints plotted closer together on the graph.
Thus, towardds the outside of the spiral the algorithm is more efficient, plotting less points (larger dt) but still capturing the behavior of the system, while towards the inside of the spiral the algortihm becomes less efficient (smaller dt) but captures more precise bahavior and small details as x and y change more in the same span of time.