Data-driven identification of vehicle dynamics using Koopman operator

TL;DR

This paper explores vehicle dynamics identification using the Koopman operator, transforming nonlinear systems into linear ones in higher-dimensional spaces through basis functions, with two approaches for basis selection demonstrated.

Contribution

It introduces two novel methods for selecting basis functions in Koopman operator-based vehicle dynamics modeling, including a data-driven approach and one leveraging known nonlinear dynamics.

Findings

Both approaches successfully model vehicle dynamics.

Eigenfunction-based basis functions improve prediction accuracy.

Discussion on basis function selection enhances Koopman operator application.

Abstract

This paper presents the results of identification of vehicle dynamics using the Koopman operator. The basic idea is to transform the state space of a nonlinear system (a car in our case) to a higher-dimensional space, using so-called basis functions, where the system dynamics is linear. The selection of basis functions is crucial and there is no general approach on how to select them, this paper gives some discussion on this topic. Two distinct approaches for selecting the basis functions are presented. The first approach, based on Extended Dynamic Mode Decomposition, relies heavily on expert basis selection and is completely data-driven. The second approach utilizes the knowledge of the nonlinear dynamics, which is used to construct eigenfunctions of the Koopman operator which are known by definition to evolve linearly along the nonlinear system trajectory. The eigenfunctions are then used as basis functions for prediction. Each approach is presented with a numerical example and discussion on the feasibility of the approach for a nonlinear

Equations

$$m_v(\left[\begin{array}{c}{\dot{v}}_x\\ {\dot{v}}_y\end{array}\right]+\dot{\psi }\left[\begin{array}{c}-v_y\\ v_x\end{array}\right])=\sum _{i=1}^4\left[\begin{array}{c}{\mathbf{F}}_{i,x}\\ {\mathbf{F}}_{i,y}\end{array}\right]-\frac{1}{2}{c}_w\rho A_w\sqrt{v_x^2+v_y^2}\left[\begin{array}{c}{v}_x\\ v_y\end{array}\right]$$

$$J_{zz}\ddot{\psi }=\sum _{i=1}^4{\mathbf{r}}_i{\mathbf{F}}_{i,y}$$

$$J_{R_i}{\ddot{\rho }}_{R_i}=M_{a,R_i}-M_{b,R_i}sign({\dot{\rho }}_{R_i})-r{F}_{R_i,x},\forall i=1,3$$

$$\mathbf{r}={\left[\begin{array}{c}{\mathbf{r}}_1\\ {\mathbf{r}}_2\\ {\mathbf{r}}_3\\ {\mathbf{r}}_4\end{array}\right]}^T=\left[\begin{array}{c}\left[\begin{array}{c}{l}_v\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{l}_v\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}-l_h\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}-l_h\\ 0\\ 0\end{array}\right]\end{array}\right]$$

$${\dot{\rho }}_{R_1}={\dot{\rho }}_{R_2}$$

$${\dot{\rho }}_{R_3}={\dot{\rho }}_{R_4}.$$

$$F=D\cos (C\arctan (Bx-E(Bx-\arctan (Bx)))).$$

$$\left[\begin{array}{c}{F}_{i,x}\\ F_{i,y}\end{array}\right]=\left[\begin{array}{cc}\cos ({\delta }_i)& -\sin ({\delta }_i)\\ \sin ({\delta }_i)& \cos ({\delta }_i)\end{array}\right]\left[\begin{array}{c}{F}_{R_i,x}\\ F_{R_i,y}\end{array}\right]$$

$$\begin{array}{rl}{x}^{+}& =f(x)\\ y& =g(x)\end{array}$$

$$(\mathcal{K}\psi )(x)=\psi (f(x))$$

$$\underset{A}{\min }\sum _{j=1}^K||\mathbf{\psi }(x_j^{+})-A\mathbf{\psi }(x_j)|{|}_2^2$$

$$\begin{array}{rl}{z}^{+}& =Az\\ \hat{y}& =Cz\end{array}$$

$$\underset{C}{\min }\sum _{j=1}^K||y_j-C\mathbf{\psi }(x_j)|{|}_2^2$$

$$\begin{array}{rl}{x}^{+}& =f(x,u)\\ y& =g(x)\end{array}$$

$$\underset{A,B}{\min }\sum _{j=1}^K||\mathbf{\psi }(x_j^{+})-(A\mathbf{\psi }(x_j)+B{u}_j)|{|}_2^2$$

