Data-driven identification of vehicle dynamics using Koopman operator

# Data-driven identification of vehicle dynamics using Koopman operator

Vít Cibulka1, Tomáš Haniš2 and Martin Hromčík3 Department of Control Engineering, Faculty of Electrical Engineering, Czech Technical University in Prague
Email: 1cibulka.vit@fel.cvut.cz 2hanis.tomas@fel.cvut.cz 3hromcik.martin@fel.cvut.cz
###### 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 vehicle system.

Koopman operator, eigenfunctions, basis functions, data-driven methods, identification

## 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 .

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 . 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 , where is longitudinal velocity, lateral velocity, yawrate and 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.

 (1)
 Jzz¨ψ=4∑i=1riFi,y (2)
 JRi¨ρRi=Ma,Ri−Mb,Risign(˙ρRi)−rFRi,x,∀i=1,3 (3)

Where

 r=⎡⎢ ⎢ ⎢⎣r1r2r3r4⎤⎥ ⎥ ⎥⎦T=⎡⎢⎣⎡⎢⎣lv00⎤⎥⎦,⎡⎢⎣lv00⎤⎥⎦,⎡⎢⎣−lh00⎤⎥⎦,⎡⎢⎣−lh00⎤⎥⎦⎤⎥⎦ (4)

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. is the vehicle mass, is a force acting on i-th wheel along x/y axis in body-fixed coordinates. is a force acting along x axis in wheel coordinate system (direct output of the tire model). The term is an approximation of air-resistance, is a drag coefficient, is air density and is the total surface exposed to the air flow. is the vehicle inertia about z-axis. the wheel inertia about y-axis. Inputs are wheel torques (throttle), (break) and steering angles .

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

 ˙ρR1=˙ρR2 (5)
 ˙ρR3=˙ρR4. (6)

The forces are calculated using the “Pacejka magic formula” 

 F=Dcos(Carctan(Bx−E(Bx−arctan(Bx)))). (7)

The same formula can be used for calculating (tire longitudinal force) and (tire lateral force) with a different set of parameters for each. The argument can be either sideslip angle or longitudinal slip (see ) for calculating or respectively. The parameters B,C,D and E are usually time-dependent. This work uses the Pacejka tire model  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

 [Fi,xFi,y]=[cos(δi)−sin(δi)sin(δi)cos(δi)][FRi,xFRi,y] (8)

### Ii-a 3 state singletrack

To further simplify the model, equation (3) can be omitted, resulting in a model with only 3 states: . 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 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)  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  will be reviewed now. Let us assume uncontrolled discrete nonlinear dynamical system with state and dynamics

 x+=f(x)y=g(x) (9)

with being the current state and the next state.

The Koopman operator is defined as

 (Kψ)(x)=ψ(f(x)) (10)

for each basis function where is the space of basis functions.

is infinite-dimensional linear operator, which can be approximated by EDMD. The approximation is done by solving optimization problem

 minAK∑j=1||ψ(x+j)−Aψ(xj)||22 (11)

where . The states and are obtained by simulating the nonlinear system model, is the cardinality of the simulated dataset . Note that the states 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 and perform a one-step simulation to obtain the .

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

 z+=Az^y=Cz (12)

where and is the output estimate. The matrix is obtained in a similar fashion

 minCK∑j=1||yj−Cψ(xj)||22 (13)

For controlled system

 x+=f(x,u)y=g(x) (14)

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

 minA,BK∑j=1||ψ(x+j)−(Aψ(xj)+Buj)||22 (15)

while (13) stays the same. The Koopman operator is approximated both by and and acts as a one-step predictor. Note that the learning dataset now consists of . See  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 and is much more difficult when the tire nonlinearities are present.

#### Iii-A2 Data selection

The learning dataset was selected as follows. The velocities and were uniformly sampled in the interval , yawrate was uniformly sampled in interval . 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 with , 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 values in the dataset .

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 ). In our case, polynomial basis functions yielded much better results. A polynomial basis consists of monomials of the elements of the state vector.

 Bk={vax⋅vby⋅˙ψc|a,b,c∈{0,1,2,…k}} (16)

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

 RMSE=100√∑k||xkoop(kTs)−xreal(kTs)||22√||∑kxreal(kTs)||22 (17)

calculated over 3375 trajectories with the same distribution of initial conditions as the learning dataset . The length of each trajectory was 30 samples, with . The results can be seen in Fig.2. Fig. 2: Comparison of polynomial bases with different orders. It is clear from the figure that large set of basis functions does not imply better prediction error. The best results were for polynomial bases of orders 6 and 7 with RMSE=11.1%.

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

 ˙ψ(x)=Aψ(x), (18)

