# Data-driven identification and reconstruction of partially observed dynamical systems

Pierre Tandeo

► To cite this version:

Pierre Tandeo. Data-driven identification and reconstruction of partially observed dynamical systems.
ISDA 2024: 10th International Symposium on Data Assimilation, RIKEN, Oct 2024, Kobe, Japan.
hal-04764759

HAL Id: hal-04764759

https://imt-atlantique.hal.science/hal-04764759v1

Submitted on 8 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Data-driven identification and reconstruction of
partially observed dynamical systems

Dr. Pierre Tandeo, associate professor

IMT Atlantique, Brest, France
RIKEN Cluster for Pioneering Research, Kobe, Japan

The 10th International Symposium on Data Assimilation (ISDA)
October 21-25, 2024, Kobe, Japan

<figure><img src="image_2.png" /></figure>

In collaboration with:

- Pierre Ailliot (Univ. Brest, France)

- Takemasa Miyoshi (RIKEN, Japan)

- Florian Sévellec (CNRS, France)

Many thanks to my colleagues at RIKEN
## Context:

* rapidly growing field of research

* still lot of developments to do

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Cheng et al., Machine Learning with Data Assimilation and Uncertainty Quantification for Dynamical Systems: a Review, IEEE/CAA Journal of Automatica Sinica, 2023</footnote>
Data assimilation formulation:

$$x(t + dt) = M(x(t), \eta(t))$$

$$y(t) = H(x(t), \epsilon(t))$$

Machine learning appears everywhere:

* state **x** and observations **y**

* operators $M$ and $H$

* errors $\eta$ and $\epsilon$

* inversion of the system $p(x|y)$
$$x(t + dt) = M(x(t), \eta(t))$$

$$y(t) = H(x(t), \epsilon(t))$$

State of the system (x):

* reduced-order modeling (principal components, VAE)

* state augmentation (parameters, latent variables)

<figure><img src="image_4.png" /></figure>

Source: Cheng et al. 2023, IEEE/CAA Journal of Automatica Sinica
$$x(t + dt) = \mathcal{M}(x(t), \eta(t))$$

$$y(t) = \mathcal{H}(x(t), \epsilon(t))$$

Dynamical operator (M):

* surrogate models (AI, RNN, analogs, AR1 processes)

* hybrid models (GCM core + ML parameterization)

<figure><img src="image_5.png" /></figure>

Source: Kochkov et al. 2024, Nature
$$ 
\begin{aligned}
\mathbf{x}(t + dt) &= \mathcal{M}(\mathbf{x}(t), \eta(t)) \\
\mathbf{y}(t) &= \mathcal{H}(\mathbf{x}(t), \epsilon(t))
\end{aligned}
$$

Model and observation errors ($\eta$ and $\epsilon$):

* innovation statistics, inflation methods, statistical approaches, NN

* joint estimation of error covariances **Q** and **R**

<figure><img src="image_6.png" /></figure>

Source: Tandeo et al. 2020, Monthly Weather Review
$$x(t + dt) = M(x(t), \eta(t))$$

$$y(t) = H(x(t), \epsilon(t))$$

Observation operator (H):

* emulation of transfer functions using ML (e.g., radiative)

* estimation of representativity error

<figure><img src="image_7.png" /></figure>

Source: Janjić et al. 2018, QJRMS
$$x(t + dt) = \mathcal{M}(x(t), \eta(t))$$

$$y(t) = \mathcal{H}(x(t), \epsilon(t))$$

Noisy and partial observations (y):

* aggregation, fusion, and data exclusion

* optimal sampling of the observations

<figure><img src="image_8.png" /></figure>

Source: Marcille et al. 2023, Wind Energy Science
$$ 
\begin{aligned}
\mathbf{x}(t+dt) &= \mathcal{M}(\mathbf{x}(t), \eta(t)) \\
\mathbf{y}(t) &= \mathcal{H}(\mathbf{x}(t), \epsilon(t))
\end{aligned}
$$

Full system resolution $p(\mathbf{x}|\mathbf{y})$:

* optimization problem (link between 4D-Var and NN)

* probabilistic ML using Bayesian approaches

<figure><img src="image_9.png" /></figure>

Source: Cheng et al. 2023, IEEE/CAA Journal of Automatica Sinica
**General framework:**

$$ \mathbf{x}(t + dt) = \mathcal{M}(\mathbf{x}(t), \eta(t)) $$

$$ \mathbf{y}(t) = \mathcal{H}(\mathbf{x}(t), \epsilon(t)) $$

**Proposed framework:**

$$ \mathbf{x}(t + 1) = \mathbf{M}\mathbf{x}(t) + \mathbf{Q}^{1/2}\eta(t) $$

$$ \mathbf{y}(t) = \mathbf{H}\mathbf{x}(t) + \mathbf{R}^{1/2}\epsilon(t) $$

**Assumptions:**

* **linear** model and observation operators

* **Gaussian** state, observations, and error terms

* **constant** operators and error covariances

<footnote n="1">¹Tandeo et al., *Data-driven Reconstruction of Partially Observed Dynamical Systems*, *Nonlinear Processes in Geophysics*, 2023</footnote>
**General framework:**

$$ \mathbf{x}(t + dt) = \mathcal{M}(\mathbf{x}(t), \eta(t)) $$

$$ \mathbf{y}(t) = \mathcal{H}(\mathbf{x}(t), \epsilon(t)) $$

**Proposed framework:**

$$ \mathbf{x}(t + 1) = \mathbf{M}\mathbf{x}(t) + \mathbf{Q}^{1/2}\eta(t) $$

$$ \mathbf{y}(t) = \mathbf{H}\mathbf{x}(t) + \mathbf{R}^{1/2}\epsilon(t) $$

**Goal:**

* make **probabilistic forecasts**

* purely **data-driven** approach

* interpretation of the retrieved model

<footnote n="1">¹Tandeo et al., *Data-driven Reconstruction of Partially Observed Dynamical Systems*, *Nonlinear Processes in Geophysics*, 2023</footnote>
<figure><img src="image_10.png" /></figure>

State-space model:
Estimate $\mathbf{x}(0), \dots, \mathbf{x}(T)$

Variables and matrices:
$\mathbf{x} = [x_1, x_2]$
<figure><img src="image_11.png" /></figure>

State-space model:
$$ \mathbf{x}(t+1) = \mathbf{M}\mathbf{x}(t) $$

Variables and matrices:
$$ \mathbf{M} = [\cdot \cdot] $$
<figure><img src="image_12.png" /></figure>

State-space model:

$$x(t+1) = Mx(t) + Q^{1/2}\eta(t)$$

Variables and matrices:

$$\eta(t) \sim N(0, I), \quad Q = \begin{bmatrix} \cdot & \cdot \\ \cdot & \cdot \end{bmatrix}$$
Noisy observations of the state

<figure><img src="image_13.png" /></figure>

State-space model:

$$ \mathbf{x}(t+1) = \mathbf{M}\mathbf{x}(t) + \mathbf{Q}^{1/2}\mathbf{\eta}(t) $$

$$ \mathbf{y}(t) = \mathbf{H}\mathbf{x}(t) $$

Variables and matrices:

$$ \mathbf{y} = [y_1, y_2], \quad \mathbf{H} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} $$
<figure><img src="image_14.png" /></figure>

State-space model:

$$ \mathbf{x}(t+1) = \mathbf{M}\mathbf{x}(t) + \mathbf{Q}^{1/2}\eta(t) $$

$$ \mathbf{y}(t) = \mathbf{H}\mathbf{x}(t) + \mathbf{R}^{1/2}\epsilon(t) $$

Variables and matrices:

$$ \epsilon(t) \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), \quad \mathbf{R} = \begin{bmatrix} \cdot \\ \cdot \\ \cdot \end{bmatrix} $$
<figure><img src="image_15.png" /></figure>

State-space model:

$$x(t + 1) = Mx(t) + Q^{1/2}\eta(t)$$

$$y(t) = Hx(t) + R^{1/2}\epsilon(t)$$

⇒ Kalman smoother

Variables and matrices:

$$\tilde{M}, \tilde{Q}, \tilde{R} = \underset{\mathbf{M}, \mathbf{Q}, \mathbf{R}}{\operatorname{argmax}} \mathcal{L}(\mathbf{y})$$