$$B_k=\{v_x^a\cdot v_y^b\cdot {\dot{\psi }}^c|a,b,c\in \{0,1,2,\dots k\}\}$$

$$RMSE=100\frac{\sqrt{{\sum }_k||x_{koop}(k{T}_s)-x_{real}(k{T}_s)|{|}_2^2}}{\sqrt{||{\sum }_k{x}_{real}(k{T}_s)|{|}_2^2}}$$

$$\dot{\mathbf{\psi }}(x)=A\mathbf{\psi }(x),$$

$$\frac{d}{dt}x=f(x)$$

$$\psi (x_t)=\psi (x_t{)}_{\lambda ,g}=e^{\lambda t}g(x_0)$$

$$\frac{d}{dt}\psi (x_t)=\frac{d}{dt}{e}^{\lambda t}g(x_0)=\lambda e^{\lambda t}g(x_0)=\lambda \psi (x_t)$$

$$\frac{d}{dt}\left[\begin{array}{c}{\psi }_1(x_t)\\ {\psi }_2(x_t)\\ ⋮\\ {\psi }_{N_{\psi }}(x_t)\end{array}\right]=\left[\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_{N_{\psi }}\end{array}\right]\left[\begin{array}{c}{\psi }_1(x_t)\\ {\psi }_2(x_t)\\ ⋮\\ {\psi }_{N_{\psi }}(x_t)\end{array}\right]$$

$$\frac{d}{dt}z=Az$$

$$\psi (x_{k{T}_s}{)}_{\lambda ,g}=e^{k{T}_s}g(x_0)$$

$$\underset{C}{\min }\sum _{i=1}^M||g(x_i)-C\hat{\mathbf{\psi }}(x_i)|{|}_2^2$$

$$g(x)=||x-x_c|{|}^2\log (||x-x_c||)$$

Full text

I Introduction

Car is a highly nonlinear dynamical system which has been used for decades without any sort of sophisticated automatic control system. Our goal is to create full-time-full-authority control system, which means that the driver will rather set the reference for the car, and the ultimate control over the car inputs (steering angle, motor torques) will lie with the control system, protecting the driver from himself. However, control of nonlinear systems is not an easy task.

The Koopman operator is becoming an increasingly popular tool for nonlinear control. The Koopman operator is a linear operator which can in theory approximate any nonlinear system. The idea lies in lifting the state space of the nonlinear system to a higher-dimensional space by means of a nonlinear transformation. The lifted system should then be evolving linearly along with the original nonlinear system. The lifted state space is then transformed back with a linear transformation to the original state space and the nonlinear states (or system outputs) are recovered.

Such an approach is suitable for any linear control-design method such as LQG, MPC and $H_{\infty}$.

In this paper, we present the results of application of the Koopman operator framework to vehicle system in an attempt to create a linear representation of a vehicle and use it to control the vehicle optimally with linear MPC.

II Singletrack model

The singletrack model is depicted in Fig.1. The model was derived from 16-state twin-track model described in [1]. The states and parameters of this model are a subset of states and parameters of the full twin-track model.

State vector of the singletrack model is $\begin{bmatrix}v_{x}(m/s)&v_{y}(m/s)&\dot{\psi}(rad/s)&\dot{\rho}_{f}(rad/s)&\dot{\rho}_{r}(rad/s)\end{bmatrix}^{T}$, where $v_{x}$ is longitudinal velocity, $v_{y}$ lateral velocity, $\dot{\psi}$ yawrate and $\dot{\rho}_{f/r}$ are front/rear wheel angular rates.

This model has 4 wheels, with two wheel always being in the same place. This allows for usage of asymmetric tire models, reduces the work needed to transition to twin-track model and is less error-prone than the standard approach (with only two tires, the user has to remember that each tire should generate twice as much force).

The vehicle body is modeled as a rigid body using Newton-Euler equations.

[TABLE]

[TABLE]

[TABLE]

Where

[TABLE]