the more basis functions are in , the more derivatives have to be expressed by . 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. Fig. 3: This figure shows the comparison between the real nonlinear system and linear system approximated by the Koopman operator with polynomial basis of order 7. The is input random with uniform distribution in the interval [−104,104]. The initial conditions were chosen randomly as [vx,vy,˙ψ]=[20.8533,−0.7222,−4.3479]. The prediction is restarted after every 50 samples, the restart is visualized by segmenting the whole trajectory into pieces, each 50 samples long.

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  can be used. The approach consists of creating the eigenfunctions of the Koopman operator from data and using them as basis functions . The matrices and 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 .

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

 ddtx=f(x) (19)

Starting in , the system will get to state in time . The state vector will be transformed with a function , defined as

 ψ(xt)=ψ(xt)λ,g=eλtg(x0) (20)

for arbitrary eigenvalue and function . The time derivative of is

 ddtψ(xt)=ddteλtg(x0)=λeλtg(x0)=λψ(xt) (21)

We see that evolves linearly with the system (19). For we would get

 ddt⎡⎢ ⎢ ⎢ ⎢ ⎢⎣ψ1(xt)ψ2(xt)⋮ψNψ(xt)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣λ1λ2⋱λNψ⎤⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢⎣ψ1(xt)ψ2(xt)⋮ψNψ(xt)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦ (22)
 ddtz=Az (23)

for . Choosing (20) as basis functions immediately yields the diagonal 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 for the system (19) and simulate the system for a time with a sampling period . Then for each sampled data point , evaluate a set of eigenfunctions

 ψ(xkTs)λ,g=ekTsg(x0) (24)

for some , , where is the initial condition of the state .

The values of can be interpolated in order to approximate their values for states which are not in the learning dataset. The approximation will be denoted as

The set should be a non-recurrent set, meaning that the sampled trajectories of the system (19) should not return to the states from set . 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 , because the system naturally loses its kinetic energy as time progresses. For more information and proofs of the above stated facts, see .

The matrix can be obtained in a similar fashion as in the previous case. To approximate the output , solve

 minCM∑i=1||g(xi)−C^ψ(xi)||22 (25)

where is the total number of sampled states.

The controlled case is not considered here, for more information on the derivation of the matrix, see .

#### Results

The simulation experiments were performed with the model described in II-A. The non-recurrent set 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 , which is an energy equivalent of a car weighting , riding at . The set is depicted in Fig. 4. Fig. 4: Set of states with constant kinetic energy equivalent to 1300kg car riding at 100km/h. Each vertex of the mesh is one of the 441 initial states from which the learning dataset was constructed.

Next, the nonlinear system was simulated from 441 initial states from the set (the set is visualized in Fig.4) for sufficiently long time in order to cover the state space with enough data, in our case . The trajectories were sampled at sampling rate .

Then the eigenfunctions were calculated. The eigenvalues were chosen in a similar fashion as in . By applying Dynamic Mode Decomposition (DMD) algorithm  to the dataset, three eigenvalues were obtained. Additional eigenvalues were constructed as linear combinations of elements from , more in .

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

 g(x)=||x−xc||2log(||x−xc||) (26)

The centers were chosen randomly with normal distribution , where and are the mean and standard deviation of the whole dataset. Note that instead of interpolating the functions , 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 (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. Fig. 5: Each ball in the figure is an initial condition. Its size and color indicates prediction error of its associated trajectory. The mean of RMSE is 23%, standard deviation 15%. Fig. 6: Trajectories used for construction of eigenfunctions are depicted as gray lines. The prediction error is large in areas with very little data. The trajectories converge to planes, forming X-shape around origin. The shape is caused by the vehicle not having enough energy to spin and going into a drift while decreasing vy and ˙ψ. Fig. 7: This figure shows a trajectory with RMSE=23% which is equal to the mean RMSE of the whole testing dataset.

Notice that that in Fig.7, the initial conditions are different. This is caused by the approximation error of 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.

## References

•  D. Schramm, M. Hiller, and R. Bardini, Vehicle Dynamics.   Springer Berlin Heidelberg, 2014.
•  H. Pacejka, Tire and Vehicle Dynamics.   Elsevier LTD, Oxford, 2012. [Online]. Available: https://www.ebook.de/de/product/18341528/hans_pacejka_tire_and_vehicle_dynamics.html
•  M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” Journal of Nonlinear Science, vol. 25, no. 6, pp. 1307–1346, jun 2015.
•  M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” Automatica, vol. 93, pp. 149–160, jul 2018.
•  M. Korda and I. Mezić, “Learning koopman eigenfunctions for prediction and control: the transient case.”
•  P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” Journal of Fluid Mechanics, vol. 656, pp. 5–28, jul 2010.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters   