# Sparse identification of nonlinear dynamical systems

###### tags: `system identification` `Newell-Franklin Car Following`

[Discovering governing equations from data by sparse identification of nonlinear dynamical systems](https://www.pnas.org/doi/10.1073/pnas.1517384113)

Consider dynamical systems:
$$
\frac{d}{dt}\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t))
$$

The vector $\mathbf{x}(t)\in\mathbb{R}^n$ denotes the state of a system at time t, and the function $\mathbf{f}$ represents the dynamic constraints that define the equations of motion of the system.

To determine the function $\mathbf{f}$ from data, we collect a time history of the state $\mathbf{x}(t)$ and either measure the derivative $\dot{\mathbf{x}}(t)$ or approximate it numerically from $\mathbf{x}(t)$.



:::info
pNEUMA dataset contains info:
[' lat', ' lon', ' speed', ' lon_acc', ' lat_acc', ' time']
- different vehicle trajectory has different time-duration, how does it impact matrix size
- only velocity magnitude
- matrix state would be changing because of new trajectories
:::

## Some Test Results
ref_ind = 0:Data Size $451\times246$:
X: 451 time steps (0.04 sec each) - 246 columns (83 vehicles $\times$ [Lat,Lon,Speed])
Y: 451 time steps (0.04 sec each) - 246 columns (83 vehicles $\times$ [Speed,Lon_acc,Lat_acc])
Train Test Split ratio = 0.7 --->training rows = 316
Sparsity Threshold = 0.0001

Method: Polynomial orders = 2: ---> 30628 columns
The training error is: 1.532031168463187e-12
The testing error is: 37.26821366653878

Method: Polynomial orders = 2 + trigonometric terms: ---> 31120 columns
The training error is: 9.681934018530355e-13
The testing error is: 33.07795254424301

Method: Polynomial orders = 1: ---> 247 columns
The training error is: 0.0003327737025303104
The testing error is: 18.34513192353467

Method: Polynomial orders = 1 + trigonometric terms: ---> 739 columns
The training error is: 7.730648522020209e-05
The testing error is: 28.31391796302265

ref_ind = 600: Data Size $2201\times48$:
Method: Polynomial orders = 1: ---> 49 columns
The training error is: 0.030748302679976674
The testing error is: 0.7877024058765171

Method: Polynomial orders = 1 + trigonometric terms: ---> 145 columns
The training error is: 0.023418644649387133
The testing error is: 2.9321339236273505

### Car following model

road length = 10km, total time = 500 s, time step = 0.018s, #vehicle = 50
Method: Polynomial orders = 2
The training error is: 0.16625532598612017
The testing error is: 2362843.5426082565

Method: Polynomial orders = 1
The training error is: 0.8637392497224802
The testing error is: 898.0514427068016

road length = 10km, total time = 1000 s, time step = 0.18s, #vehicle = 50, split_ratio = 0.7

Method: Polynomial orders = 2
X_train shape = (3888, 5151)
Y_train shape = (3888, 100)
The training error is: 0.03050586955449847
The testing error is: 742401.2507981361

Method: Polynomial orders = 1
X_train shape = (3888, 101)
Y_train shape = (3888, 100)
The training error is: 0.7308672634802567
The testing error is: 1.1689105160687008

Method: Polynomial orders = 1 + trigonometric terms
X_train shape = (3888, 303)
Y_train shape = (3888, 100)
The training error is: 0.594687679804249
The testing error is: 17.46102769927894

Newell-fanklin

# Implementing Newell car-following model
Newell’s car-following model entails a basic assumption: that the trajectory of following vehicle n is consistent with that of leading vehicle n - 1 in a homogeneous space. The linear correlation between the spacing headway and speed is described by
$$
\Delta x_n(t) = \tau_n v_n(t) + \delta_n
$$

where $\Delta x_n(t)$ is the spacing headway (in unit of m), $v_n(t)$ is the speed (in unit of m/s) of the vehicle n, $\tau_n$ and $\delta_n$ are referred to as the response time and critical jam spacing. The parameter ($\tau_n$ , $\delta_n$) of the NCM is assumed to be constant for each driver, independent of driving conditions, but to vary from vehicle to vehicle.

