Reduced Order Data-Driven Twin Models for Nonlinear PDEs by Randomized Koopman Orthogonal Decomposition and Explainable Deep Learning

1. Introduction

In recent years, digital twin models have emerged as a transformative tool for understanding, simulating, and predicting the behavior of complex physical systems across various domains [1,2,3]. In a broader sense, a digital twin is a virtual representation of a physical system that is continuously updated with real-time data and is capable of capturing the system’s evolution over time [4].

While digital twins have demonstrated considerable utility in the context of engineered systems with well-defined models and sensing infrastructures, their application in replicating nonlinear dynamical systems remains limited by the inherent complexity and nonlinearity of the underlying dynamics. Furthermore, traditional digital twins often require full physical modeling and high computational resources, which limits their scalability and adaptability to purely data-driven or partially observed systems.

To address these challenges, this paper introduces a novel framework, which is termed the data-driven twin model (DTM). This new paradigm is inspired by the digital twin concept but diverges fundamentally in its construction and purpose: rather than embedding a full physics-based model, the data-driven twin model is designed to learn and replicate the nonlinear behavior of a complex dynamical system directly from data with very high fidelity and reduced complexity. Unlike standard machine learning models, which often operate as black boxes, the data-twin model aims to preserve interpretability by decoupling the system’s evolution into a stationary component and a dynamic component in the mathematical framework of modal decomposition [5].

The concept of a data-driven twin model is, to our knowledge, entirely novel and particularly significant in the context of unsupervised learning for dynamical systems. By separating stationary features from time-varying dynamics, the data-twin model facilitates more efficient learning, better generalization to various nonlinear dynamics, and offers insights into the intrinsic structure of the system dynamics. This decomposition is not only beneficial for reducing the learning burden but also opens up new possibilities for transfer learning, system identification, and control.

In this research, the data-driven twin model concept is introduced as a reduced order model that has the main feature to mirror the original process behavior with significantly reduced computational complexity. This advancement addresses a fundamental gap in reduced order modeling: the reconciliation of real-time adaptability with time simulation of nonlinear feature. Associating a nonlinear dynamical process with a data-twin of reduced complexity, assisted by explainable deep learning mathematical tools, has the significant advantage to map the dynamics with highest accuracy and reduced costs in CPU time and hardware to timescales over which that suffers significantly changes and so it could be difficult to explore.

1.1. Modal Decomposition Literature Review and Previous Work

Modal decomposition [5] represents a generic mathematical framework for the characterization of nonlinear dynamical systems in the context of reduced-order modeling [6,7,8,9,10,11,12].

Among the prominent modal decomposition algorithms to produce models of reduced computational complexity, Proper Orthogonal Decomposition and Dynamic Mode Decomposition have become widely adopted in modeling complex dynamical systems.

Proper Orthogonal Decomposition (POD) has its roots in the foundational works of Loéve (1945) and Karhunen (1946) [13,14], who independently developed what is now known as the Karhunen-Loéve (K–L) expansion, a statistical method for representing stochastic processes through orthogonal basis functions. This framework was later extended and formalized in the context of fluid mechanics by Lumley (1970) [15], who introduced POD as a numerical tool for extracting dominant coherent structures from turbulent flow data.

POD has shown effectiveness across a wide range of fields, from turbulent and convection-dominated flows [16,17,18,19] to engineering and oceanography [20,21,22], particularly when combined with data assimilation [23], control [24], or machine learning approaches especially artificial neural networks [25,26,27,28,29]. POD’s reliance on energy-based criteria for selecting basis functions and using truncation-based approximation introduces significant limitations, particularly in the context of highly complex or nonlinear systems. These methods often fail to capture low-energy yet dynamically important modes, leading to reduced accuracy and insufficient representation of critical system behaviors [30,31]. POD suffers also from limited precision in forecasting because it is inherently a projection-based method. By projecting complex, high-dimensional dynamics onto a reduced set of orthogonal modes derived from snapshot data, POD captures only the dominant patterns of variability. As a result, finer scale or transient features, especially those not well represented in the training data, can be lost or poorly reconstructed, leading to reduced accuracy in long-term or highly nonlinear forecasts.

Originally introduced by Koopman (1931), Koopman Operator Theory [32,33] provides a linear, operator-theoretic framework for analyzing and deriving reduced-order models of complex nonlinear dynamical systems. Its modern resurgence began with the pioneering work of Mezić and Banaszuk (2004) in fluid dynamics applications [34,35], laying the foundation for practical, data-driven modal decomposition analysis. Subsequently, Schmid and Sesterhenn (2010) introduced the Dynamic Mode Decomposition (DMD) technique [36,37], a numerical method inspired by Koopman theory that has since been widely adopted across various scientific and engineering disciplines [6,10,38,39]. The research group led by J. Nathan Kutz further advanced the field by developing robust Koopman-based methods for data-driven modeling and system identification [8,40,41,42,43].

Today, Koopman operator theory associated with DMD stands as a central tool in modern dynamical systems analysis and is increasingly recognized as a powerful and influential mathematical framework of the 21^st century. The development of numerous DMD variants reflects the ongoing efforts to enhance its robustness and applicability: optimized DMD [44], exact DMD [40], sparsity promoting DMD [45], multi-resolution DMD [42], extended DMD [46], recursive dynamic mode decomposition [47], DMD with control [48], randomized low-rank DMD [49], dynamic mode decomposition with core sketch [50], bilinear dynamic mode decomposition [51], higher order dynamic mode decomposition [52].

DMD captures both spatial and temporal dynamics through modes associated with specific frequencies and growth rates. Yet, standard DMD techniques can suffer from reduced precision if the modes are not selected properly. The accuracy and reliability of the reconstructed or forecasted dynamics heavily depend on identifying the most relevant modes, which is a nontrivial task. Establishing additional mode selection criteria often requires manual tuning, prior knowledge, or heuristic methods, making the process complex and computationally expensive. This not only increases CPU time but also undermines the self-contained nature of the algorithm, limiting its applicability in fully unsupervised learning scenarios.

Author’s contributions have been particularly significant in advancing algorithms based on Koopman operator theory, especially for challenging scenarios involving high-dimensional data and complex temporal behavior. The author has made significant contributions to the development of mode selection strategies in DMD, introducing selection criteria based on amplitude-weighted growth rates and energy-Strouhal number weighting, aimed at systematically extracting the most dynamically relevant modes [53,54].

Moreover, to address the inherent limitations of standard Dynamic Mode Decomposition, specifically, the computational burden of mode selection, the adaptive randomized dynamic mode decomposition (ARDMD) was introduced in [55], offering a more efficient and scalable alternative for extracting dynamic features from high-dimensional data. This innovative algorithm improves efficiency and scalability in reduced-order modeling, particularly in fluid dynamics applications [56].

The techniques introduced in these previous works addressed the persistent challenge of identifying dominant modes in nonlinear and unsteady systems, crucial for accurate reconstruction and forecasting. Author’s work has also extended DMD’s utility into emerging areas like epidemiological modeling [57,58], showcasing its broader societal relevance.

1.2. Novelty and Key Advantages of the Present Research

Recognizing the limitations of standard DMD, particularly the lack of orthogonality and computational cost associated with mode selection, this work proposes a novel framework for constructing reduced-order models of nonlinear dynamical systems. The central objective is to replicate the full system dynamics with significantly lower model complexity, enabling efficient, interpretable, and scalable modeling, i.e. a data-driven twin modeling approach based on modern Koopman operator theory.

In recent years, explainable deep learning has emerged as a preferable alternative to traditional AI approaches that rely on latent variables with unclear or opaque mechanisms. Unlike classical models, whose internal representations often lack interpretability, explainable deep learning provides a framework in which the model’s behavior can be understood, analyzed, and trusted. This work presents an application of explainable deep learning in a controlled, well-understood, and mathematically rigorous setting. The approach not only ensures transparency in the modeling process but also lays a solid foundation for future extensions to unsupervised learning tasks, where the need for reliable and interpretable representations is even more critical.

The core contributions and key advantages include:

1. First mathematical framework for a compact data-driven twin model combining high dynamic fidelity with low computational cost.

2. Fully self-consistent algorithm requiring no heuristic tuning for the reduced-order model parameters.

3. Integration of explainable deep learning for transparent, adaptive modeling in a rigorous, mathematical and controlled setting.