⇒ EM algorithm
# The Expectation-Maximization algorithm<fnref n="1" />

Initialization: M, Q, R

<footnote n="1">¹Shumway and Stoffer, *An approach to time series smoothing and forecasting using the EM algorithm*, Journal of Time Series Analysis, 1982</footnote>
**Initialization:** **M**, **Q**, **R**

**Expectation (DA step):**

* fix **M**, **Q**, **R**

* estimation of $\mathbf{x}_{0:T}^s$ and $\mathbf{P}_{0:T}^s$

* $\Rightarrow$ Kalman smoother

<footnote n="1">¹Shumway and Stoffer, *An approach to time series smoothing and forecasting using the EM algorithm*, Journal of Time Series Analysis, 1982</footnote>
**Initialization:** **M**, **Q**, **R**

**Expectation (DA step):**

* fix **M**, **Q**, **R**

* estimation of $\mathbf{x}_{0:T}^s$ and $\mathbf{P}_{0:T}^s$

* $\Rightarrow$ Kalman smoother

**Maximization (ML step):**

* fix $\mathbf{x}_{0:T}^s$ and $\mathbf{P}_{0:T}^s$

* estimation of **M**, **Q**, **R**

* $\Rightarrow$ analytic solution

<footnote n="1">¹Shumway and Stoffer, *An approach to time series smoothing and forecasting using the EM algorithm*, Journal of Time Series Analysis, 1982</footnote>
# The Expectation-Maximization algorithm<fnref n="1" />

**Initialization:** **M**, **Q**, **R**

**Expectation (DA step):**

* fix **M**, **Q**, **R**

* estimation of $\mathbf{x}_{0:T}^s$ and $\mathbf{P}_{0:T}^s$

* $\Rightarrow$ Kalman smoother

**Maximization (ML step):**

* fix $\mathbf{x}_{0:T}^s$ and $\mathbf{P}_{0:T}^s$

* estimation of **M**, **Q**, **R**

* $\Rightarrow$ analytic solution

Repeat E- and M-steps until convergence of $\mathcal{L}(\mathbf{y})$

<figure><img src="image_16.png" /></figure>

<footnote n="1">¹Shumway and Stoffer, *An approach to time series smoothing and forecasting using the EM algorithm*, Journal of Time Series Analysis, 1982</footnote>
$$ \mathcal{L}(\mathbf{y}) = \prod_{t=0}^{T} \mathcal{N} \left( 0; \mathbf{y}(t) - \mathbf{H}\mathbf{x}^f(t), \mathbf{H}\mathbf{P}^f(t)\mathbf{H}^\top + \mathbf{R} \right) $$
$$ \mathcal{L}(\mathbf{y}) = \prod_{t=0}^{T} \mathcal{N} \left( 0; \mathbf{y}(t) - \mathbf{H}\mathbf{x}^f(t), \mathbf{H}\mathbf{P}^f(t)\mathbf{H}^{\top} + \mathbf{R} \right) $$

Graphical representation of the innovation likelihood:

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
$$ \mathcal{L}(y) = \prod_{t=0}^{T} \mathcal{N} \left( 0; y(t) - \mathbf{H}\mathbf{x}^f(t), \mathbf{H}\mathbf{P}^f(t)\mathbf{H}^{\top} + \mathbf{R} \right) $$

Graphical representation of the innovation likelihood:

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

$\mathcal{L}(y)$ depends on $\mathbf{M}, \mathbf{Q}, \mathbf{R}$:

* $\mathbf{x}^f(t) = \mathbf{M}\mathbf{x}^a(t-1) \rightarrow$ forecast precision

* $\mathbf{P}^f(t) = \mathbf{M}\mathbf{P}^a(t-1)\mathbf{M}^{\top} + \mathbf{Q} \rightarrow$ forecast reliability

* $\mathbf{R} \rightarrow$ observation error
Introduction of a latent variable

<figure><img src="image_23.png" /></figure>

State-space model:

$$x(t + 1) = Mx(t) + Q^{1/2}\eta(t)$$

$$y(t) = Hx(t) + R^{1/2}\epsilon(t)$$

