Discover governing differential equations from evolving systems

We start by collecting the spatiotemporal series data at $m$ time points and $n$ spatial locations of the state variables. For each state variable, the captured state measurements are represented as a single column state vector $u\in C^{nm\times 1}$. Based on all the observables, e.g., variables $x$, $y$, and $z$ in Lorenz systems, a series of functional terms associated with these quantities can be calculated and then reshaped into a single column as well, such as $x^{2}$, $xy$, $xyz$, $x^{2}y$, etc. Likewise, partial differential terms should be considered in the candidate library if PDEs govern the dynamics. Furthermore, “1” denotes the constant term that possibly appears in equations. The compositive function library $\Theta\left(U\right)\in C^{nm\times p}$ is a matrix that contains $p$ designed functional terms. Note that the computed time derivative $u_{t}\in C^{nm\times 1}$ is also a single-column vector presented on the left-hand side of Eq. (2). The constant term in the library $\Theta$ sufficiently considers the bias term in the governing differential equation so that the model can be regarded as an unbiased representation of dynamics. For O-SINDy, an important prerequisite for revealing the correct governing differential equation is that the candidate function library contains all the members which constitute the concise dynamics expression, so as to ensure that the exact sparse coefficient $\Xi$ is computed under iterations.

Formulating a hypothesis that the governing differential equations are evolving at any moment, we model each example as a dynamic process to simulate the arrival of streaming data. Given a set of time-series state measurements at a fixed number of spatial locations in $x$, the goal is to construct the form of $N(\bullet)$ online. To fix this issue, the core of our innovation is utilizing real-time streaming data to denote the loss function. Considering a coefficient vector in $\Xi$, $\xi_{j}(j=1,2,\ldots,d)$, which corresponds to the specific state variable $u$, the fitness function is designed as follows:

where $\lambda_{1}\geq 0$ denotes the regularization coefficient, $n$ is the number of positions, $m$ denotes the length of the time series, $\xi_{j}=\left[\xi_{j1},\xi_{j2},\ldots,\xi_{jp}\right]^{T}$, and $L_{i}\left(\xi_{j}\right)$ is defined as follows:

where ${\dot{U}}_{i}$ is the time derivative of the $i$th state observation, $\Theta_{i}$ is the $i$th row of the common library $\Theta$ that represents all candidate function values for the ith data. The sparsity constraints mean that the coefficient vector $\xi_{j}$ is sparse with only a few non-zero entries, each representing a subsistent item in the function library. Subsequently, we exploit the follow-the-regularized-leader (FTRL)-Proximal [37, 38, 39] style methodology to optimize the outcome by considering solely one instance each time.

Leveraging the fact that each instance is considered individually, we rewrite the loss function and define a single loss term (see Eq. (5)). For example, if the first state measurement $u_{1}$ is available, the loss is defined as: $L_{1}\left(\xi_{j}\right)+\lambda_{1}\|\xi_{j}\|_{1}=0.5(\Theta_{1}\xi_{j}-{\dot{U}}_{1})^{2}+\lambda_{1}\|\xi_{j}\|_{1}$. Next, if the second instance of $u_{2}$ arrives, the loss is defined as: $L_{1}\left(\xi_{j}\right)+L_{2}\left(\xi_{j}\right)+\lambda_{1}\|\xi_{j}\|_{1}=0.5(\Theta_{1}\xi_{j}-{\dot{U}}_{1})^{2}+0.5(\Theta_{2}\xi_{j}-{\dot{U}}_{2})^{2}+\lambda_{1}\|\xi_{j}\|_{1}$. In this way, the computer only needs to store information about a single example, thus relieving memory stress. On the $i$th sample, the gradient is calculated as follows:

Normally, the online gradient descent algorithm can be used to update the coefficient vector $\xi_{j}^{i}$ after the arrival of the $i$th data by using:

However, this method has been proven to lack general applicability [38]. Correspondingly, the FTRL-proximal style approach is introduced to solve the online problem. We use the following equation to update the coefficient vector $\xi_{j}^{i}$:

where $\lambda_{1}$ and $\lambda_{2}$ both denote the regularization coefficient, which is the positive constant, while $\sigma_{k}$ is on the $i$th learning rate $\eta_{i}$ aspect, defined by:

Given the gradient ${\nabla L}_{i}\left(\xi_{j}\right)$ a shorthand $g_{i}$, we set $g_{1:i}=\sum_{k=1}^{i}g_{k}$. Then the update in Eq. (8) can be rewritten as follows:

where the last term $0.5\sum_{k=1}^{i}\sigma_{k}\|\xi_{j}^{k}\|_{2}^{2}$ is a constant with respect to $\xi_{j}$ and negligibly impacts the update process. Let $Z_{i}=g_{1:i}-\sum_{k=1}^{i}{\sigma_{k}\xi_{j}^{k}}=Z_{i-1}+g_{i}-\sigma_{i}\xi_{j}^{i}$, and we ignore the last term and rewrite Eq. (II.2) as:

Let $\xi_{js}$ and $Z_{i,s}(s=1,2,\ldots,p)$ represent the $s$th element of the vector $\xi_{j}$ and $Z_{i}$, respectively. Equation (9) is generalized as

Aiming to simplify the expression, we suppose $w^{\ast}$ to be the optimal solution of $\xi_{js}^{i}$ and $\Phi\in\nabla w^{\ast}$ to be its gradient. Then, Eq. (13) is satisfied. Equation (17) shows the subdifferential of $w=sgn(w)$.

The solution to $w^{\ast}$ can be discussed in three cases [40]:

1. 1.

   If $|Z_{i,s}|\leq\lambda_{1}$,

   1. (a)

      If $w^{\ast}=0$, then $sgn(0)\in(-1,1)$ and Eq. (13) is satisfied.

   2. (b)

      If $w^{\ast}>0$, then $Z_{i,s}+\lambda_{1}sgn(w^{\ast})=Z_{i,s}+\lambda_{1}\geq 0$ and $(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k})w^{\ast}>0$. Equation (13) does not hold.

   3. (c)

      If $w^{\ast}<0$, then $Z_{i,s}+\lambda_{1}sgn(w^{\ast})=Z_{i,s}-\lambda_{1}\leq 0$ and $(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k})w^{\ast}<0$. Equation (13) does not hold.

2. 2.

   If $|Z_{i,s}|<-\lambda_{1}$,

   1. (a)

      If $w^{\ast}=0$, then $sgn(0)\in(-1,1)$ and Eq. (13) does not hold.

   2. (b)

      If $w^{\ast}>0$, then $Z_{i,s}+\lambda_{1}sgn(w^{\ast})=Z_{i,s}+\lambda_{1}<0$ and $(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k})w^{\ast}>0$. Equation (13) holds and the solution is:

      $$w^{\ast}=-\left(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k}\right)^{-1}\left(Z_{i,s}+\lambda_{1}\right).$$

      (18)

   3. (c)

      If $w^{\ast}<0$, then $Z_{i,s}+\lambda_{1}sgn(w^{\ast})=Z_{i,s}-\lambda_{1}<0$ and $(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k})w^{\ast}<0$. Equation (13) does not hold.

3. 3.

   If $|Z_{i,s}|>\lambda_{1}$,

   1. (a)

      If $w^{\ast}=0$, then $sgn(0)\in(-1,1)$ and Eq. (13) does not hold.

   2. (b)

      If $w^{\ast}>0$, then $Z_{i,s}+\lambda_{1}sgn(w^{\ast})=Z_{i,s}+\lambda_{1}>0$ and $(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k})w^{\ast}>0$. Equation (13) does not hold.

   3. (c)

      If $w^{\ast}<0$, then $Z_{i,s}+\lambda_{1}sgn(w^{\ast})=Z_{i,s}-\lambda_{1}>0$ and $(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k})w^{\ast}<0$. Equation (13) holds and the solution is:

      $$w^{\ast}=-\left(\lambda_{2}+\sum_{k=1}^{i}\sigma_{k}\right)^{-1}\left(Z_{i,s}-\lambda_{1}\right).$$

      (19)

The above discussion in different estimation scenarios is summarized in Eq. (22), where $w^{\ast}=\xi_{js}^{i}$.

Moreover, we set a remarkable learning rate

where $\alpha$ and $\beta$ are both positive constants, $g_{k,s}$ denotes the $s$th element of the gradient $g_{k}$ and $\eta_{is}$ is the learning rate of $g_{k,s}$. Thus, Eq. (9) can be rewritten as follows:

The final form of the closed solution is obtained and expressed as

where $\alpha$ and $\beta$ are both positive constants, $g_{k,s}$ denotes the $s$th element of the gradient $g_{k}$.