The NCM is described mathematically in the discrete-time car-following context by
$$
v_n(t+\tau_n) = \min{\{v_{\max},\frac{\Delta x_n(t)-\delta_n}{\tau_n}\}} \\
x_n(t+\Delta t) = x_n(t) + v_n(t+\Delta t)\Delta t
$$

where $v_{\max}$ is the maximum speed (in unit of m/s) under the free-flow condition, and $\Delta t$ is the numerical update time interval.

:::info
$$
\Delta t \text{ or } \tau_n = 1 s \\
v_{\max} = 120 \text{km/h for high way} \\
v_{\max} = 54 \text{km/h for city traffic}
$$
:::

:::info
Newell-Frankin speed-density relation estimates traffic densities from estimated speeds in each cell.
$$
\rho(x_i,t_i) = \frac{\rho_{jam}}{1-\frac{v_f}{\omega}\ln{(1-\frac{v(x_i,t_i)}{v_f})}}
$$
where $\rho_{jam}$ is the jam density (143 veh/km), $v_f$ is the free-flow speed, $\omega$ is the shock wave speed.
:::

:::warning
Different eqations for Newell model
$$
\tau = 1/(wk_j) \\
\delta = 1/k_j
$$
:::

# Implementing Heterogeneous car-following model with Newell-Franklin (stationary) speed-spacing relation

Newell-Franklin (stationary) speed-spacing relation:
$$
v_n(s) = v_{n,f} - v_{n,f}\exp{\frac{-c_n}{v_{n,f}}(s-d_n)}
$$

where $v_{n,f}$ is driver n’s desired (free-flow) speed, $d_n$ is minimum safety distance, $c_n$ is the inverse of the reaction time of driver n when their speed is restricted by the trajectory of their leader. The spacing between vehicle n and their leader, n − 1, by
$$
s_n(t) = x_{n-1}(t) -x_n(t)
$$

The position dynamics are given, for any n, by
$$
x_n(t+\Delta t) = x_n(t) + v_n(t)\Delta t
$$

The spacing dynamics can be simulated using
$$
s_n(t+\Delta t) = s_n(t) + (v_{n-1}(s_{n-1}(t))-v_n(s_n(t)))\Delta t
$$

:::info
$x_n(t) \to s_n(t) \to v_n(t) \to x_n(t+\Delta t) \to s_n(t+\Delta t) \to v_n(t+\Delta t)$
:::

$\mathbf{V} = \Theta(\mathbf{S})\Xi$

$$
\mathbf{V} = \begin{bmatrix}
v_1(t = 0) & ... & v_{50}(t=0) \\
\vdots & \ddots & \vdots \\
v_1(t = 800) & ... & v_{50}(t=800)
\end{bmatrix}
$$

$$
\mathbf{S} = \begin{bmatrix}
s_1(t = 0) & ... & s_{50}(t=0) \\
\vdots & \ddots & \vdots \\
s_1(t = 800) & ... & s_{50}(t=800)
\end{bmatrix}
$$

2nd order, Without force sparsity, V-S model:
The training error is: 2.95501085437511e-14
The testing error is: 2.716796706169583e-06
2nd order, With force sparsity, threshold = 10**(-6), V-S model:
The training error is: 2.9636902661293757e-14
The testing error is: 2.716796656790843e-06
['s27 s29', 's28 s30', 's27 s31', 's25 s31', 's26 s30', 's26 s28', 's29 s31', 's26 s32', 's29^2', 's25 s29', 's3 s49', 's28^2', 's0 s9',...]

:::warning
$$
\mathbf{V} \neq \dot{\mathbf{S}}
$$
:::

2nd order, Without force sparsity, V-X model:
The training error is: 0.022489543906529928
The testing error is: 0.017399802487414916
2nd order, With force sparsity, threshold = 10**(-6), V-X model:
The training error is: 0.02248954390653043
The testing error is: 0.0173998024873622