is the vector describing position of each wheel with respect to the center of gravity. The wheels are numbered in this order: front-left, front-right, rear-left, rear-right. $m_{v}$ is the vehicle mass, $\mathbf{F}_{i,x/y}$ is a force acting on i-th wheel along x/y axis in body-fixed coordinates. $F_{R_{i},x}$ is a force acting along x axis in wheel coordinate system (direct output of the tire model). The term $-\frac{1}{2}c_{w}\rho A_{w}\sqrt{v_{x}^{2}+v_{y}^{2}}\begin{bmatrix}v_{x}\\ v_{y}\end{bmatrix}$ is an approximation of air-resistance, $c_{w}$ is a drag coefficient, $\rho$ is air density and $A_{w}$ is the total surface exposed to the air flow. $J_{zz}$ is the vehicle inertia about z-axis. $J_{R_{i}}$ the wheel inertia about y-axis. Inputs are wheel torques $M_{a,R_{i}}$ (throttle), $M_{b,R_{i}}$ (break) and steering angles $\delta_{f/r}$.

Since this is a 4-wheel singletrack model, the following holds:

[TABLE]

[TABLE]

The forces $\begin{bmatrix}F_{R_{i},x}\\ F_{R_{i},y}\\ \end{bmatrix}$ are calculated using the “Pacejka magic formula” [2]

[TABLE]

The same formula can be used for calculating $F_{x}$ (tire longitudinal force) and $F_{y}$ (tire lateral force) with a different set of parameters for each. The argument $x$ can be either sideslip angle $\alpha$ or longitudinal slip $\kappa$ (see [2]) for calculating $F_{y}$ or $F_{x}$ respectively. The parameters B,C,D and E are usually time-dependent. This work uses the Pacejka tire model [2] with coefficients from the Automotive challenge 2018 organized by Rimac Automobili.

The transformation of tire forces from wheel-coordinate system to car coordinate system is done as follows

[TABLE]

II-A 3 state singletrack

To further simplify the model, equation (3) can be omitted, resulting in a model with only 3 states: $\begin{bmatrix}v_{x}&v_{y}&\dot{\psi}\end{bmatrix}^{T}$. The inputs are then longitudinal slips (which were previously derived from (3)) and steering angles. This model will be used for the Eigenfunction approach in III-B.

II-B 3 state singletrack without tire model

To simplify the model even more, one can omit the tire model and use the tire forces $\begin{bmatrix}F_{i,x}\\ F_{i,y}\\ \end{bmatrix}$ as input, assuming the existence of a higher level control system controlling the tire forces and thus securing the assumption that the car can be controlled directly by force reference.

This model will be used for the Extended Dynamic Mode Decomposition (EDMD) [3] approach in III-A.

III The Koopman Operator

III-A Extended Dynamic Mode Decomposition approach

III-A1 Framework description

The Koopman operator framework for controlled systems, as described in [4] will be reviewed now. Let us assume uncontrolled discrete nonlinear dynamical system with state $x$ and dynamics $f(\cdot)$

[TABLE]

with $x$ being the current state and $x^{+}$ the next state.

The Koopman operator $\mathcal{K}:\mathcal{F}\rightarrow\mathcal{F}$ is defined as

[TABLE]

for each basis function $\psi:\mathbb{R}^{N}\rightarrow\mathbb{R}$ where $\mathcal{F}$ is the space of basis functions.

$\mathcal{K}$ is infinite-dimensional linear operator, which can be approximated by EDMD. The approximation is done by solving optimization problem

[TABLE]

where $\boldsymbol{\psi}=\begin{bmatrix}\psi_{1}(x)&\psi_{2}(x)&\ldots\psi_{N_{\psi}}(x)\end{bmatrix}^{T}$. The states $x_{j}$ and $x_{j}^{+}$ are obtained by simulating the nonlinear system model, $K$ is the cardinality of the simulated dataset $\mathcal{D}=\begin{bmatrix}x_{j}&x_{j}^{+}&y_{j}\end{bmatrix}_{j=1}^{K}$. Note that the states $x_{j}$ do not have come from some trajectory of the system. When the system model is available, it is sufficient to sample the state space with $x_{j}$ and perform a one-step simulation to obtain the $x_{j}^{+}$.

The matrix $A$ now defines a discrete linear system which approximates the nonlinear dynamics from (9)

[TABLE]

where $z=\boldsymbol{\psi}(x)$ and $\hat{y}$ is the output estimate. The matrix $C$ is obtained in a similar fashion

[TABLE]

For controlled system

[TABLE]

the approach is similar. The optimization problem (11) changes to

[TABLE]