This framework not only improves the interpretability and precision of modal decompositions but also enhances the practical applicability of Koopman operator theory in data-driven modeling.

The structure of the article is organized as follows: Section 2 outlines the mathematical foundations of reduced-order modeling grounded in Koopman operator theory. Section 3 addresses the computational aspects and introduces a novel modal decomposition algorithm, accompanied by a rigorous mathematical justification. Section 4 demonstrates the application of the proposed technique to shock wave phenomena through a series of three experiments of increasing complexity. Section 5 presents the corresponding numerical results. Finally, Section 6 concludes the study with a summary of findings and potential future directions.

2. Mathematical Background

2.1. Computational Complexity Reduction of PDEs by Reduced-Order Modeling

Partial differential equations (PDEs) serve as a fundamental framework for accurately modeling nonlinear dynamical phenomena, owing to their capacity to represent spatial-temporal variations in complex systems. Nevertheless, the numerical solution of PDEs is often computationally demanding, particularly in high-dimensional and strongly nonlinear settings. This computational burden necessitates strategies for complexity reduction to enable practical simulation and analysis. In this context, by reducing the dimensionality of the system while retaining its essential dynamics, reduced-order modeling significantly lowers computational cost and facilitates efficient simulation and analysis.

Suppose that $\Omega \subset {\mathbb{R}}^d$ represents the bounded computational domain and let the Hilbert space $L^2\left(\Omega \right)$ of square integrable functions on $\Omega$:

$$L^2\left(\Omega \right)=\left\{\varphi :\Omega \to \mathbb{R}\left|{\int }_{\Omega }{\left|\varphi \right|}^{2}d\Omega <\infty \right.\right\}$$

(1)

to be endowed with the inner product:

$$〈{\varphi }_i,{\varphi }_j〉={\int }_{\Omega }{\varphi }_i\overline{{\varphi }_j}d\Omega for{\varphi }_i,{\varphi }_j\in L^2\left(\Omega \right)$$

(2)

and the induced norm $∥\varphi ∥=\sqrt{〈\varphi ,\varphi 〉}$ for $\varphi \in L^2\left(\Omega \right)$.

Consider a nonlinear PDE for state variable $u\left(x,t\right)$, $x\in {\Omega }_x\subset {\mathbb{R}}^d$, $t\in {\Omega }_t=\left[t_0,T\right]$, on domain $D={\Omega }_x\times {\Omega }_t$, of the form:

$${\displaystyle \frac{\partial u(x,t)}{\partial t}}=\mathcal{N}\left[u(x,t)\right],$$

(3)

where $\mathcal{N}$ is a nonlinear differential operator, with the initial condition and Robin boundary conditions given in the following form:

$$IC:u\left(x,t_0\right)=u_0\left(x\right),$$

(4)