2nd order, With force sparsity, threshold = 10**(-4), V-X model:
The training error is: 0.022489570400765314
The testing error is: 0.017396874216536854

['1', 'x9', 'x0', 'x2', 'x3', 'x7', 'x8', 'x4', 'x1', 'x5', 'x6', 'x5 x9', 'x0 x5', 'x6 x9', 'x0 x9', 'x0 x6', 'x0 x4', 'x4 x9', 'x7 x9', 'x0 x7', 'x0^2', 'x2 x5', 'x2 x4', 'x4 x7', 'x0 x2',...]

3rd order, Without force sparsity, V-X model:
The training error is: 7.046342321560165e-13
The testing error is: 0.0003719990214855108
3rd order, With force sparsity, threshold = 10**(-6), V-X model:
The training error is: 2.3113518244171475e-12
The testing error is: 0.0003720309595927169
3rd order, With force sparsity, threshold = 10**(-4), V-X model:
The training error is: 0.00445702339073122
The testing error is: 0.004908670528638443

['x0 x9', 'x0 x4 x6', 'x1 x9', 'x4 x5 x7', 'x2 x9', 'x7 x9', 'x4 x8 x9', 'x0^2', 'x8 x9', 'x4 x6 x9', 'x0 x1', 'x4 x8^2', 'x0 x4 x7', 'x6^2', 'x2 x8', 'x2 x6 x9', 'x1 x8', 'x3 x4 x8', 'x3 x5', 'x3 x4 x9', 'x3 x8', 'x4 x6^2', 'x3 x9', 'x3^2', 'x0 x7 x8',...]

5th order, Without force sparsity, V-X model:
The training error is: 6.383434894184935e-14
The testing error is: 0.0008670534978616565
5th order, With force sparsity, threshold = 10**(-10), V-X model:
The training error is: 1.5395409276796716e-13
The testing error is: 0.000867115521032239
['x0^4 x9', 'x0^4 x2', 'x0^5', 'x0^3 x4 x9', 'x0^3 x9^2', 'x0^3 x2 x9', 'x0^3 x1 x2', 'x0^2 x1^2 x8', 'x0^3 x2^2', 'x0^3 x1 x9', 'x0 x1 x2^3', 'x0^2 x1 x9^2', 'x8 x9^4', 'x0^2 x2^3', 'x0^3 x2 x5', 'x0^2 x8^3', 'x0^2 x1^2 x3', 'x0^2 x1 x5 x9', 'x0^2 x4^2 x9', 'x0 x3^4', 'x3^3 x4 x9', 'x0 x1^2 x9^2', 'x4 x9^4',...]

2nd order, With force sparsity, LASSO alpha = 0.1, V-X model:
The training error is: 0.13213539016741258
The testing error is: 0.05916988952570645

['x1', 'x7', 'x5', 'x3', 'x9', 'x0 x1', 'x0^2', 'x1 x3', 'x7 x9', 'x5 x9', 'x2 x3', 'x8 x9', 'x0 x8', 'x1 x7', 'x1 x8', 'x2^2', 'x0 x3', 'x3 x7',...]


2nd order, With force sparsity, LASSO alpha = 1, V-X model:
The training error is: 0.13581649737376694
The testing error is: 0.056821847055094264
['x1^2', 'x0^2', 'x0 x9', 'x5 x9', 'x0 x1', 'x1 x8', 'x0 x8', 'x6^2', 'x1 x7', 'x3^2', 'x0 x3', 'x7^2', 'x7 x9', 'x6 x9', 'x5^2', 'x8^2', 'x8 x9', 'x0 x7', 'x1 x6', 'x2^2', 'x1 x3', 'x2 x3',...]


standardzation + threshold 10

2nd order, Without force sparsity, V-X model, standardization:
The training error is:
0.022489543906529918
The testing error is:
45242.07310310795
2nd order, With force sparsity, threshold = 10, V-X model, standardization:
The training error is:
0.9939046741040938
The testing error is:
30691.020734232574