while (13) stays the same. The Koopman operator $\mathcal{K}$ is approximated both by $A$ and $B$ and acts as a one-step predictor. Note that the learning dataset now consists of $\mathcal{D}=\begin{bmatrix}x_{j}&x_{j}^{+}&y_{j}&u_{j}\end{bmatrix}_{j=1}^{K}$. See [4] for more information.

This approach will be showed on the model which does not include the tire nonlinearities, described in II-B, because the selection of $\boldsymbol{\psi}$ and $\mathcal{D}$ is much more difficult when the tire nonlinearities are present.

III-A2 Data selection

The learning dataset was selected as follows. The velocities $v_{x}$ and $v_{y}$ were uniformly sampled in the interval $\begin{bmatrix}-30&30\end{bmatrix}m/s$, yawrate was uniformly sampled in interval $\begin{bmatrix}-10&10\end{bmatrix}rad/s$. The number of samples was 15 for each state, resulting in 3 1-by-15 vectors.

The input forces to the model were sampled uniformly in the interval $\begin{bmatrix}-F_{m}&F_{m}\end{bmatrix}N$ with $F_{m}\in\{1,5,10,100\}$, each interval with 15 samples, resulting in 4 4-by-15 vectors. The columns of sampled vectors were combined into all possible combinations, resulting in $4\cdotp 15^{4}$ values in the dataset $\mathcal{D}$.

Note that the tire force of a regular car is usually in order of thousands of Newtons. However, the results were much better with relatively small inputs forces. Large forces were too dominant and their influence on the system overshadowed the system dynamics and since there are no tire nonlinearities, the system can be sufficiently excited with small input forces.

III-A3 Basis selection

There is unfortunately no general recommendation to determine which basis functions to choose for a given system. Some systems show satisfactory results with thin plate spline radial basis functions (see [5]). In our case, polynomial basis functions yielded much better results. A polynomial basis consists of monomials of the elements of the state vector.

[TABLE]

where $k$ is the order of the basis $B_{k}$. Bases with various orders were compared, the comparison was made using averaged root mean square error (RMSE)

[TABLE]

calculated over 3375 trajectories with the same distribution of initial conditions as the learning dataset $\mathcal{D}$. The length of each trajectory was 30 samples, with $T_{s}=0.01$. The results can be seen in Fig.2.

One might think the error would decrease with adding more basis functions but error grows exponentially as the order increases, which means that one cannot simply add more basis functions and expect the error to diminish. When considering linear Koopman system

[TABLE]

the more basis functions are in $\boldsymbol{\psi}$, the more derivatives $\boldsymbol{\dot{\psi}}$ have to be expressed by $\boldsymbol{\psi}$. The derivatives of the polynomial basis grow in complexity with increasing order, which makes it harder to express them, so the result makes sense, although it might seem counter-intuitive at first glance.

Results

In our experiments, a polynomial basis with order 7 was used. The Koopman operator was approximated using the dataset described above. The results can be seen in Fig. 3.

The yawrate is tracked without error because it is simply an integral of the input force as seen in (2).

III-B Basis functions from data

To deal with the problems of basis function selection, the approach from [5] can be used. The approach consists of creating the eigenfunctions of the Koopman operator from data and using them as basis functions $\boldsymbol{\psi}$. The matrices $A$ and $B$ are calculated separately, contrary to the previous approach, which allows for optimizing the resulting system for a multiple-step prediction as opposed to one-step prediction of the previous approach. More in [5].

The basic idea of an eigenfunction will now be established. Let us consider uncontrolled nonlinear continuous system

[TABLE]

Starting in $x_{0}$, the system will get to state $x_{t}$ in time $t$. The state vector $x_{t}$ will be transformed with a function $\psi(x_{t})$, defined as

[TABLE]

for arbitrary eigenvalue $\lambda$ and function $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$. The time derivative of $\psi(x_{t})$ is

[TABLE]

We see that $\psi(x_{t})$ evolves linearly with the system (19). For $\boldsymbol{\psi}(x_{t})=\begin{bmatrix}\psi_{1}(x_{t})&\psi_{2}(x_{t})&\ldots\psi_{N_{\psi}}(x_{t})\end{bmatrix}^{T}$ we would get

[TABLE]

[TABLE]