$$BC:\left\{\begin{array}{c}{\alpha }_{1}u\left(a,t\right)+{\alpha }_2{u}_x\left(a,t\right)=g_a\left(t\right),\hfill \\ {\beta }_{1}u\left(b,t\right)+{\beta }_2{u}_x\left(b,t\right)=g_b\left(t\right).\hfill \end{array}\right.$$

(5)

Over the past half-century, partial differential equations have been predominantly solved using numerical techniques such as the finite difference and finite element methods. While these computational tools have significantly advanced our ability to approximate complex PDE solutions, they also introduce substantial implementation complexity and often require a solid understanding of the underlying analytical structure of the equations. In modern computational mathematics, reduced-order modeling has emerged as a powerful approach to mitigate computational demands by simplifying high-dimensional PDE systems, enabling more efficient and tractable solutions without sacrificing essential accuracy.

Definition 1. 

The principle of reduced-order modeling aims to finding an approximation solution of the form

$$\left\{\begin{array}{c}u\left(x,t\right)=\underset{p\to \infty }{lim}{\displaystyle \sum _{j=1}^p}{a}_j\left(t\right){\psi }_j\left(x\right),t\in {\mathbb{R}}_{\ge 0}\hfill \\ {\displaystyle \frac{d{a}_j\left(t\right)}{dt}}=\mathcal{M}\left(a_j\left(t\right),t\right),a_j\left(t_0\right)=a_j^0\hfill \end{array}\right.$$

(6)

which assumes the existence of a basis functions $\Psi =\left\{{\psi }_j\left(x\right)\right\},j=1,..p$ and a locally Lipschitz continuous map $\mathcal{M}$, assuring preservation of small global approximation error compared to the true solution of (3) and dynamic stability and computational stability of the reduced-order solution. Expression (6) is called the reduced-order model of the PDE (3).

The basis functions $\Psi$ are called modes or shape-functions, therefore the representation (6) is called modal decomposition. Two primary goals are achieved through this modal decomposition approach: (1) the decoupling of spatial structures from the unsteady (time-dependent) components of the solution, and (2) the modeling of the original system’s nonlinearity through a lower-dimensional temporal dynamical system. This separation not only reduces computational complexity but also preserves the essential dynamics of the original PDE system in a more efficient and interpretable form.

2.2. Koopman Operator Framework for Modal Decomposition

Koopman operator theory [32,33,59] offers a rigorous mathematical framework for the modal decomposition of nonlinear dynamical systems by enabling their representation in terms of linear, infinite-dimensional operators acting on observables. In this paper, the Koopman framework is adapted to the context of nonlinear partial differential equations, extending its application to infinite-dimensional dynamical systems governed by PDEs.

Let $\mathcal{U}\left(t\right)$ denote the evolution operator of (3) such that:

$$u(x,t)=\mathcal{U}\left(t\right)u(x,t_0),$$

(7)

i.e., $\mathcal{U}\left(t\right)$ maps the initial condition $u(x,t_0)$ to the state at time t.

Definition 2. 

An observable of the PDE (3) is defined as a possibly nonlinear functional on the state space (e.g., a function space such as $\mathcal{H}=L^2\left(\Omega \right)$)

$$g:\mathcal{H}\to {\mathbb{C}}^m,$$

(8)

that maps the infinite-dimensional state u to a measurable quantity.

Remark 1. 

The choice of observable g in the Koopman framework is flexible and can range from simple to highly structured representations. A common example is theidentity observable, where $g\left(u\right)=u$, mapping the system state directly to itself. Alternatively, one may consider an infinite-dimensional observable vector, comprising a set of functions of the state, such as

$$g\left(u\right)={[u,\kappa u,u^2,u^3,sin\left(u\right),\dots ]}^T,$$

(9)

where κ denotes a linear operator (e.g., a differential operator), and the higher-order terms capture nonlinear features of the system. This extended observable space facilitates a richer linear representation of the underlying nonlinear dynamics through the Koopman operator.

Definition 3. 

The semigroup of Koopman operators ${\mathcal{K}}^t$ acts on observable functions by composition with the evolution operator of the states:

$$\left({\mathcal{K}}^{t}g\right)\left(u_0\right)=g\left(\mathcal{U}\left(t\right)u_0\right),$$

(10)

where $u_0=u(x,t_0)$ is the initial condition.

As a result, the Koopman operator is also known as a composition operator.

In practice, the Koopman operator is considered to be the dual, or more precisely, the left adjoint, of the Ruelle-Perron-Frobenius transfer operator [60,61,62], which acts on the space of probability density functions.

Remark 2. 

The Koopman operator is linear, even though the underlying dynamics may be nonlinear:

$${\mathcal{K}}^t(a{g}_1+b{g}_2)=a{\mathcal{K}}^t{g}_1+b{\mathcal{K}}^t{g}_2.$$

(11)

Remark 3. 

The infinitesimal generator $\mathcal{L}$ of the Koopman semigroup satisfies:

$${\displaystyle \frac{d}{dt}}g\left(u\right)=\mathcal{L}g\left(u\right),$$

(12)

where

$$\mathcal{L}g\left(u\right)=〈{\displaystyle \frac{\delta g}{\delta u}},\mathcal{N}\left[u\right]〉$$

(13)

and ${\displaystyle \frac{\delta g}{\delta u}}$ denotes the functional (Fréchet) derivative of g with respect to u, and the inner product is taken over the spatial domain.

Since its introduction in 1931, the operator proposed by Koopman [33,59] remained a purely theoretical construct. In 1985, Lasota and Mackey brought renewed scholarly attention to the Koopman operator [63]. Notably, they were the first to explicitly introduce and formalize the term "Koopman operator" in the literature, thereby establishing a standardized nomenclature for this concept. This terminology was further reinforced in the second edition of their book [64] published in 1994, where the use of the Koopman operator as the adjoint to the Frobenius-Perron operator was maintained and elaborated upon functional analysis in space $L^{\infty }$.

However, at that time, computational resources were still insufficient for the numerical treatment of large-scale complex systems. The practical potential of Koopman Operator Theory began to be realized in the early 2000s. Notably, Mezić and Banaszuk applied it in fluid dynamics [34]. The paper of Mezić [35] is the key modern reference that introduced a rigorous spectral framework for the analysis of nonlinear dynamical systems based on the Koopman operator theory. He formalized the Koopman spectral decomposition, and was the first to introduce the decomposition of observables into Koopman eigenfunctions and associated modes, laying the foundation for a spectral analysis of nonlinear dynamical systems through the Koopman operator framework. The spectral properties of the Koopman operator have been the subject of ongoing investigation, as exemplified in [41,44,65] .

Following Koopman operator theory, the solution $u(x,t)$ can be expressed as a modal decomposition of the form:

$$u(x,t)=\sum _{j=1}^{\infty }{a}_j\left(t\right){\Phi }_j\left(x\right),$$

(14)

where ${\Phi }_j\left(x\right)$ are spatial Koopman modes, and $a_j\left(t\right)$ are time-dependent coefficients.

The evolution of the modal coefficients $\mathbf{a}\left(t\right)={[a_1\left(t\right),a_2\left(t\right),\dots ]}^T$ is governed by a possibly nonlinear system:

$${\displaystyle \frac{d\mathbf{a}\left(t\right)}{dt}}=\mathcal{M}(\mathbf{a}\left(t\right),t),$$

(15)

where $\mathcal{M}$ is a nonlinear operator that captures the dynamics of the time-dependent coefficients. Thus, the reduced-order solution of form (6) to the PDE is achievable form the perspective of Koopman mode decomposition.

Koopman modes are generally not unique, and their characterization depends on several interrelated factors. First, the definition of Koopman modes is inherently linked to the choice of observables, as the Koopman operator acts on a space of functions rather than directly on the state space. Different selections of observables can lead to different spectral representations, thereby affecting the resulting modes.

This intrinsic non-uniqueness is further compounded in practical settings where Koopman modes are often approximated using numerical techniques such as Dynamic Mode Decomposition (DMD), introduced by P. J. Schmid and J. Sesterhenn in 2008 [36,37]. Such methods introduce additional variability due to factors including data noise, finite sample sizes, and specific algorithmic implementations. Moreover, in cases where the Koopman operator exhibits a degenerate spectrum [10,66], the associated eigenspaces are typically multidimensional, allowing for multiple valid decompositions and further contributing to the non-uniqueness of the Koopman modes.

This paper further introduces a novel algorithm derived from the principles of Dynamic Mode Decomposition, yet fundamentally distinct in its formulation, specifically designed for the construction of a data-driven twin model. This new method computes Koopman modes and their associated temporal coefficients by leveraging the Koopman propagator operator projected onto a Krylov subspace of optimally selected rank. This approach enables a more structured and efficient extraction of spectral components by constraining the dynamics within a finite-dimensional subspace tailored to capture the most relevant features of the system’s evolution.

3. Computational Aspects

This section introduces a novel methodology for constructing a data-driven twin model of a partial differential equation using Koopman mode decomposition. The proposed approach comprises two primary stages. The first stage, referred to as the offline phase, involves the extraction of Koopman modes and the corresponding modal decomposition coefficients that characterize the state evolution of the system. The second stage, known as the online phase, focuses on identifying a nonlinear continuous-time dynamical model that governs the temporal behavior of the data-twin, thereby enabling real-time simulation and prediction. A detailed account of the computational implementation of the proposed methodology is presented hereinafter.

3.1. Offline Phase: Koopman Randomized Orthogonal Decomposition Algorithm

The data $u_i\left(x,t\right)=u\left(x,t_i,\right),t_i=i\Delta t,i=0,...,N_t$, represent measurements of the PDE solution at the constant sampling time $\Delta t$, x representing the Cartesian spatial coordinate.

The data matrix whose columns represent the individual data samples is called the snapshot matrix

$$V=\left[\begin{array}{cccc}{u}_0& u_1& ...& u_{N_t}\end{array}\right]\in {\mathbb{R}}^{N_x\times (N_t+1)}.$$

(16)

Each column $u_i\in {\mathbb{R}}^{N_x}$ is a vector with $N_x$ components, representing the spatial measurements corresponding to the $N_t+1$ time instances.

We aim to find the data-driven twin representation of the solution at every time step $\left\{{\mathrm{t}}_1,...,{\mathrm{t}}_{N_t}\right\}$, according to the following relation:

$$u_i=u\left(x,t_i\right)=\sum _{j=1}^{N_{DTM}}\underset{Modalamplitudes}{\underbrace{a_j\left(t_i\right)}}\stackrel{LeadingKoopmanmodes}{\overbrace{{\varphi }_j\left(x\right)}},t_i\in \left\{{\mathrm{t}}_1,...,{\mathrm{t}}_{N_t}\right\},i=1,...,N_t,$$

(17)

where $\Phi ={\left\{{\varphi }_j\right\}}_{j=1}^{N_{DTM}}$ represent the extracted Koopman modes base functions, which we call the leading Koopman modes, $N_{DTM}\ll min\left(N_x,N_t+1\right)$ represents the number of terms in the representation (17) which we impose to be minimal and $a_j\left(t_i\right)$ represent the modal growing amplitudes.

In classical Dynamic Mode Decomposition the Koopman modes are not orthogonal, which often necessitates a large number of modes to achieve an accurate modal decomposition. This lack of orthogonality represents a significant drawback, as it may lead to redundancy and reduced computational efficiency.

It is more convenient to seek an orthonormal base of Koopman modes

$$\Phi =\left\{{\varphi }_1,{\varphi }_2,...\right\},〈{\varphi }_i\left(x\right),{\varphi }_j\left(x\right)〉={\delta }_{ij},∥\varphi ∥=1,$$

(18)

where ${\delta }_{ij}$ is the Kronecker delta symbol, consisting of a minimum number of Koopman modes ${\varphi }_i\left(x\right)$, such that the approximation of $u\left(x,t\right)$ through this base is as good as possible, in order to create a data-twin model to the PDE, of reduced computational complexity.

The objective of the present technique is to represent with highest precision the original data through the data-twin model (17) having as few terms as possible. To achieve this goal, several major innovative features are introduced in the proposed algorithm:

• The method produces and extracts orthonormal Koopman modes, ensuring mutual orthogonality among the modes. As a result, a more compact representation is achieved, requiring a smaller number of modes to accurately capture the system dynamics.
• To avoid the computational burden associated with traditional high-dimensional algorithms, the proposed method employs a randomized singular value decomposition (RSVD) technique for dimensionality reduction. The incorporation of randomization offers a key advantage: it eliminates the need for additional selection criteria to identify shape modes, which is typically required in classical approaches such as DMD or POD. The resulting algorithm efficiently identifies the optimal reduced-order subspace that captures the dominant Koopman mode basis, thereby ensuring both computational efficiency and representational fidelity.
• The methodology aims to achieve maximal correlation and minimal reconstruction error between the data-driven twin model and the exact solution of the governing PDE.

The theoretical foundation of randomized singular value decomposition (RSVD) technique was laid in the seminal work of Halko, Martinsson, and Tropp (2011) [67] as part of a broader class of probabilistic algorithms for matrix approximation. Subsequent study of Erichson and Donovan (2016), demonstrated its practical efficiency in motion detection [49]. Bistrian and Navon (2017) extended its use for the first time in fluid dynamics [55,56]. Recently, the theory of randomization has been further employed in the development of reduced-order modeling techniques, particularly in conjunction with projection learning methods [68], block Krylov iteration schemes [69] or deep learning [70].

Proposition 1. (k-RSVD: Randomized Singular Value Decomposition of rank k)

Let $V_0\in {\mathbb{R}}^{N_x\times N_t}$ be a real-valued data matrix. Having imposed a target rank $k<N_t$, the Randomized Singular Value Decomposition of rank k (k-RSVD) produces k left/right singular vectors of ${\mathrm{V}}_0{{\mathrm{V}}_0}^H$ and ${{\mathrm{V}}_0}^H{\mathrm{V}}_0$, respectively, and has the following steps:

1. 

Generate a Gaussian random test matrix M of size $N_t\times k$.

2. 

Compute a compressed sampling matrix by multiplication of data matrix with random matrix $Q=V_{0}M$.

3. 

Project the data matrix to the smaller space $P=Q^H{V}_0$, where H denotes the conjugate transpose.

4. 

Produce the economy-size singular value decomposition of low-dimensional data matrix $\left[T,\Sigma ,W\right]=svd\left(P\right)$.

5. 

Compute the right singular vectors $U=QT$, $U\in {\mathbb{R}}^{N_x\times k}$, $\Sigma \in {\mathbb{R}}^{k\times k}$, $W\in {\mathbb{R}}^{N_t\times k}$.

Theorem 1. (Koopman Randomized Orthogonal Decomposition)

Let the snapshot matrix (16) containing the PDE solution states be recast into two time-shifted data matrices, whose columns are data snapshots:

$$\begin{array}{c}{V}_0=\left[\begin{array}{cccc}{u}_0& u_1& ...& u_{N_t-1}\end{array}\right]\equiv {\left\{u_i^0\right\}}_{i=0}^{N_t-1}\in {\mathbb{R}}^{N_x\times N_t},\hfill \\ V_1=\left[\begin{array}{cccc}{u}_1& u_2& ...& u_{N_t}\end{array}\right]\equiv {\left\{u_i^1\right\}}_{i=1}^{N_t}\in {\mathbb{R}}^{N_x\times N_t},\hfill \end{array}$$

(19)

For a target rank $k<N_t$, the state solution of PDE has the following reduced k-order representation:

$$u_i^{\left(k\right)}=u^{\left(k\right)}\left(x,t_i\right)=\sum _{j=1}^k{a}_j\left(t_i\right){\varphi }_j\left(x\right),t_i\in \left\{t_1,...,t_{N_t}\right\},i=1,...,N_t,$$

(20)

on the subspace spanned by the sequence of orthonormal Koopman modes:

$$span\left\{{\varphi }_1,{\varphi }_2,...,{\varphi }_k\right\},{\left\{{\varphi }_j\right\}}_{j=1}^k=〈U,X_{.,j}〉/∥〈U,X_{.,j}〉∥,$$

(21)

where U represents the matrix of left singular vectors produced by k-RSVD of data matrix $V_0$ and X denotes the eigenvectors to the Gram matrix of discrete Koopman propagator operator $\mathcal{K}$, i.e. $u_{N_t}={\mathcal{K}}^{N_t}{u}_0$.

The modal temporal coefficients are obtained as scalar products between data snapshots and Koopman modes:

$$a_j\left(t\right)=〈V_0,{\varphi }_j〉,j=1,...,k.$$

(22)

The proposed modal representation method is hereafter referred to as Koopman Randomized Orthogonal Decomposition (KROD).

Proof. 

Following the Koopman decomposition assumption [32], a Koopman propagator operator $\mathcal{K}$ exists, that maps every column vector onto the next one, i.e.

$$\left\{u_0,u_1=\mathcal{K}{u}_0,u_2=\mathcal{K}{u}_1={\mathcal{K}}^2{u}_0,...,u_{N_t}=\mathcal{K}{u}_{N_t-1}={\mathcal{K}}^{N_t}{u}_0\right\}.$$

(23)

The sequence (23) lies in the Krylov subspace

$${\mathcal{K}}_{N_t+1}\left(\mathcal{K},u_0\right)=span\left\{u_0,\mathcal{K}{u}_0,{\mathcal{K}}^2{u}_0,...,{\mathcal{K}}^{N_t}{u}_0\right\},$$

(24)

that captures how the observable $u_0$ evolves under the Koopman dynamics.

The Krylov subspace (24) gives a finite-dimensional proxy to study the action of the infinite-dimensional Koopman operator. For a sufficiently long sequence of the snapshots, suppose that the last snapshot $u_{N_t}$ can be written as a linear combination of previous $N_t$ vectors, on the Krylov subspace (24), such that:

$$u_{N_t}=c_0{u}_0+c_1{u}_1+...+c_{N_t-1}{u}_{N_t-1}+\mathcal{R},$$

(25)

in which ${\mathrm{c}}_i\in \mathbb{R},\mathrm{i}=0,...,{\mathrm{N}}_{\mathrm{t}}-1$ and $\mathcal{R}$ is the residual vector. The following relations are true:

$$\left\{u_1,u_2,...u_{N_t}\right\}=\mathcal{K}\left\{u_0,u_1,...u_{N_t-1}\right\}=\left\{u_1,u_2,...,V_{0}c\right\}+\mathcal{R},$$

(26)

where $c={\left(\begin{array}{cccc}{c}_0& c_1& ...& c_{N_t-1}\end{array}\right)}^T$ is the unknown column vector.

Assume a diagonal matrix $\mathcal{S}$ exists that asymptotically approximates the eigenvalues of $\mathcal{K}$ as ${∥\mathcal{R}∥}_2$ approaches zero. Solving the eigenvalue problem:

$$V_1=\mathcal{K}{V}_0=V_0\mathcal{S}+\mathcal{R},$$

(27)

is equivalent to solve the minimization problem:

$$\underset{\mathcal{S}}{min}\mathcal{R}={∥V_1-V_0\mathcal{S}∥}_2,$$

(28)

where ${∥·∥}_2$ is the $L_2$-norm of ${\mathbb{R}}^{N_x\times N_t}$.

For $1\le k\le \min \left(N_x,N_t\right)$, we identify the k-RSVD of $V_0$, that yields the factorization:

$$V_0=U\Sigma W^H,$$

(29)

where $U\in {\mathbb{R}}^{N_x\times k}$ and $W\in {\mathbb{R}}^{N_t\times k}$ are orthogonal matrices that contain the eigenvectors of ${\mathrm{V}}_0{{\mathrm{V}}_0}^H$ and ${{\mathrm{V}}_0}^H{\mathrm{V}}_0$, respectively, $\Sigma =diag\left({\sigma }_1,...,{\sigma }_k\right)\in {\mathbb{R}}^{k\times k}$ is a square diagonal matrix containing the singular values of $V_0$ and H means the conjugate transpose.

Relations $\mathcal{K}{V}_0=V_1=V_{0}S+\mathcal{R},{∥\mathcal{R}∥}_2\to 0$ and $V_0=U\Sigma W^H$ yield:

$$\mathcal{K}U\Sigma W^H=V_1=U\Sigma W^H\mathcal{S}\Rightarrow U^H\mathcal{K}U\Sigma W^H=U^{H}U\Sigma W^H\mathcal{S}\Rightarrow \mathcal{S}=U^H\mathcal{K}U.$$

As a consequence, the solution to the minimization problem (28) is the matrix operator:

$$\mathcal{S}=U^H\left(V_{1}W{\Sigma }^{-1}\right).$$

(30)

As a direct result, the eigenvalues of $\mathcal{S}$ will converge toward the eigenvalues of the Koopman propagator operator $\mathcal{K}$.

Let $X\in {\mathbb{C}}^{k\times k}$, $\Lambda \in {\mathbb{C}}^{k\times k}$ be the eigenvectors, respectively the eigenvalues of the Hermitian positive semi-definite Gram matrix $\mathcal{G}={\mathcal{S}}^H\mathcal{S}$:

$$\mathcal{G}X=X\Lambda .$$

(31)

From the k-RSVD of $V_0$, it follows that the column space of $V_0$ is the same as the column space of U, i.e.,

$$\mathrm{Range}\left(V_0\right)=\mathrm{Range}\left(U\right).$$

(32)

Using the definition $V_1=V_0\mathcal{S}$ and substituting into the expression for $\mathcal{S}$:

$$\begin{array}{cc}\hfill \mathcal{S}& =U^H{V}_{1}W{\Sigma }^{-1}\hfill \\ \hfill & =U^H\left(V_0\mathcal{S}\right)W{\Sigma }^{-1}\hfill \\ \hfill & =U^H\left(U\Sigma W^H\mathcal{S}\right)W{\Sigma }^{-1}\hfill \\ \hfill & =\left(U^{H}U\right)\Sigma \left(W^H\mathcal{S}W\right){\Sigma }^{-1}\hfill \\ \hfill & =\Sigma \left(W^H\mathcal{S}W\right){\Sigma }^{-1},\hfill \end{array}$$

since $U^{H}U=I$.

Since $\mathcal{S}$ is diagonal and W is unitary, $W^H\mathcal{S}W$ represents a unitary similarity transform [71] of matrix $\mathcal{S}$, i.e. $\mathcal{S}$ acts in a reduced coordinate system defined by the k-RSVD.

Let

$$\Phi =UX.$$

(33)

Since the columns of X are vectors in ${\mathbb{C}}^k$ and U is an orthonormal basis for the column space of $V_0$, the columns of $\Phi$ are linear combinations of the columns of U.

Thus,

$$\mathrm{Range}(\Phi )\subseteq \mathrm{Range}\left(U\right)=\mathrm{Range}\left(V_0\right),$$

(34)

i.e. the set of vectors $\Phi =UX$ lies in the column space of $V_0$. Therefore, the span of the eigenvectors $\Phi$ forms a subspace of the space containing the data in $V_0$.

It follows that there exists a subspace of the space containing the data in $V_0$, spanned by the sequence of Koopman modes:

$$span\left\{{\varphi }_1,{\varphi }_2,...,{\varphi }_k\right\},{\left\{{\varphi }_j\right\}}_{j=1}^k=〈U,X_{.,j}〉,$$

(35)

where U represents the matrix of left singular vectors produced by k-RSVD of data matrix $V_0$.

We proceed to demonstrate that the columns of $\Phi =UX$ are orthonormal. Since U comes from the k-RSVD of $V_0$, it has orthonormal columns:

$$U^{H}U=I.$$

(36)

Now compute the inner product of $\Phi$:

$${\Phi }^H\Phi ={\left(UX\right)}^H\left(UX\right)=X^H{U}^{H}UX=X^{H}X.$$

(37)

So the orthonormality of $\Phi$ depends on $X^{H}X$. The matrix X contains the eigenvectors of the Hermitian Gram matrix $\mathcal{G}={\mathcal{S}}^H\mathcal{S}$. Since $\mathcal{G}$ is Hermitian, it has a complete set of orthonormal eigenvectors. Thus, by choosing the eigenvectors to be orthonormal, we have:

$$X^{H}X=I\Rightarrow {\Phi }^H\Phi =X^{H}X=I,$$

(38)

i.e. columns of $\Phi$ are orthonormal.

It follows that:

$$〈{\varphi }_i,{\varphi }_j〉={\delta }_{ij},1\le i\le j\le k$$

(39)

i.e., Koopman modes $\Phi$ form an orthnormal base to the data space.

We continue by establishing the formula for the modal coefficients.

Assume that the data evolves linearly in reduced k-order representation (20), written in matrix formulation:

$$V_1=\Phi A.$$

(40)

Multiplying both sides on the left by ${\Phi }^H$ results:

$${\Phi }^H{V}_1={\Phi }^H\Phi A.$$

(41)

Since $\Phi$ has orthonormal columns, we get:

$${\Phi }^H\Phi =I\Rightarrow {\Phi }^H{V}_1=A.$$

(42)

Thus:

$$A={\Phi }^H{V}_1.$$

(43)

If $V_1=V_0\mathcal{S}$, and $\Phi =UX$, then we can alternatively write:

$$A={\Phi }^H{V}_1={\Phi }^H{V}_0\mathcal{S}.$$

(44)

We choose to project $V_0$ onto the subspace spanned by the orthonormal basis $\Phi$ to reduce computational complexity. Consequently, the coefficient matrix

$$A={\Phi }^H{V}_0$$

(45)

provides a representation with a smaller approximation error in the least-squares sense. □

Proposition 2. 

Given a set of integer target ranks $\mathcal{K}=\left\{1<k_1<k_2<k_3<...<k_N\right\},N\le N_t$ and N triples:

$$\mathcal{T}=\left\{T_1,T_2,...,T_N\right\},T_i=\left({\Phi }_i,A_i,k_i\right),i=1,...,N\le N_t$$

each triple encodes a Koopman modes basis, a modal amplitude set, and a corresponding target rank for the reduced k-order representation $u^{\left(T_i\right)}$ of the PDE solution, in form of Eq.(20).

Let $u^{\left(\mathcal{T}\right)}$ be the set of all N reduced-order solutions and u be the true PDE solutions.

Determination of the optimal triple to generate the data-twin model is achieved through a constrained multi-objective optimization framework designed to balance model fidelity with consistency to the reference dynamics:

$$\begin{array}{c}\mathrm{Given}:\mathcal{T}={\left\{T_i=\left({\Phi }_i,A_i,k_i\right)\right\}}_{i=1}^N\hfill \\ \mathrm{Find}:T^{*}\in \mathcal{T}\hfill \\ \mathrm{to}minE\left(\mathcal{T}\right)={∥u-u^{\left(\mathcal{T}\right)}∥}_2\hfill \\ \mathrm{and}maxC\left(\mathcal{T}\right)={\displaystyle \frac{{∥u^H{u}^{\left(\mathcal{T}\right)}∥}_2^2}{{∥u^{H}u∥}_2{∥u^{\left(\mathcal{T}\right)H}{u}^{\left(\mathcal{T}\right)}∥}_2}}\hfill \end{array}$$

(46)

where H denotes the conjugate transpose, and the notation ${∥·∥}_2$ typically denotes the Euclidean norm for vectors, or the Frobenius norm when the state solution of the partial differential equation resides in a space of dimension higher than two.

The optimal triple $T^{*}=\left({\Phi }_{DTM},A_{DTM},N_{DTM}\right)$ consists of the Koopman leading modes, the associated modal amplitudes and the rank of the Koopman basis, i.e. the number of terms in Koopman randomized orthogonal decomposition in form:

$$u_i^{DTM}=u^{DTM}\left(x,t_i\right)=\sum _{j=1}^{N_{DTM}}{a}_j\left(t_i\right){\varphi }_j\left(x\right),t_i\in \left\{{\mathrm{t}}_1,...,{\mathrm{t}}_{N_t}\right\},i=1,...,N_t.$$

(47)

The optimization problem (46) is resolved numerically, with the methodology and implementation details provided in a subsequent section.

3.2. Online Phase: Modeling the Data-Twin Temporal Dynamic by Explainable Deep Learning

The calculation step of the modes in the offline phase is performed only once. The data-twin model is not complete without the temporal component. These computations are carried out in the second step, during which the mathematical model is integrated over time. Additionally, if needed, the dynamics of the reduced-order model can be evaluated at any specific time point. For this reason, this step is referred to as the online phase, as it can be executed on demand to simulate the behavior of the original process.

Unlike conventional machine learning models, which typically function as opaque black boxes operating on entire datasets without disentangling underlying factors, a key innovation of this study lies in the explicit modeling of temporal dynamics decoupled from spatial dependencies. This is achieved through the use of explainable deep learning techniques, which offer interpretable insights into the data-twin model’s parameters.

In this work, Nonlinear AutoRegressive models with eXogenous inputs (NLARX) are employed as a novel deep learning-based approach within the domain of nonlinear system identification. Since the emergence of artificial neural networks as numerical tools [72,73], NLARX models have been used for various purposes, ranging from simulation [74], to nonlinear predictive control [75] or higher order nonlinear optimization problems [76,77].

The objective is to identify a dynamical system that effectively captures the nonlinear temporal evolution of the underlying process, expressed in the form:

$${\displaystyle \frac{d{\widehat{a}}_j\left(t\right)}{dt}}=\mathcal{M}\left({\widehat{a}}_j\left(t\right),a_j\left(t\right),t\right),a_j\left(t_0\right)=a_j^0,t\in {\mathbb{R}}_{\ge 0},$$

(48)

based on modal temporal coefficients $a_j\left(t_i\right),j=1,...,N_{DTM},t_i\in \left\{{\mathrm{t}}_1,...,{\mathrm{t}}_{N_t}\right\},$ computed by KROD algorithm in the offline phase.

To approximate the system dynamics, the evolution of the state $\widehat{a}\left(t\right)$ is modeled using the NLARX model structure, which is expressed as:

$$\widehat{a}\left(t\right)=\mathcal{F}\left[\widehat{a}\left(t-1\right),...,\widehat{a}\left(t-n_a\right),a\left(t-n_k\right),...,a\left(t-n_k-n_b+1\right)\right]+e\left(t\right),$$

(49)

where $\mathcal{F}$ is a nonlinear estimator realized by a cascade forward neural network, capturing the underlying nonlinear temporal dependencies and input effects. Here, $n_a$ denotes the number of past output terms, $n_b$ represents the number of past input terms used to predict the current output, and $n_k$ corresponds to the input delay, while $e\left(t\right)$ captures the modeling error. Each output of the NLARX model described in (49) is determined by regressors, which are transformations of previous input and output values. Typically, this mapping consists of a linear component and a nonlinear component, with the final model output being the sum of these two contributions.

Thus, the NLARX model serves as an empirical approximation of the continuous dynamics $\mathcal{M}$, with the neural network mapping $\mathcal{F}$ implicitly representing the evolution operator of the system, such that:

$$\widehat{a}\left(t\right)\approx \widehat{a}\left(t-1\right)+\Delta t·\mathcal{M}\left(\widehat{a}\left(t-1\right),a\left(t-1\right),t-1\right)\approx \mathcal{F}\left[\widehat{a}\left(t-1\right),...,a\left(·\right)\right]+e\left(t\right),$$

(50)

where $\Delta t$ is the discrete time step and $e\left(t\right)$ denotes the modeling error.

Training the NLARX model can be formulated as a nonlinear unconstrained optimization problem, of form:

$$\theta \left(n_a,n_b,n_k\right)=argmin{\displaystyle \frac{1}{2N_t}}\sum _{i=1}^{N_t}{∥a\left(t_i\right)-\widehat{a}\left(t_i\right)∥}_2.$$

(51)

In this formulation, the training dataset comprises the measured input $a\left(t\right)$, while $\widehat{a}\left(t\right)$ denotes the output generated by the NLARX model. The symbol ${∥·∥}_2$ refers to the $L_2$ norm and $\theta \left(n_a,n_b,n_k\right)$ represents the parameter vector associated with the nonlinear function $\mathcal{F}$.

The NLARX architecture is well-suited for capturing system dynamics by feeding previous network outputs back into the input layer. It also allows for the specification of how many past input and output time steps are necessary to accurately represent the system’s behavior. A critical aspect of effectively applying an NLARX network lies in the careful selection of inputs, input delays, and output delays. The objective is to adjust the network parameters $\theta \left(n_a,n_b,n_k\right)$ across the entire trajectory to minimize the objective function defined in (51).

The reduced-order data-driven twin model for the PDE solution is obtained in the form:

$$\left\{\begin{array}{c}{u}^{DTM}\left(x,t\right)={\displaystyle \sum _{j=1}^{N_{DTM}}}{\widehat{a}}_j\left(t_i\right){\varphi }_j\left(x\right),t\in R_{\ge 0},\hfill \\ \widehat{a}\left(t\right)=\mathcal{F}\left[\widehat{a}\left(t-1\right),...,\widehat{a}\left(t-n_a\right),a\left(t-n_k\right),...,a\left(t-n_k-n_b+1\right)\right],a_j\left(t_0\right)=a_j^0.\hfill \end{array}\right.$$

(52)

3.3. Qualitative Analysis of the Data-Driven Twin Model

A qualitative assessment of the orthogonality between the Koopman modes derived using the KROD algorithm proposed in this study can be carried out through the Modal Assurance Criterion (MAC) [78]. The MAC value corresponding to a pair of shape modes is defined as follows:

$$MAC\left({\varphi }_i,{\varphi }_j\right)={\displaystyle \frac{{〈{\varphi }_i,{\varphi }_j〉}^2}{〈{\varphi }_i,{\varphi }_i〉·〈{\varphi }_j,{\varphi }_j〉}}.$$

(53)

The $MA{C}_{ij}$ values lie within the interval $\left[0,1\right]$ as the leading koopman modes obtained from the KROD algorithm are normalized vectors. In this context, Koopman modes are considered orthogonal if:

$$MA{C}_{ij}=\left\{\begin{array}{c}1,i=j\hfill \\ 0,i\ne j.\hfill \end{array}\right.$$

(54)

The Pearson correlation coefficient is employed to quantify the correlation between the data-twin model solution and the true solution of the partial differential equation, providing a measure of how closely the data-twin model approximates the true dynamics across the space-time domain.

Let $u(x,t)$ be the true solution and $u^{DTM}\left(x,t\right)$ be the data-twin model solution of PDE, sampled at $N_x$ spatial points and $N_t$ time instances. Define the flattened vectors:

$$\mathbf{U}=[u(x_1,t_1),\dots ,u(x_{N_x},t_1),u(x_1,t_2),\dots ,u(x_M,t_{N_t})]\in {\mathbb{R}}^{N_x{N}_t},$$

$${\mathbf{U}}^{\mathbf{DTM}}=[u^{DTM}(x_1,t_1),\dots ,u^{DTM}(x_{N_x},t_1),u^{DTM}(x_1,t_2),\dots ,u^{DTM}(x_M,t_{N_t})]\in {\mathbb{R}}^{N_x{N}_t}.$$

The Pearson correlation coefficient $\rho$ between $\mathbf{U}$ and ${\mathbf{U}}_{\mathbf{DTM}}$ is given by:

$$\rho ={\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_x{N}_t}}\left(U_k-\overline{U}\right)\left(U_k^{DTM}-\overline{U^{DTM}}\right)}{\sqrt{{\displaystyle \sum _{k=1}^{N_x{N}_t}}{\left(U_k-\overline{U}\right)}^2}\sqrt{{\displaystyle \sum _{k=1}^{N_x{N}_t}}{\left(U_k^{DTM}-\overline{U^{DTM}}\right)}^2}}},$$

(55)

where

$$\overline{U}={\displaystyle \frac{1}{N_x{N}_t}}\sum _{k=1}^{N_x{N}_t}{U}_k,\overline{U^{DTM}}={\displaystyle \frac{1}{N_x{N}_t}}\sum _{k=1}^{N_x{N}_t}{U}_k^{DTM}.$$

In addition to the Pearson correlation coefficient, the mean absolute error (MAE) is employed as a complementary metric to assess the average magnitude of pointwise discrepancies over the space-time domain, providing a more direct measure of predictive accuracy irrespective of correlation structure:

$$MAE={\displaystyle \frac{1}{N_x{N}_t}}\sum _{k=1}^{N_x{N}_t}\left|U_k-U_k^{DTM}\right|.$$

(56)

In addition to these metrics, to provide fine-grained, pointwise insight into how well the data-twin model approximates the true solution of PDE at specific locations in space and time, we define the absolute local error matrix $E\in {\mathbb{R}}^{N_x\times N_t}$ as:

$$E_{ij}=\left|u\left(x_i,t_j\right)-u^{DTM}\left(x_i,t_j\right)\right|$$

(57)

for all $i=1,...,N_x$ and $j=1,...,N_t$.

3.4. Time Simulation and Validation of the Data-Driven Twin Model

Data-driven twin model is subjected to time-domain simulation, followed by a two-step input-output validation of its predicted response. In the first validation step, the model is trained using the initial two-thirds of the available snapshots, with the remaining one-third reserved for validation. In the subsequent step, the training dataset is further reduced to the first one-third of the snapshots, and the remaining two-thirds are employed to assess the model’s predictive performance.

4. Data-Driven Twin Modeling of Shock Wave Phenomena Using KROD

This section demonstrates the application of the proposed KROD algorithm to shock wave modeling using the viscous Burgers’ equation.

4.1. Governing Equations of the Mathematical Model

The viscous Burgers equation model is considered, of the form:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial }{\partial t}}u\left(x,t\right)+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{u{\left(x,t\right)}^2}{2}}\right)=\nu {\displaystyle \frac{{\partial }^2}{\partial x^2}}u\left(x,t\right),t>0,\nu >0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),x\in \mathbb{R},\hfill \end{array}\right.$$

(58)

where $u\left(x,t\right)$ is the unknown function of time t, $\nu$ is the viscosity parameter.

Three experiments are considered, exhibiting progressively more complex dynamics, each characterized by an initial condition of the following form:

$$Experiment1:u_0\left(x\right)=-sin\left(\pi x\right),$$

(59)

$$Experiment2:u_0\left(x\right)=\left\{\begin{array}{c}{u}_L,x\le 0\hfill \\ u_R,x>0\hfill \end{array}\right.,$$

(60)

$$Experiment3:u_0\left(x\right)=-cos{\left(1.5\pi x\right)}^2.$$

(61)

Homogeneous Dirichlet boundary conditions are imposed for all three examples, specified as follows:

$$u\left(0,t\right)=u\left(L,t\right)=0$$

(62)

Equation (59), which defines the initial condition for Experiment 1, produces a sinusoidal pulse characterized by an abrupt change in slope at the boundaries of the domain. The initial value problem with the discontinuous initial condition given by (60) is referred to as a Riemann problem, which leads to the formation of a shock wave. For Experiment 3, the initial condition described by (61) induces a more complex evolution of the pulse. All three experiments involve phenomena that pose significant challenges for standard numerical methods to accurately capture.

The Cole-Hopf transformation, independently introduced by Hopf (1950) [79] and Cole (1951) [80], converts the Burgers’ equation into the heat equation, thereby enabling the derivation of its exact solution. The Cole-Hopf transformation is presented in the following subsection.

4.2. Derivation of the Analytical Solution Using the Cole-Hopf Transformation

The Cole-Hopf transformation is defined by

$$u=-2\nu {\displaystyle \frac{1}{\phi }}{\displaystyle \frac{\partial \phi }{\partial x}}.$$

(63)

Through analytical treatment, it is found that:

$${\displaystyle \frac{\partial u}{\partial t}}={\displaystyle \frac{2\nu }{{\phi }^2}}\left({\displaystyle \frac{\partial \phi }{\partial t}}{\displaystyle \frac{\partial \phi }{\partial x}}-\phi {\displaystyle \frac{{\partial }^2\phi }{\partial x\partial t}}\right),u{\displaystyle \frac{\partial u}{\partial x}}={\displaystyle \frac{4{\nu }^2}{{\phi }^3}}{\displaystyle \frac{\partial \phi }{\partial x}}\left(\phi {\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}-{\displaystyle \frac{\partial \phi }{\partial x}}{\displaystyle \frac{\partial \phi }{\partial x}}\right),$$

(64)

$$\nu {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}=-{\displaystyle \frac{2{\nu }^2}{{\phi }^3}}\left(2{\left({\displaystyle \frac{\partial \phi }{\partial x}}\right)}^3-3\phi {\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}{\displaystyle \frac{\partial \phi }{\partial x}}+{\phi }^2{\displaystyle \frac{{\partial }^3\phi }{\partial x^3}}\right).$$

(65)

Substituting these expressions into (58) it follows that:

$${\displaystyle \frac{\partial \phi }{\partial x}}\left({\displaystyle \frac{\partial \phi }{\partial t}}-\nu {\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}\right)=\phi \left({\displaystyle \frac{{\partial }^2\phi }{\partial x\partial t}}-\nu {\displaystyle \frac{{\partial }^3\phi }{\partial x^3}}\right)=\phi {\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial \phi }{\partial t}}-\nu {\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}\right).$$

(66)

Relation (66) indicates that if $\phi$ solves the heat equation, then $u\left(x,t\right)$ given by the Cole-Hopf transformation (Section 4.2) solves the viscous Burgers equation (58). Thus the viscous Burgers equation (58) is recast into the following one:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial \phi }{\partial t}}-\nu {\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}=0,x\in R,t>0,\nu >0,\hfill \\ \phi \left(x,0\right)={\phi }_0\left(x\right)=e^{-{\int }_0^x{\displaystyle \frac{u_0\left(\xi \right)}{2\nu }}d\xi },x\in \mathbb{R}.\hfill \end{array}\right.$$

(67)

Taking the Fourier transform with respect to x for both heat equation and the initial condition (67) it is found that:

$$\phi \left(x,t\right)={\displaystyle \frac{1}{2\sqrt{\pi \nu t}}}\underset{-\infty }{\overset{\infty }{\int }}{\phi }_0\left(\xi \right)e^{-{\displaystyle \frac{{(x-\xi )}^2}{4\nu t}}}d\xi .$$

(68)

Using the Cole-Hopf transformation (Section 4.2), the analytic solution to the problem (58) is obtained in the following form:

$$u\left(x,t\right)={\displaystyle \frac{{\int }_{-\infty }^{\infty }\frac{x-\xi }{t}{\phi }_0\left(\xi \right)e^{-\frac{{\left(x-\xi \right)}^2}{4\nu t}}d\xi }{{\int }_{-\infty }^{\infty }{\phi }_0\left(\xi \right)e^{-\frac{{\left(x-\xi \right)}^2}{4\nu t}}d\xi }}.$$