$$H = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}$$

Variables and matrices:

$$\text{Augmented state } \mathbf{x} = [x_1, x_2, z_1]$$

M and Q are 3 × 3

R is 2 × 2
<figure><img src="image_24.png" /></figure>

State-space model:

$$x(t + 1) = Mx(t) + Q^{1/2}\eta(t)$$

$$y(t) = Hx(t) + R^{1/2}\epsilon(t)$$

$$H = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}$$

Variables and matrices:

$$\text{Estimate } \hat{M}, \hat{Q}, \hat{R}$$

$$\text{Estimate } z_1(0), \dots, z_1(T)$$
<figure><img src="image_25.png" /></figure>

**State-space model:**

$$x(t + 1) = Mx(t) + Q^{1/2}\eta(t)$$

$$y(t) = Hx(t) + R^{1/2}\epsilon(t)$$

$$H = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}$$

**Variables and matrices:**

New state $x = [x_1, x_2, z_1, z_2]$

$M$ and $Q$ are $4 \times 4$

$R$ is $2 \times 2$
<figure><img src="image_26.png" /></figure>

**State-space model:**

$$x(t + 1) = Mx(t) + Q^{1/2}\eta(t)$$

$$y(t) = Hx(t) + R^{1/2}\epsilon(t)$$

$$H = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}$$

**Variables and matrices:**

$$\text{Estimate } \hat{M}, \hat{Q}, \hat{R}$$

$$\text{Estimate } z_2(0), \dots, z_2(T)$$
## Latent variables:

* how many are needed?

* what is there meaning?

* why do we need them?

<figure><img src="image_27.png" /></figure>

## Quality of the forecasts:

* are they precise?

* are they reliable?

* how they compare?

<figure><img src="image_28.png" /></figure>
Application to Lorenz-63: formulation and OSSE

<figure><img src="image_29.png" /></figure>

Observations, state, and latent variables:

• $\mathbf{y} = [y_2, y_3] \rightarrow \textbf{partial observations}$ of the system

• $\mathbf{x} = [x_2, x_3, z] \rightarrow \textbf{augmented state}$ with latent variables z

• $\mathbf{z} = [z_1, \dots, z_k] \rightarrow \textbf{latent variables}$

OSSE:

• $dt = 0.001 \rightarrow \text{observation sampling}$

• $\text{Var}(\epsilon) = 1 \rightarrow \text{observation error variance}$
Estimation of $\hat{M}, \hat{Q}, \hat{R}$ and $z_1, \dots, z_k$ using training data

<figure><img src="image_30.png" /></figure>

Time series are divided in three parts:

- train, test, and validation datasets

- EM algorithm is applied to the training data
# Estimation and meaning of the latent variables

## Estimation of $z_1$ and $z_2$:

<figure><img src="image_31.png" /></figure>
# Estimation and meaning of the latent variables

## Estimation of $z_1$ and $z_2$:

<figure><img src="image_32.png" /></figure>

## Meaning of $z_1$:

$$z_1 = a_2 \dot{x}_2 + a_3 \dot{x}_3$$
# Estimation and meaning of the latent variables

## Estimation of $z_1$ and $z_2$:

<figure><img src="image_33.png" /></figure>

### Meaning of $z_1$:

$$z_1 = a_2 \dot{x}_2 + a_3 \dot{x}_3$$

### Meaning of $z_2$:

$$ \begin{aligned} z_2 &= b_1 \dot{z}_1 + b_2 \dot{x}_2 + b_3 \dot{x}_3 \\ &= b_2 \dot{x}_2 + b_3 \dot{x}_3 + b_1 a_2 \ddot{x}_2 + b_1 a_3 \ddot{x}_3 \end{aligned} $$
# Estimation and meaning of the latent variables

## Estimation of $z_1$ and $z_2$:

<figure><img src="image_34.png" /></figure>

### Meaning of $z_1$:

$$z_1 = a_2 \dot{x}_2 + a_3 \dot{x}_3$$

### Meaning of $z_2$:

$$\begin{aligned} z_2 &= b_1 \dot{z}_1 + b_2 \dot{x}_2 + b_3 \dot{x}_3 \\ &= b_2 \dot{x}_2 + b_3 \dot{x}_3 + b_1 a_2 \ddot{x}_2 + b_1 a_3 \ddot{x}_3 \end{aligned}$$

⇒ latent variables take into account **time delays** of $x_2$ and $x_3$
<figure><img src="image_35.png" /></figure>

Different models are tested:

- from $\mathbf{x} = [\mathbf{x}_2, \mathbf{x}_3, \mathbf{z}_1]$

- to $\mathbf{x} = [\mathbf{x}_2, \mathbf{x}_3, \mathbf{z}_1, \dots, \mathbf{z}_{20}]$
<figure><img src="image_36.png" /></figure>

Optimal model in terms of likelihood:

- with $\mathbf{x} = [\mathbf{x}_2, \mathbf{x}_3, \mathbf{z}_1, \dots, \mathbf{z}_6]$

- six latent variables are necessary
<figure><img src="image_37.png" /></figure>

Recursive predictions, starting from an initial time t:

$$E[x_{t+1}|y_1, \dots, y_t] = \text{ME}[x_t|y_1, \dots, y_t]$$

$$\text{Cov}[x_{t+1}|y_1, \dots, y_t] = \text{MCov}[x_t|y_1, \dots, y_t]M^\top + Q$$
Example of probabilistic forecasts (data- VS model-driven)

<figure><img src="image_38.png" /></figure>
Performance metrics (data- VS model-driven)

<figure><img src="image_39.png" /></figure>
Performance metrics (data- VS model-driven)

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>
## Context:

* pure data-driven predictions (i.e., without physical equations)

* missing components of the system (e.g., $x_1$ of Lorenz-63)

* simple assumptions (linearity, Gaussianity, constant operators)

<figure><img src="image_42.png" /></figure>
Key messages:

* fruitful combination of DA + ML

* inference of latent variables (interesting physical meaning)

* precise and reliable probabilistic forecasts

<figure><img src="image_43.png" /></figure>

Source: Brajard et al. 2020, Journal of Computational Science
Perspectives:

* application to climate indices (e.g., ENSO, AMO, NAO, etc.)

* application to high-dimensional data

* learn time-dependent $M(t)$ and $Q(t)$ operators

<figure><img src="image_44.png" /></figure>
Interested in this topic? Collaborations are welcome!

Contact:

• pierre.tandeo@imt-atlantique.fr

• www.tandeo.wordpress.com

<figure><img src="image_45.png" /><figcaption>Figure: Nice places to visit near Brest, Bretagne, France</figcaption></figure>
# Meaning of $z_1$

## Meaning of $z_1$:

* revealed by **symbolic regression**

* between $z_1$ and usual transformations of $x_2$ and $x_3$

<table><tr><td>$z_1$</td><td>$x_2$</td><td>$x_3$</td><td>$x_2x_3$</td><td>$x_2^2$</td><td>$x_3^2$</td><td>$\sqrt{x_2}$</td><td>$\sqrt{x_3}$</td><td>...</td><td>$\dot{x}_2$</td><td>$\dot{x}_3$</td></tr></table>

<figure><img src="image_46.png" /></figure>
Link between latent variables and $x_1$

Expression of $z_1$:

$$z_1 = a_2 \dot{x}_2 + a_3 \dot{x}_3$$

Lorenz-63 system:

$$\dot{x}_2 = x_1(\rho - x_3) - x_2$$

$$\dot{x}_3 = x_1x_2 - \beta x_3$$

Approximation of $x_1$:

$$\hat{x}_1 = \frac{z_1 + a_2x_2 + a_3\beta x_3}{a_2\rho - a_2x_3 + a_3x_2}$$

<figure><img src="image_47.png" /></figure>

⇒ $x_1$ is partially retrieved (division by zero in specific regions)
Meaning of $z_1$:

$$z_1 = a_2 \dot{x}_2 + a_3 \dot{x}_3$$

<figure><img src="image_48.png" /></figure>

⇒ **two clusters** of solutions ($a_2 \approx 0$ and $a_3 \approx 0$)

⇒ $x_3$ **more important** than $x_2$ in the Lorenz-63 system