Notably, $\xi_{js}$ is always updated according to the current input instance, thereby guiding the detection of change points if there exists any variation in the evolving system. Additionally, the mechanism frees the computer from storing the whole large-scale data. It should be highlighted that the method suggested above is entirely different from the batch learning methods, which update gradients depending on all available measurements. O-SINDy ¹¹1Available code: https://github.com/xiaoyuans/O-SINDy. takes advantage of a single available instance and updates gradients based on the current loss $L_{i}\left(\xi_{j}\right)(i=1,2,\ldots,mn)$ in each iteration, thereby being easily extended to handle evolving systems. Note that the optimum selection of hyper-parameters in each case is determined by grid search ($10^{i},i=-8,-7,…,1$).

We first consider an evolving two-dimensional damped harmonic oscillator that changes from cubic to linear.

A particle in simple harmonic vibration is known as a harmonic oscillator, whose motion is the simplest ideal vibration model. Fig. 2 illustrates the dynamic data and the detection of change points. Firstly, we create the solution to Eq. (32), $U_{1}$, with 2,500 timesteps and the initial data $[2,0]$. Then, we utilize the last instance in $U_{1}$ as the initial data and create the solution to Eq. (31) with 2,500 timesteps as well, recorded as $U_{2}$. Specifically, the change point arises when the system changes, and then the original distribution of generated data varies. Therefore, the corresponding loss value gradually decreases and stabilizes unless the change point appears. In this context, the coefficients in the former governing equation and the change location are simultaneously identified according to the sudden increase of loss value.

With an augmented nonlinear library including polynomials up to the fifth order, the change point is spotted in evolving circumstances. The loss curve tends to stabilize faster after resetting the training coefficient, thereby accelerating the convergence of O-SINDy because the former result may mislead the training process.

The transition from the KdV to the Kuramoto-Sivashinsky equation is simulated as an illustrative example that exhibits qualitatively different dynamics as the system evolves. Summary of O-SINDy and the state-of-the-art methods in SINDy for this identification in Table 2 demonstrates that our approach works well. Most data-driven methods regard batch measurements as a whole and update gradients integrally to obtain an ephemeral fitting. Owing to the ignorance of variations, existing methodologies are biased to strike a compromise solution, so that strenuous to identify the changes from consecutively generated data. Quite the contrary, O-SINDy succeeds in simultaneously estimating forms and identifying changes of the governing equations.

Moreover, we investigate the Susceptible-Infected-Recovered (SIR) disease model with time-varying transmission rates, widely studied in epidemiology [45, 46]. The following equation is a mathematical description of the SIR model.

It is a time-dependent hybrid dynamical system, where $v=d=0.0027$, $N=1000$ is the total population, $\gamma=0.2$, $b=0.8$, and $\hat{\beta}=9.336$ is a base transmission rate. We simulate the model for three years within the initial condition at $[S_{0},I_{0},R_{0}]=[12,13,975]$, recording at a daily interval. A random perturbation is added to the start of each season with the same probability.

To illustrate the dynamics of the SIR model, we display the recovered population by subtracting 920 and scaling the loss in O-SINDy by the maximum value. As shown in Fig. 3, O-SINDy captures the transition between states by the dramatically increasing loss. In this example, the perturbation is able to be identified utilizing solely the information of the recovered population, which is ignored in Hybrid-Sparse Identification of Nonlinear Dynamics (Hybrid-SINDy)[46]. It is an offline batch learning method and outputs the governing equation of a single day by repeatedly clustering the whole time-series data during the model selection step. The identified governing equations of four random days in the four seasons are summarized in Table 3. The results show that O-SINDy provides more precise reconstruction accuracy for a single day. Specifically, instances of the same coefficients share the hyper-parameters in O-SINDy. At the same time, Hybrid-SINDy selects a model of the highest frequency from results recovered by SINDy with different hyper-parameters for each certain instance. Moreover, the effect of cluster size is particularly vital to solutions in Hybrid-SINDy.

We also test O-SINDy on two switching linear systems, including synthetic and real-world datasets. Define $x_{t}\in R^{N}$ as the sampled trajectory of a dynamical system, $N$ is the feature dimension. Time-varying autoregressive model with low-rank tensors (TVART) [47] holds the assumption, $x_{t+1}=A_{t}x_{t}$, where $A_{t}$ is constant within a time window but varies in the switching system. On the contrary, O-SINDy discards the time window and admits that $A_{t}$ may change at every time point.