(69)

4.3. Derivation of the Numerical Exact Solution Using the Gauss-Hermite Quadrature

To compute the exact solution (69) corresponding to the initial conditions of the three test cases (59)–(61), the Gauss-Hermite quadrature technique [81] is employed.

Gauss–Hermite quadrature approximates the value of integrals of the following kind:

$$\underset{-\infty }{\overset{\infty }{\int }}f\left(z\right)e^{-z^2}dz\approx \sum _{i=1}^n{w}_{i}f\left(x_i\right),$$

(70)

where n represents the number of sample points used, $x_i$ are the roots of the Hermite polynomial $H_n\left(x\right)$ and the associated weights $w_i$ are given by

$$w_i={\displaystyle \frac{2^{n-1}n!\sqrt{\pi }}{n^2{\left(H_{n-1}\left(x_i\right)\right)}^2}},i=1,...,n.$$

(71)

In Experiment 1, the initial condition is defined by Eq.(59):

$${\phi }_0\left(x\right)=e^{-{\displaystyle \frac{1}{2\nu }}{\int }_0^x{u}_0\left(\xi \right)d\xi }=e^{-{\displaystyle \frac{1}{2\nu }}{\int }_0^x-sin\left(\pi \xi \right)d\xi }=e^{{\displaystyle \frac{1}{2\nu \pi }}}·e^{-{\displaystyle \frac{cos\left(\pi x\right)}{2\nu \pi }}}.$$

(72)

and the exact solution (69) is in form:

$$u\left(x,t\right)={\displaystyle \frac{{\int }_{-\infty }^{\infty }\frac{x-\xi }{t}{e}^{-\frac{cos\left(\pi \xi \right)}{2\nu \pi }}·e^{-{\left(\frac{x-\xi }{\sqrt{4\nu t}}\right)}^2}d\xi }{{\int }_{-\infty }^{\infty }{e}^{-\frac{cos\left(\pi \xi \right)}{2\nu \pi }}·e^{-{\left(\frac{x-\xi }{\sqrt{4\nu t}}\right)}^2}d\xi }}.$$

(73)

With the variable change:

$$z={\displaystyle \frac{x-\xi }{\sqrt{4\nu t}}},$$

(74)

the exact solution to the Burgers equation model (58) with the initial condition (59) is equivalent with the following form:

$$u\left(x,t\right)={\displaystyle \frac{{\int }_{-\infty }^{\infty }4\nu z{e}^{-\frac{1}{2\nu \pi }cos\left[\pi \left(x-z\sqrt{4\nu t}\right)\right]}{e}^{-z^2}dz}{{\int }_{-\infty }^{\infty }\sqrt{4\nu t}{e}^{-\frac{1}{2\nu \pi }cos\left[\pi \left(x-z\sqrt{4\nu t}\right)\right]}{e}^{-z^2}dz}}.$$

(75)

Similarly, by applying the variable transformation (74), the exact solutions for the subsequent two experiments are derived as follows.

In the case of Experiment 2 the exact solution to the Burgers equation model (58) with the initial condition (60) is:

$$u\left(x,t\right)={\displaystyle \frac{{\int }_{-\infty }^{\infty }4\nu z{e}^{-\frac{u_R}{2\nu }\left(x-z\sqrt{4\nu t}\right)}{e}^{-z^2}dz}{{\int }_{-\infty }^{\infty }\sqrt{4\nu t}{e}^{-\frac{u_R}{2\nu }\left(x-z\sqrt{4\nu t}\right)}{e}^{-z^2}dz}}.$$

(76)

In the case of Experiment 3 the exact solution to the Burgers equation model (58) with the initial condition (61) is:

$$u\left(x,t\right)={\displaystyle \frac{{\int }_{-\infty }^{\infty }4\nu z{e}^{\frac{1}{4\nu }\left(\frac{1}{3\pi }sin\left(3\pi (x-z\sqrt{4\nu t})\right)+x-z\sqrt{4\nu t}\right)}{e}^{-z^2}dz}{{\int }_{-\infty }^{\infty }\sqrt{4\nu t}{e}^{\frac{1}{4\nu }\left(\frac{1}{3\pi }sin\left(3\pi (x-z\sqrt{4\nu t})\right)+x-z\sqrt{4\nu t}\right)}{e}^{-z^2}dz}}.$$