standardzation + threshold 5

### related works
governing equations identifications

symbolic regression
[PySR](https://github.com/MilesCranmer/PySR)
[differentiable CGP (Cartesian Genetic Programming)](https://github.com/darioizzo/dcgp/)
[SISSO (Sure Independence Screening and Sparsifying Operator) method](https://github.com/Matgenix/pysisso)

[Physics-informed learning of governing equations from scarce data: ](https://www.nature.com/articles/s41467-021-26434-1)
calculate candidate terms/function from DNN, present a new PINN-SR paradigm to simultaneously model the system response and identify the parsimonious closed form of the governing PDE(s). T

## KMD experiments result
delay = 5 (Hankel)
print(LA.norm(KMD_estimates - temp,'fro')) = 1057.6711076812996
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.27987708017489304


print(LA.norm(KMD_estimates - temp,'fro')) = 717.0728858126741
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.18974921797174485

print(LA.norm(KMD_estimates - temp,'fro')) = 403.2049282361322
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.10669462104741505

delay > 17: overflow encountered in exp

Extended DMD: (order = 2)
print(LA.norm(KMD_estimates - temp,'fro')) = 708.9028372004465
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.2244964014149628

Plot the Koopman modes (delay=5):

Plot the Koopman modes (delay=17):

NGSIM:
print(LA.norm(KMD_estimates - temp,'fro')) = 2.267145951926545
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.0007290897529683445

Koopman mode extracted:




### prediction power
latency = 10
print(LA.norm(KMD_estimates - temp,'fro')) = 15955.912491682442
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 8.943063195497011

latency = 20
print(LA.norm(KMD_estimates - temp,'fro')) = 1002.4409039887216
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.561853943407846

latency = 30
print(LA.norm(KMD_estimates - temp,'fro')) = 1017.4028307035385
print(LA.norm(KMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 0.5702398916390582

print(LA.norm(DMD_estimates - temp,'fro')) = 1.186576308978214e-09
print(LA.norm(DMD_estimates - temp,'fro')/LA.norm(temp,'fro')) = 4.151391805194544e-15

960.5111081472776
0.5383528860925847

:::warning
lack of forcasting power of long-term forecasting/high resolution?
:::

### fourier & Koopman (& Vandermonde)
[From Fourier to Koopman: Spectral Methods for Long-term Time Series Prediction](https://jmlr.csail.mit.edu/papers/volume22/20-406/20-406.pdf)

[vandermonde decomposition](https://www.ethanepperly.com/index.php/2022/06/27/the-vandermonde-decomposition/)

F:expanding the signal in the basis of sines and cosines
K:nonlinear oscillatory basis
observable functions assumed to be invertible, can be approximated by a function (as opposed to a functional).

Fourier transform expands a signal into a possibly infinite sum of sines and cosines.
Koopman theory explains the signal as a nonlinear function of sines and cosines.


print(LA.norm(K_future - x_hat[:,1200:],'fro')) = 812.4414496627545
print(LA.norm(K_future - x_hat[:,1200:],'fro')/LA.norm(K_future,'fro')) = 0.47928736905877867

### parsimonious model
$$
\mathbf{x}_{kmd}(t) =\sum_{i=1}^l b_{0i}\mathbf{v}_i e^{\omega_i t}=\mathbf{V}\mathbf{e}^{\omega t}\mathbf{b}_0
$$

Find the prominent frequency?
$$
\arg\min \|\mathbf{e}^{\omega t}\|_0 \quad s.t. |\mathbf{x}(t)-\mathbf{V}\mathbf{e}^{\omega t}\mathbf{b}_0|=0 \quad \text{for all }t
$$

Sparsify the 'amplitude'?
$$
\arg\min \|\mathbf{V}\|_0 \quad s.t. |\mathbf{x}(t)-\mathbf{V}\mathbf{e}^{\omega t}\mathbf{b}_0|=0 \quad \text{for all }t
$$