for $z=\boldsymbol{\psi}(x)$. Choosing (20) as basis functions immediately yields the diagonal $A$ matrix. Notice that in this case, the Koopman system is derived as continuous system, contrary to the discrete-time derivation in III-A. This is done simply because of the continuous-time definition of (20), the system can be discretized after the whole procedure.

The idea is to select a set of initial conditions $\Gamma$ for the system (19) and simulate the system for a time $T$ with a sampling period $T_{s}$. Then for each sampled data point $x_{kT_{s}}$, evaluate a set of eigenfunctions

[TABLE]

for some $\lambda\in\Lambda$, $g\in\mathcal{G}$, where $x_{0}\in\Gamma$ is the initial condition of the state $x_{kT_{s}}$.

The values of $\psi(x)_{\lambda,g}$ can be interpolated in order to approximate their values for states $x$ which are not in the learning dataset. The approximation will be denoted as $\hat{\psi}(x)_{\lambda,g}$

The set $\Gamma$ should be a non-recurrent set, meaning that the sampled trajectories of the system (19) should not return to the states from set $\Gamma$. This condition is not difficult to fulfill with a dissipative system such a car. If, for example, the set consisted of states with the same kinetic energy, none of the trajectories would return to the set $\Gamma$, because the system naturally loses its kinetic energy as time progresses. For more information and proofs of the above stated facts, see [5].

The matrix $C$ can be obtained in a similar fashion as in the previous case. To approximate the output $y=g(x)$, solve

[TABLE]

where $M$ is the total number of sampled states.

The controlled case is not considered here, for more information on the derivation of the $B$ matrix, see [5].

Results

The simulation experiments were performed with the model described in II-A. The non-recurrent set $\Gamma$ was chosen as set of states with the same kinetic energy, the reasoning for this selection was mentioned above.

The level of the kinetic energy was set to $500kJ$, which is an energy equivalent of a car weighting $1300kg$, riding at $100km/h$. The set is depicted in Fig. 4.

Next, the nonlinear system was simulated from 441 initial states from the set $\Gamma$ (the set $\Gamma$ is visualized in Fig.4) for sufficiently long time $T$ in order to cover the state space with enough data, in our case $T=1s$. The trajectories were sampled at sampling rate $T_{s}=0.01s$.

Then the eigenfunctions $\psi_{\lambda,g}$ were calculated. The eigenvalues were chosen in a similar fashion as in [5]. By applying Dynamic Mode Decomposition (DMD) algorithm [6] to the dataset, three eigenvalues $\Lambda_{DMD}$ were obtained. Additional eigenvalues were constructed as linear combinations of elements from $\Lambda_{DMD}$, more in [5].

For the functions $g$, thin plate spline radial basis functions were used:

[TABLE]

The centers $x_{c}$ were chosen randomly with normal distribution $N(\mu,\sigma^{2})$, where $\mu$ and $\sigma$ are the mean and standard deviation of the whole dataset. Note that instead of interpolating the functions $\psi_{\lambda,g}$, nearest neighbour was chosen. This was done for simplicity and it yielded good results.

The approximated system was evaluated using 2000 randomly selected initial conditions, all within the ellipsoid in Fig. 4. Both the Koopman system and the real one were then simulated for $0.5s$ (which is unnecessarily long for MPC control, but it demonstrates the capabilities of the approach). The error of each trajectory is again calculated using RMSE, defined in (17).

Results can be seen in Fig.5. There are two clusters of initial conditions with large prediction error. These are the areas with very little data as can be seen in Fig.6. Trajectory with RMSE equal to the mean RMSE can be seen in Fig.7.

Notice that that in Fig.7, the initial conditions are different. This is caused by the approximation error of $\boldsymbol{\hat{\psi}}$ which was minimized in (25).

IV Conclusion

This paper presented two algorithms based on the Koopman operator for identifying vehicle dynamics.

The EDMD approach in III-A showed good result, but only with the model II-B which does not include the tire nonlinearities which are the greatest challenge in vehicle dynamics.

The eigenfunction approach in III-B yielded very promising results on a model with the tire nonlinearities II-A. This approach will be the focus of our future work including the controlled case, which wasn’t discussed in this work, and its implementation with linear MPC.

Acknowledgment

This research was supported by the Czech Science Foundation (GACR) under contract No. 19-16772S. We would like to thank Milan Korda for valuable discussions on the topic of eigenfunctions and Koopman theory.