(77)

To obtain the numerical exact solutions for the three test cases, both the numerator and denominator are evaluated using Gauss-Hermite quadrature (70).

5. Numerical Results

The following section presents numerical results demonstrating the computational performance of the Koopman Randomized Orthogonal Decomposition (KROD) method introduced in this study, in conjunction with Deep Learning (DL) for time simulation. It is shown that the resulting output is a high-fidelity, reduced-complexity model that effectively replicates the behavior of the original system and enables real-time simulation with high accuracy, thereby producing a data-driven twin model.

The computational domain is defined as $\left[0,L\right]$, with $L=2$, and the time interval considered is $\left[0,T\right]$, where $T=3$. The viscosity parameter in the Burgers equation is set to $\nu ={10}^{-2}$ for all three experiments. The domain is uniformly discretized using $N=100$ grid points, resulting in a mesh size of $\Delta x=0.02$.

The Riemann problem parameters are assumed to be $u_L=0.1$ and $u_R=0.5$.

By applying the Cole-Hopf transformation (Section 4.2), the analytical solutions to the Burgers problem (58) are obtained for all three experiments, corresponding respectively to the forms (75), (76), and (77), as detailed in Section 4.2. Subsequently, the exact solutions for the three test cases are computed numerically using the Gauss-Hermite quadrature method (70) with $n=100$ nodes (see Section 4.3).

Figure 1 presents the exact solutions for the three experiments under consideration.

The data-driven twin models of the form (47) are generated using the Koopman Randomized Orthogonal Decomposition (KROD) algorithm as high-fidelity surrogates for the three Burgers problems examined in this study. For this purpose, the training dataset consists of a total of $N_t=300$ snapshots collected at uniformly spaced time intervals of $\Delta t=0.01$, and $N_x=101$ spatial measurements recorded for each time snapshot.

KROD algorithm determines the optimal dimension $N_{DTM}$ of the subspace spanned by the dominant Koopman modes by solving the multiobjective optimization problem with nonlinear constraints stated in (46). This problem is solved by employing a genetic algorithm to identify the Pareto front corresponding to the two fitness functions.

Figure 2 depicts the optimization objectives defined in (46) along with the Pareto front solutions for each of the three experiments.

For Experiment 1, the dimension of the subspace spanned by the leading modes is $N_{DTM}=10$. In Experiment 2, the data-twin model representing the Riemann problem dynamics has a subspace dimension of $N_{DTM}=20$. For Experiment 3, the increased complexity of the pulse is captured by a data-twin model with a higher dimension of $N_{DTM}=45$. In all three cases, the computational time required by the KROD algorithm in offline phase remains extremely low, as illustrated in Figure 3.

KROD algorithm identifies the optimal dimension of the modal basis and computes the optimal Koopman modes along with their associated time-dependent amplitudes. As a result, the data-driven twin model (47) is constructed, which-by definition-exhibits minimal approximation error and maximal correlation with respect to the original data model.

To facilitate a qualitative assessment of the DTMs, Table 1 reports the outcomes of the multiobjective optimization problem (46) as solved by the KROD. Specifically, the table includes the dimensionality of the leading Koopman modes, the absolute error as defined in Eq. (56), and the correlation coefficient as given by Eq. (55) between the exact solution and the corresponding data-twin model for each of the three experiments considered.

An analysis of the results presented in Table 1 reveals that KROD algorithm produces models exhibiting perfect correlation with the original data, i.e., $\rho =1$, across all three test cases. Furthermore, the absolute errors remain remarkably low, on the order of $MAE<\mathcal{O}\left({10}^{-7}\right)$.

The selection of modes solely based on their energy ranking proves effective only under specific conditions [47,82]. Some modes, despite contributing minimally to the overall energy, may exhibit rapid temporal growth at low amplitudes or represent strongly damped high-amplitude structures. Such modes can play a significant role in enhancing the fidelity of the data-twin model. A key advantage of the KROD algorithm over traditional DMD-based approaches lies in its ability to generate a substantially lower-dimensional subspace that retains the most dynamically significant Koopman modes. As a result, the KROD method avoids the need for an additional-often subjective-mode selection criterion. Moreover, the orthogonality of the KROD modes contributes to their qualitative superiority, enabling accurate modeling with a reduced number of modes.

To assess the orthogonality of the Koopman modes computed by KROD algorithm, the orthogonality matrix introduced in Section 3.3 is evaluated for each of the three experiments considered in this study. The following figures display the computed $\mathrm{MAC}$ values for each pair of leading shape modes, as defined previously in Eq. (53).

Figure 4. The orthogonality matrices corresponding to the three investigated experiments are shown, respectively. The $\mathrm{MAC}$ values, computed according to Eq. (53), confirm the mutual orthogonality of the Koopman modes in each case.

Figure 4. The orthogonality matrices corresponding to the three investigated experiments are shown, respectively. The $\mathrm{MAC}$ values, computed according to Eq. (53), confirm the mutual orthogonality of the Koopman modes in each case.

The KROD algorithm facilitates the identification of the leading Koopman modes and their corresponding temporal coefficients. To achieve high-fidelity time simulation of the data-driven twin model’s temporal parameters, Nonlinear AutoRegressive models with eXogenous inputs (NLARX) are employed (revise Subsection (Section 3.4). Figure 5, Figure 6 and Figure 7 present the two-fold input-output validation results for three randomly selected temporal coefficients, simulated using the optimal NLARX models for Experiments 1, 2, and 3, respectively. These results demonstrate that the deep learning-based NLARX models accurately replicate the coefficient dynamics, even when the training dataset is reduced to only the first third of the available snapshots. Relevant performance metrics for the deep learning simulations are also provided in the figures.

The time-simulated data-driven twin models are assembled in the form specified by Eq. (52) for each of the three experiments considered. The output responses are then compared to the corresponding exact solutions in terms of absolute local error given by Eq. (57), as illustrated in Figure 8.

Figure 9 presents the CPU time required during the online phase of the algorithm, in which the deep learning process is executed, as a function of the optimized modal space dimension for the three experimental cases analyzed. The reported CPU times correspond to two successive computational runs. All simulations and visualizations were performed using custom scripts developed in Matlab R2022a and executed on a system equipped with an Intel i7-7700K processor. The observed computational times are considered acceptable and fall within a reasonable range for practical implementation.

6. Conclusions

This study introduces a novel framework for constructing data-driven twin models of reduced computational complexity for partial differential equations models, a topic of significant relevance in data science. To our knowledge, this work is the first to introduce the Koopman Randomize Orthogonal Decomposition (KROD) algorithm within the context of nonlinear PDEs.

The proposed approach consists of two main phases. Initially, a new numerical algorithm is developed to generate a reduced-order data-twin model based on the KROD decomposition of the data. Subsequently, explainable NLARX Deep Learning technique is employed to enable real-time, adaptive calibration of the data-twin model and to compute its dynamic response.

The resulting twin data model exhibits minimal error and maximal correlation with the original dataset, offering a highly accurate and computationally efficient representation of the underlying physical system.

This study introduces several methodological advancements that distinguish it from classical modal decomposition techniques such as Dynamic Mode Decomposition and Proper Orthogonal Decomposition, which often depend on heuristic or ad hoc criteria for mode selection. In contrast, the proposed Randomized Koopman Orthogonal Decomposition leverages randomized singular value decomposition to substantially reduce computational cost while fully automating the selection of dominant modes.

Building upon the author’s earlier work [55,56,70], it is shown that the stochastic nature of RSVD inherently favors the most relevant modes, thereby enhancing robustness without introducing additional algorithmic complexity.

A key novelty lies in the first quantitative framework for computing orthogonal Koopman modes using randomized orthogonal projections, which ensures optimal representation in the data space. Moreover, the proposed method integrates explainable deep learning techniques within the modal decomposition process, enabling interpretable and real-time calibration of the reduced-order data-twin model. This integration also provides a pathway for future extensions into unsupervised learning paradigms.

Importantly, the framework is inherently compatible with unsupervised learning, functioning effectively on raw observational data without requiring predefined models. This characteristic promotes broader applicability and generalization across a wide range of fluid dynamics problems and complex dynamical systems. The specific application of the KROD methodology to fluid dynamics will be addressed in future investigations.