High-fidelity Digital Twin Data Models by Randomized Dynamic Mode Decomposition and Deep Learning with Applications in Fluid Dynamics Diana A. Bistrian
arXiv:2609.17101v1 [cs.LG] 15 Sep 2026
University Politehnica Timisoara, Department of Electrical Engineering and Industrial Informatics, Romania [email protected]
Originally published in Modelling 3 (2022), no. 3, 314–332. Received 9 June 2022; accepted 19 July 2022; published 21 July 2022. https://doi.org/10.3390/modelling3030020 Abstract The purpose of this paper is the identification of high-fidelity digital twin data models from numerical code outputs by non-intrusive techniques (i.e. not requiring Galerkin projection of the governing equations onto the reduced modes basis). In this paper we define the concept of digital twin data model (DTM) as a model of reduced complexity that has the main feature to mirror the original process behavior. The significant advantage of a DTM is to map the dynamics with high accuracy and reduced costs in CPU time and hardware, to timescales difficult to explore because of the complexity of the dynamics over time. This paper introduces a new framework for creating efficient digital twin data models by combining two state-of-the-art tools: randomized dynamic mode decomposition and deep learning. We show that the outputs are consistent with the original source data with the advantage of a reduced complexity. The DTMs are investigated in the numerical simulation of three shock wave phenomena with increasing complexity. We perform a thorough assessment of the performance of the new digital twin data models in terms of numerical accuracy and computational efficiency.
Keywords: digital twin data model; randomized dynamic mode decomposition; shock wave phenomena; deep learning
1
Introduction
Data-driven algorithms for analyzing complex systems are of growing interest. Among fluid dynamics researchers, many efforts were directed in recent years to identification of a reliable approximation of the complex flow dynamics by models of low complexity, i.e. reduced order models (ROM). In order to reveal underlying physical processes in an appropriate way, a dynamic analysis should be made. Modal decomposition techniques are superior to other techniques, because they are able to link physical behaviour to a spatial pattern. Among several modal decomposition methods, Proper Orthogonal Decomposition (POD) [1–6] and Dynamic Mode Decomposition (DMD) [7–13] have been widely applied to study the physics of the dynamics of the flows in different applications. The POD and its variants are also known as Karhunen-Loeve expansions in feature selection and signal processing, empirical orthogonal functions in atmospheric science or principal component analysis in statistics. The strong point of POD is that it can be applied to non-linear partial differential equations, especially for smooth systems in which the energetics can be hierarchically
1
ranged and characterized by the first few modes. The applicability of POD to complex systems is limited mainly due to errors associated with the truncation of the POD modes [14–16]. Being rooted in the Koopman mode theory [17], a recent decomposition technique, namely Dynamic Mode Decomposition (DMD) [7, 8, 18], has the significant advantage of linking a spatial structure (coherent structure) to a single oscillating frequency and growing/decay rate. Therefore, DMD is promising, especially for hydrodynamic research like flow field analysis. Since the first application of the Koopman theory for the purposes of reduced order modelling by Igor Mezić [7], a considerable amount of work has focused on understanding and improving the method of dynamic mode decomposition and several DMD procedures have been released: optimized DMD [19], exact DMD [20], sparsity promoting DMD [21], multi-resolution DMD [22, 23], extended DMD [24], recursive DMD [25, 26], DMD with control [27], randomized DMD [28, 29], dynamic mode decomposition with core sketch [30]. Under assumptions on the underlying dynamics, Mezić [31] provided the first result on the convergence rate under sample size increase in the case of finite-section approximation and introduced a discussion on the choice of observables in the context of finite-section approximations. A comparison of DMD vs. POD for model reduction was illustrated in our previous paper [16], for the study of shallow water equations model. A procedure of coupling POD and DMD for nonlinear model order reduction was introduced in [32]. Modelling fluid dynamics data is even more difficult when we handle with discrete or so called non-intrusive data, especially when there is no mathematical model associated with the data. In this paper we introduce the concept of digital twin data model (DTM) as a model of reduced complexity that has the main feature to mirror the original process behavior. The significant advantage of a DTM is to map the dynamics with high accuracy and reduced costs in CPU time and hardware, even to timescales difficult to explore because of the rapidly changing dynamics over time. This paper introduces a new framework for creating efficient digital twin data models from non-intrusive data by combining two state-of-the-art tools: randomized dynamic mode decomposition introduced in [28] and deep learning artificial intelligence. We show that the outputs are consistent with the original source data with the advantage of a reduced complexity. The DTMs are investigated in the numerical simulation of three shock wave phenomena with increasing complexity. We perform a thorough assessment of the performance of the new digital twin data models in terms of numerical accuracy and computational efficiency. The remainder of the article is organized as follows. In Section 2 the test problem consisting of the nonlinear viscous Burgers equation model is presented. In Section 3 we recall the principles governing the dynamic mode decomposition and we provide the description of the randomized dynamic mode decomposition algorithm in Section 4. In Section 5 we outline the technique of fast digital twin data model identification using deep learning Nonlinear Autoregressive Estimators. Section 6 presents the numerical results together with a computational efficiency study. Summary and conclusions are drawn in the final section.
2
Shock Wave Phenomena: Full-order Model of Nonlinear Viscous Burgers Equation
We consider that the experimental data are provided by the simulation of the nonlinear viscous Burgers equation model [33]) (
∂ ∂ ∂t u (x, t) + ∂x
u(x,t)2 2
2
∂ = ν ∂x 2 u (x, t) ,
u (x, 0) = u0 (x) ,
2
x ∈ [0, L] , t ∈ [0, T ] ,
(1)
where u (x, t) is the unknown function of time t, ν = 1/Re is the viscosity term and Re is the Reynolds number. We consider the discontinuous initial condition of the form (
u0 (x) =
uL , uR ,
x ≤ 0, x > 0.
(2)
This setting yields a shock wave phenomenon. The the initial value problem (1)-(2) is solved using a finite difference discretization of the conservative form of the equation (1), and then carrying out a parabolic integration scheme [34]. The constants used for the test model are L = 2,
T = 3,
uL = 0.1,
uR = 0.5.
The training data comprises of Nt + 1 = 301 total number of snapshots taken in time at regularly spaced time intervals ∆t = 0.01, Nx = 101 number of spatial measurements per time snapshot. The nonlinear evolution governed by Burgers equation can be obtained also by the Cole-Hopf transformation [35]. The Cole-Hopf transformation is defined by u = −2ν
1 ∂φ . φ ∂x
(3)
Through an analytical handling we find that ∂u 2ν = 2 ∂t φ
!
!
∂φ ∂φ ∂2φ −φ , ∂t ∂x ∂x∂t
∂2u 2ν 2 ν 2 =− 3 ∂x φ
∂φ 2 ∂x
∂u 4ν 2 ∂φ ∂ 2 φ ∂φ ∂φ u = 3 φ 2 − , ∂x φ ∂x ∂x ∂x ∂x 3
(4)
!
∂ 2 φ ∂φ ∂3φ − 3φ 2 + φ2 3 . ∂x ∂x ∂x
(5)
Substituting these expressions into (1) it follows that ∂φ ∂x
∂2φ ∂φ −ν 2 ∂t ∂x
!
∂2φ ∂3φ =φ −ν 3 ∂x∂t ∂x
!
∂ =φ ∂x
!
∂2φ ∂φ −ν 2 . ∂t ∂x
(6)
Relation (6) indicates that if φ solves the heat equation, then u (x, t) given by the Cole-Hopf transformation (3) solves the viscid Burgers equation (1). Thus we have reduced the viscid Burgers equation (1) to the following one ∂φ − ν ∂ 2 φ = 0,
x ∈ R, t > 0, ν > 0,
∂t ∂x2 R x u (ξ) φ (x, 0) = φ (x) = e− 0 02ν dξ , x ∈ R.
(7)
0
Taking the Fourier transform with respect to x for both heat equation and the initial condition (7) we obtain the analytic solution 1 φ (x, t) = √ 2 πνt
Z∞
(x−ξ)2
φ0 (ξ) e− 4νt dξ.
(8)
−∞
From the Cole-Hopf transformation (3) we obtain the analytic solution to the problem (1) in the following form (x−ξ)2 R ∞ x−ξ − 4νt dξ −∞ t φ0 (ξ) e u (x, t) = R . (x−ξ)2 ∞ −
−∞ φ0 (ξ) e
3
4νt
dξ
(9)
Three solution types will be discussed in the paper, corresponding to Reynolds number of Re = 102 , Re = 103 and Re = 104 , respectively (see Figure 1). We notice that the model solution exhibits some oscillations for all three experiments. The unphysical oscillation originates due to high Reynolds numbers. We experienced that when the Reynolds number gets higher, the numerical solution exhibits more oscillations, so that the fluid dynamics become more and more complex.
Figure 1: Dynamics of shock wave phenomena as the exact solution of viscous Burgers equation model at Re = 102 , Re = 103 and Re = 104 , respectively. We aim in this paper to identify a reduced order model of the nonlinear viscous Burgers equation model to approximate as faithful the true solution and to create a digital twin data model of low complexity for the three shock wave phenomena, respectively. An efficient numerical technique is provided in the following sections.
3
Reduced order modeling based on Dynamic Mode Decomposition
3.1
The key steps of Dynamic Mode Decomposition and snapshots collection
We proceed by collecting data ui (t, x) = u (ti , x) , ti = i∆t, i = 0, ..., Nt , at the constant sampling time ∆t, x representing the spatial coordinate. We form a data matrix whose columns represent the individual data samples, called the snapshot matrix: h i V = u0 u1 ... uNt ∈ RNx ×(Nt +1) . (10) Each column ui is a vector with Nx components, representing the numerical measurements. For simplicity of description, we consider here real data ui ∈ RNx . Following the Koopman decomposition assumption [17], we consider that a propagator matrix A exists, that maps every column vector onto the next one. The DMD algorithm constructs the best approximation of the propagator matrix A, i.e. n
o
u0 , u1 = Au0 , u2 = Au1 = A2 u0 , . .., uNt = AuNt −1 = ANt u0 .
(11)
The next computational step consists in forming two time–shifted data matrices from the snapshot sequence. A matrix V0 is formed with the first Nt columns and the matrix V1 contains the last Nt columns of V : V0 =
h
u0 u1 ... uNt −1
i
∈ RNx ×Nt , V1 =
h
u1 u2 ... uNt
i
∈ RNx ×Nt .
(12)
For a sufficiently long sequence of the snapshots, we suppose that the last snapshot uNt can be written as a linear combination of previous Nt vectors, such that uNt = c0 u0 + c1 u1 + ... + cNt −1 uNt −1 + R, 4
(13)
in which ci ∈ R, i = 0,..., Nt − 1 and R is the residual vector. We assemble the following relations {u1 , u2 , ...uNt } = A {u0 , u1 , ...uNt −1 } = {u1 , u2 , ..., V0 c} + R,
(14)
T
where c = c0 c1 ... cNt −1 is the unknown column vector. In matrix notation form, Eq. (14) reads
...
0 1 S= .. .
AV0 = V0 S + R,
0 0 .. .
c0 c1 .. .
.. . 0 . . . 1 cNt −1
,
(15)
where S is the companion matrix. Relation (15) is true when the residual R = uNt − V0 c,
(16)
is minimized when c is chosen such that R is orthogonal to span {u0 , ..., uNt −1 }. The goal of DMD algorithm is to solve the eigenvalue problem of the companion matrix S V1 = AV0 = V0 S + R,
(17)
where S approximates the eigenvalues of A when ∥R∥2 → 0, which is equivalent to solve the minimization problem M inimize R = ∥V1 − V0 S∥2 . (18) S
In our previous work [16], we estimate the solution to the minimization problem (18) multiplying V1 by the Moore-Penrose pseudoinverse [36] of V0 : S = (V0 )+ V1 .
(19)
As we previously pointed out in [16], the Moore-Penrose pseudoinverse approach might not be feasible when dealing with high dimensional data. Following Schmid [37] who was the first to introduce the DMD as a numerical tool to compute the Koopman modes, we developed an alternate algorithm based on Singular Value Decomposition (SVD) of snapshot matrix V0 . This approach is helpful especially when the matrix V0 is rank deficient (Nx >Nt ). In the following, we describe this technique. We first identify a singular value decomposition of V0 : V0 = U ΣW H ,
(20)
where U contains the proper orthogonal modes of V0 , Σ is a square diagonal matrix containing the singular values of V0 and W H is the conjugate transpose of W . A direct consequence of solving the minimization problem (18) is that decreasing the residual increases overall convergence and therefore the eigenvalues λj and the eigenvectors ϕj , j = 1, ..., Nt of S will converge toward the eigenvalues and the eigenvectors of the Koopman operator A, respectively. More specifically, every column vector ui , i = 1, ..., Nt can be written as a linear combination of its predecessor: ui = Aui−1 = . . . = Ai−1 u1 ,
i = 1, ..., Nt .
(21)
The eigenvectors of S form a basis for the span of A, therefore, we can write every column vector as a linear combination of the eigenvectors ui =
Nt X
e j ϕj Ai−1 a
⇔
ui =
j=1
Nt X j=1
5
ej λi−1 a j ϕj ,
i = 1, ..., Nt .
(22)
A straightforward interpretation of relations (22) brings the data snapshots at every time step {t1 ,...,tNt } as a linear combination of DMD modes according to V1 =
h
u1 u2 ... uNt
i
= λ11 λ12 .. 1 . .. . . . .. .. e Nt a 1 λ1Nt
=
h
ϕ1 ϕ2
i ... ϕNt
1 1
e1 a e2 a
λ21 λ22 .. . .. . 2 λNt
t −1 . . . λN 1 Nt −1 . . . λ2 , .. .. . . .. .. . . Nt −1 . . . λ Nt
(23)
where the right eigenvectors of S, ϕj ∈ C are dynamic shape (or Koopman) modes, the eigenvalues ej ∈ C are denoted as amplitudes or Koopman of S, λj are called Ritz values [38] and coefficients a eigenfunctions. Each Ritz value λj = e(σj +iωj )∆t is associated with the growth rate σj and the frequency ωj , where log (|λj |) arg (|λj |) σj = , ωj = . (24) ∆t ∆t The superposition of all Koopman modes, weighted by their amplitudes and complex frequencies, approximates the entire data sequence, but there are also modes that have a weak contribution. Our goal is to produce the ROM involving only the most significant modes, having a strong contribution to the data representation, which we are calling leading modes. Thus, the data snapshots at every time step {t1 ,...,tNt } will be represented as a linear combination of the leading DMD modes according to uDM D (ti , x) =
NX DM D
ej ϕj (x) λi−1 a j ,
i ∈ {1, ..., Nt } ,
ti ∈ {t1 , ..., tNt } ,
(25)
j=1
where NDM D represents the number of leading DMD modes involved in reconstruction of data snapshots. One advantage of DMD is that each mode is associated with a pulsation, a growth rate and each mode oscillates at a single frequency, as seen from (25). Representation (25) is suitable when one wants to isolate a mode with a certain frequency, or to identify a maximum or a minimum amplitude and for hydrodynamic stability analysis also. In the seminal article [31], Mezić provides a characterization of Koopman modes in Banach spaces using Generalized Laplace Analysis. For the purpose of model order reduction, in our paper we adopt the following form uDM D (ti , x) =
NX DM D
aj (ti ) ϕj (x),
ti ∈ {t1 , ..., tNt } ,
(26)
j=1
ej λi−1 where ϕj ∈ C are dynamic leading modes and aj (ti ) = a j , i ∈ {1, ..., Nt } , j ∈ {1, ..., NDM D } are modal amplitudes. Here we point out that the NDM D leading modes involved in ROM representation of data (26) are not the first NDM D modes from representation (23). The leading modes represent a subset of DMD modes that will be selected from all computed DMD modes via numerical algorithm presented in the next section.
4
Offline stage: Randomized Dynamic Mode Decomposition
The modes’ selection plays a central role in model reduction and constitutes also the source of many discussions among modal decomposition practitioners [19, 21, 39, 40]. Several procedures for selecting the most influential modes in dynamic mode decomposition can be found in our 6
previous papers [13, 16, 41]. We have introduced in [28] the procedure of randomization of data prior to singular value decomposition (SVD). Thus, we endow the DMD algorithm with a randomized SVD function, aiming to improve the accuracy of the reduced order linear model and to reduce the CPU time. The major advantage of this method is that does not require an additional selection algorithm of the DMD modes. The randomized DMD produces a reduced order subspace of Ritz values, having the same dimension as the rank of randomized SVD function, where the leading modes live. The second advantage consists in reducing the problem dimension to avoid a computationally expensive SVD. The objective of the DMD-based ROM is to represent, as accurately as possible, the high fidelity solution using the dynamics given by the DMD modes. It is then natural to seek the leading DMD modes and their temporal eigenfunctions that minimize the error EDM D = ⟨∥u (x, t) − uDM D (x, t)∥2 ⟩T ,
(27)
where ⟨·⟩T is a time average operator over [t1 , T ] and ∥ · ∥2 is the L2 -norm of RNx . In this paper, ⟨·⟩T corresponds to the arithmetic time-average on Nt equally spaced elements of the interval [t1 , T ]: Nt 1 X f (ti ), ti ∈ {t1 , t2 , ..., tNt = T } . (28) ⟨f (t)⟩T = Nt i=1 Determination of the optimal rank NDM D of the ROM then amounts to finding the solution to the following constrained optimization problem: N
F ind
∈N,N
DM D DM D Subject to N
≥2
uDM D (ti , x) =
NDM PD j=1
aj (ti ) ϕj (x),
ti ∈ {t1 , ..., tNt } ,
(29)
DM D = arg min {EDM D } ,
where EDM D is the error of the low-rank model defined by Eq. (27). Generally, DMD does not produce orthogonal modes, therefore ROMs produced via DMD require a closure model consisting of additional regularization techniques [42–44] especially when applying Galerkin or Petrov–Galerkin projection based techniques. In this paper we propose a variant of randomized dynamic mode decomposition introduced in [28] augmented with deep learning artificial intelligence that confers multiple advantages to the ROM, which will be presented in the following. The algorithm proceeds as follows: Algorithm 1: Randomized Dynamic Mode Decomposition Initial data: V0 ∈ RNx ×Nt , V1 ∈ RNx ×Nt , integer target rank k ≥ 2 and k < Nt . (1:) For k = 2 to Nt − 1. (2:) Produce the randomized singular value decomposition of rank k [U, Σ, W ] = k-RSVD (V0 , k) , where U contains the proper orthogonal modes of V0 and Σ contains the singular values. The RSVD function is described in continuation of this algorithm. (3:) Solve the minimization problem (18). 4. Compute dynamic modes solving the eigenvalue problem SX = XΛ and obtain dynamic modes as Φ = U X. The diagonal entries of Λ represent the eigenvalues λ. 5. Project dynamic modes onto the first snapshot to calculate the vector containing dynamic rank(Λ) modes amplitudes Ampl = (aj )j=1 . 7
6. The DMD model of rank k is given by the product VDM D = Φ · diag (Ampl) · V an,
(30)
where the Vandermonde matrix is −2 . . . λN 1 N −2 . . . λ2 .. .. V an = . . . 1 ... ... ... ... ... −2 1 λ1k λ2k . . . λN k
1 1
λ11 λ12 .. .
λ21 λ22 .. .
7. Solve the optimization problem (29) and obtain the optimal low rank k and associated VDM D . Output: k, VDM D .
The following routine is used to produce the randomized singular value decomposition. Algorithm 2: Randomized Singular Value Decomposition of Rank k (k-RSVD) Initial data: V0 ∈ RNx ×Nt , integer target rank k ≥ 2 and k < Nt . (1:) Generate random test matrix M = rand (Nt , r), r = min (Nt , 2k). (2:) Compute sampling matrix by multiplication of snapshot matrix with random matrix Q = V0 M . (3:) Orthonormalization of sampling matrix via Gram–Schmidt orthonormal method Q ← GramSchmidt (Q). (4:) Projection of snapshot matrix to smaller space V = QH V0 , where H denotes the conjugate transpose. (5:) Produce the economy-size singular value decomposition of low-dimensional snapshot matrix [T, Σ, W ] = SV D (V ). (6:) Compute the right singular vectors U = QT . Output: Procedure returns U ∈ RNx ×k , Σ ∈ Rk×k , W ∈ RNt ×k .
An intuitive understanding of k-RSVD is illustrated in Figure 2. We avoid a computationally expensive algorithm and we reduce the problem dimension by using the randomized singular value decomposition (RSVD) technique.
8
Intuitive Understanding of k-RSVD
V0
1.
M
Q =
2.
QH
k<min(Nx,Nt)
V0
P =
3.
P
WH
T SVD(P)
4.
Q
T
U =
Figure 2: An intuitive understanding of k-RSVD.
5
Online stage: Fast Digital Twin Data Model Identification Using Deep Learning Nonlinear Autoregressive Estimators
The algorithm previously described allows the identification of the leading dynamic modes and their associated temporal coefficients in discrete form. The goal in this section is the identification of the reduced order digital twin data model (DTM) of the form: uROM DT M (t, x) =
NX DM D
bj (t) ϕj (x), a
t ∈ [0, T ] ,
(31)
j=1
bj (t), j = 1, ..., NDM D represent the where ϕj , j = 1, ..., NDM D are the DMD modes and a temporal coefficients of the DTM. Nonlinear AutoRegressive models with eXogenous inputs (NLARX) represent a novel approach in the field of nonlinear system identification [45–47]. Since the emergence of artificial neural networks as numerical tools, NLARX models have been used for various purposes, ranging from simulation [48], to nonlinear predictive control [49] or higher order nonlinear optimization problems [50, 51, 54, 55]. We will investigate in this paper the application of NLARX models to a high-fidelity approximation of temporal coefficients of the DMD-ROM model (31). Let a (t) be the system input represented by the DMD computed amplitudes at discrete time instances t ∈ {t1 , ..., tNt } and b (t) be the output. The formulation of the NLARX model can be described as: a b (t) = f [a b (t − 1) , ..., a b (t − na ) , a (t − nk ) , ..., a (t − nk − nb + 1)] + e (t) , a
(32)
where the na is the integer number of past output terms, nb is the number of past input terms used to predict the current output, nk is the pure input delay, f is a nonlinear function (typically 9
implemented by an artificial neural network) and e (t) represents the modeling error. Each output of NLARX model (32) is a function of regressors that are transformations of past inputs and past outputs. Usually this function has a linear block and a nonlinear block. The model output is the sum of the outputs of the two blocks. The NLARX model training can be cast as a non-linear unconstrained optimization problem: θ (na , nb , nk ) = arg min
Nt 1 X b (ti )∥2 , ∥a (ti ) − a 2Nt i=1
(33)
b (t) is the NLARX output, ∥ · ∥2 is where the training set consists of the measured input a (t), a the L2 norm and θ (na , nb , nk ) represents the parameter vector of the nonlinear function f . The NLARX structure can accommodate the dynamics of the system by feeding previous network outputs back into the input layer. It also enables the user to define how many previous output and input time steps are required for a best representation of the systems dynamics. One of most important points of an application of NLARX network is a proper selection of inputs, input delays and output delays. The task will be to modify the network parameters θ (na , nb , nk ) over the complete trajectory to achieve the minimal value of (33). We have implemented the nonlinear estimator f in form of a cascade forward neural network with 10 hidden layer sizes, see Figure 3. Cascade-Forward Neural Network
Figure 3: The cascade forward neural network with 10 hidden layer sizes, used as nonlinear estimator for the NLARX models. In the next section, we will detail the numerical results.
6
Numerical Results: Computational Efficiency of the Algorithm
In the following, we present numerical results demonstrating the computational performance of the algorithm, considering the nonlinear viscous Burgers equation model (1) generating three shock wave phenomena with increasing complexity. The randomization of input data has been leveraged to accelerate DMD computations. The optimal rank of the reduced DMD model is the unique solution to the optimization problem (29). We have tested several global optimization methods like genetic algorithm combined with sequential quadratic programming (GA-SQP) [52] and simulated annealing (SA) [53], to solve the optimization problem (29), with similar computational efforts. A major advantage that comes from application of randomized DMD algorithm is that this leads to the optimal low rank NDM D and associated DMD subspace where the most influential DMD modes are identified.
10
The correlation coefficient defined below is used as additional metric to validate the quality of the low-rank DMD model: CDM D = D
⟨∥u (x, t) · uDM D (x, t)∥2 ⟩T 2 u(x, t)H · u (x, t)
E D 2 T
uDM D (x, t)H · uDM D (x, t)
E ,
(34)
2 T
where u (t, x) means the numerical data, uDM D (t, x) represent the computed solution by means of the reduced order DMD model, (·) represents the Hermitian inner product, H denotes the conjugate transpose and ⟨·⟩T is the norm defined by Eq. (28). Figures 4–6 present the process of evaluation of DMD model target rank. The error computed as a function of retained number of dynamic modes and the correlation coefficient are presented, respectively, in the cases Re = 102 , Re = 103 and Re = 104 . Table 1 presents the order of DMD subspace obtained in the three test cases, next to the error defined by Eq. (27) and correlation coefficient defined by Eq. (34).
Figure 4: Case of Re = 102 : a) The relative error computed as a function of retained number of dynamic modes, b) The correlation coefficient computed as a function of retained number of dynamic modes. NDM D = 15 leading modes have been selected.
Table 1: Comparison of the numerical results returned by the randomized dynamic mode decomposition algorithm. 300 data snapshots have been processed. Test case
Model rank
Error
Correlation coefficient
Re = 102 Re = 103 Re = 104
NDM D = 15 NDM D = 20 NDM D = 20
EDM D = 3.1984 × 10−7 EDM D = 4.4247 × 10−7 EDM D = 5.4416 × 10−8
CDM D = 1.0000 CDM D = 1.0000 CDM D = 1.0000
The algorithm introduced in this paper confers the best correlation coefficient to the DMD model (see Table 1), thus we have identified a digital twin data model. The DMD leading modes are illustrated in Figures 7-9, next to the representation of the modal growth rates and the associated frequencies of the eigenvectors of the Koopman matrix S for the three test cases, respectively. bj (t), j = 1, ..., NDM D of the reduced order model (31) have been estimated The coefficients a for the entire time window by considering the DMD computed coefficients as inputs of the 11
Figure 5: Case of Re = 103 : a) The relative error computed as a function of retained number of dynamic modes, b) The correlation coefficient computed as a function of retained number of dynamic modes. NDM D = 20 leading modes have been selected.
Figure 6: Case of Re = 104 : a) The relative error computed as a function of retained number of dynamic modes, b) The correlation coefficient computed as a function of retained number of dynamic modes. NDM D = 20 leading modes have been selected. NLARX model (32). The numbers of the input terms, output terms and the value of delay are presented in Tables 2-4, for the three test cases, respectively. The following metrics have been used to perform a qualitative analysis of the digital twin data models (DTMs): D E EDT M = u (x, t) − uROM (x, t) , (35) DT M 2 T
D
CDT M = D
H
u (x, t) · uROM DT M (x, t)
u(x, t) · u (x, t)
E D 2 T
E 2
2 T E , H ROM uDT M (x, t) · uROM DT M (x, t) 2 T
(36)
where EDT M measures the error of the digital twin data model, CDT M is the correlation coefficient of the digital twin data model, u (t, x) means the numerical data, uROM DT M (t, x) represent the 12
Figure 7: Case of Re = 102 : a) The DMD leading modes, b) Growth rates and associated frequencies (σ, ω) of the eigenvectors of the Koopman matrix S.
Figure 8: Case of Re = 103 : a) The DMD leading modes, b) Growth rates and associated frequencies (σ, ω) of the eigenvectors of the Koopman matrix S.
Figure 9: Case of Re = 104 : a) The DMD leading modes, b) Growth rates and associated frequencies (σ, ω) of the eigenvectors of the Koopman matrix S. computed solution by means of the DMD-ROM model, (·) represents the Hermitian inner product,
13
H denotes the conjugate transpose and ⟨·⟩T is the norm defined by Eq. (28). Table 2: Initial data for NLARX estimator of temporal coefficients, case of Re = 102 , NDM D = 15. Index
Outputs, inputs, delay
DTM Error and Correlation coefficient
j = 1, 3 − 12, 14 j=2 j = 13, 15
na = 1, nb = 2, nk = 1 na = 2, nb = 2, nk = 2 na = 1, nb = 1, nk = 1
EDT M = 8.0806 × 10−4 CDT M = 1.0000
Table 3: Initial data for NLARX estimator of temporal coefficients, case of Re = 103 , NDM D = 20. Index
Outputs, inputs, delay
DTM Error and Correlation coefficient
j = 1, 4, 10 − 12 j = 2, 3, 13 j =5−9 j = 14 − 19 j = 20
na = 2, nb = 1, nk = 1 na = 1, nb = 1, nk = 2 na = 1, nb = 1, nk = 5 na = 2, nb = 1, nk = 2 na = 1, nb = 1, nk = 1
EDT M = 1.3080 × 10−4 CDT M = 1.0000
Table 4: Initial data for NLARX estimator of temporal coefficients, case of Re = 104 , NDM D = 20. Index
Outputs, inputs, delay
DTM Error and Correlation coefficient
j = 1, 3, 4, 7 j = 2, 5, 6, 8, 16 − 20 j = 9 − 15
na = 1, nb = 3, nk = 2 na = 1, nb = 2, nk = 2 na = 1, nb = 1, nk = 1
EDT M = 1.4000 × 10−4 CDT M = 1.0000
Solution of the digital twin data models are illustrated in Figures 10-12, for the three test cases, respectively. The very good correlation coefficients and the low value of errors presented in Tables 2-4 confirm the computational efficiency of the digital twin data models. The CPU time required in the offline-online stage is presented in Figure 13, for the three test cases. The required offline CPU time does not exceed two seconds and does not present large variations depending on the case study. The online CPU time fall between 8 and 17 seconds, depending on the index of the temporal coefficient which is estimated along the entire time window. It is obvious that the NLARX estimator requires more time to estimate the temporal behaviour in the case of very high Reynolds number. Figures 14-16 illustrate the validation for the first temporal coefficient, as simulated response of the optimal NLARX estimator, in the case of the three experiments, respectively.
7
Conclusions
The present investigation has focused on the identification of high-fidelity digital twin data models from numerical code outputs by non-intrusive techniques (i.e. not requiring Galerkin projection of the governing equations onto the reduced modes basis). In this paper we define the
14
a.
b.
Figure 10: Solution of the digital twin data model in the case of experiment Re = 102 : a) 3D view; b) Projection view of the shock wave.
a.
b.
Figure 11: Solution of the digital twin data model in the case of experiment Re = 103 : a) 3D view; b) Projection view of the shock wave.
a.
b.
Figure 12: Solution of the digital twin data model in the case of experiment Re = 104 : a) 3D view; b) Projection view of the shock wave. concept of digital twin data model (DTM) as a model of reduced complexity that has the main feature to mirror the original process behavior. We developed an algorithm that utilizes a variant of adaptive randomized dynamic mode
15
Figure 13: The CPU time required in the offline-online stage, for the three test cases.
a.
b.
Figure 14: The validation for the first temporal coefficient, as simulated response of the optimal NLARX estimator, in the case of experiment Re = 102 . decomposition introduced in [28] to obtain a reduced basis in the offline stage, combined with a fast digital twin identification using neural network based nonlinear autoregressive estimators in the online stage. To overcome the inconveniences of developing and implementing a mode selection criterion associated with dynamic mode decomposition, we developed a technique based on randomized dynamic mode decomposition as a fast and accurate option in model order reduction. The rank of the ROM is given as the unique solution of an optimization problem whose constraint consists in the smallest error of DTM. Solving the optimization problem (29) using a hybrid simulated annealing [53] we gain a fast and accurate randomized DMD algorithm, with a low rank for the ROMs. The DTMs have been investigated in the numerical simulation of three shock wave phenomena with increasing complexity, with Reynolds number varying from 102 to 104 . We showed that the significant advantage of DTM is to map the dynamics with high accuracy and reduced costs in CPU time and hardware, even to timescales difficult to explore because of the rapidly changing dynamics over time. The procedure of online estimation of the DTM temporal coefficients by employing neural network based nonlinear autoregressive estimators leads to a fast and accurate identification of the digital twin data models, as seen in Tables 2-4. We investigated the computational efficiency 16
a.
b.
Figure 15: The validation for the first temporal coefficient, as simulated response of the optimal NLARX estimator, in the case of experiment Re = 103 .
a.
b.
Figure 16: The validation for the first temporal coefficient, as simulated response of the optimal NLARX estimator, in the case of experiment Re = 104 . of the proposed algorithm and we provided a qualitative analysis of the DTM in the three experiments investigated.
References [1] Y. Cao, J. Zhu, I. Navon, Z. Luo, A reduced order approach to four-dimensional variational data assimilation using proper orthogonal decomposition, International Journal for Numerical Methods in Fluids 53 (10) (2007) 1571–1583. [2] G. Dimitriu, R. Stefanescu, I. M. Navon, POD-DEIM approach on dimension reduction of a multi-species host-parasitoid system, Ann. Acad. Rom., Sci. Ser. Math. Appl. 7 (1) (2015) 173–188. [3] J. Du, F. Fang, C. C. Pain, I. M. Navon, J. Zhu, D. Ham, POD reduced order unstruc-
17
tured mesh modelling applied to 2d and 3d fluid flow, Computers and Mathematics with Applications 65 (2013) 362–379. [4] R. Stefanescu, I. M. Navon, POD/DEIM nonlinear model order reduction of an ADI implicit shallow water equations model, Journal of Computational Physics 237 (2013) 95–114. [5] Z. Wang, I. Akhtar, J. Borggaard, T. Iliescu, Proper orthogonal decomposition closure models for turbulent flows: A numerical comparison, Computer Methods in Applied Mechanics and Engineering 237-240 (2012) 10–26. [6] N.K. Bellam Muralidhar, N. Rauter, A. Mikhaylenko, R. Lammering, D.A. Lorenz, Parametric Model Order Reduction of Guided Ultrasonic Wave Propagation in Fiber Metal Laminates with Damage, Modelling 2 (2021) 591—608. doi:10.3390/modelling2040031. [7] I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics 41 (1-3) (2005) 309–325. [8] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, D. S. Henningson, Reduced-order models for flow control: balanced models and Koopman modes, in: Seventh IUTAM Symposium on Laminar-Turbulent Transition, IUTAM Bookseries, Vol. 18, 2010, pp. 43–50. [9] S. Bagheri, Koopman-mode decomposition of the cylinder wake, Journal of Fluid Mechanics 726 (2013) 596–623. [10] P. J. Schmid, D. Violato, F. Scarano, Decomposition of time-resolved tomographic PIV, Springer-Verlag, 2012. [11] O. Frederich, D. M. Luchtenburg, Modal analysis of complex turbulent flow, in: The 7th International Symposium on Turbulence and Shear Flow Phenomena (TSFP-7), Ottawa, Canada„ 2011. [12] M. J. Balajewicz, E. H. Dowell, B. R. Noack, Low-dimensional modelling of high-Reynoldsnumber shear flows incorporating constraints from the Navier–Stokes equation, Journal of Fluid Mechanics 729 (2013) 285–308. [13] D. A. Bistrian, I. M. Navon, The method of dynamic mode decomposition in shallow water and a swirling flow problem, International Journal for Numerical Methods in Fluids 83 (2016) 73–89. doi:10.1002/fld.4257. [14] X. Chen, I. M. Navon, F. Fang, A dual-weighted trust-region adaptive POD 4D-VAR applied to a finite-element shallow-water equations model, International Journal for Numerical Methods in Fluids 65 (2011) 250–541. [15] G. Dimitriu, R. Stefanescu, I. Navon, Comparative numerical analysis using reduced-order modeling strategies for nonlinear large-scale systems, Journal of Computational and Applied Mathematics 310 (2017) 32–43. [16] D. A. Bistrian, I. M. Navon, An improved algorithm for the shallow water equations model reduction: Dynamic mode decomposition vs POD, International Journal for Numerical Methods in Fluids 78 (9) (2014) 552–580. doi:10.1002/fld.4029. [17] B. Koopman, Hamiltonian systems and transformations in Hilbert space, Proc. Nat. Acad. Sci. 17 (1931) 315–318. [18] P. J. Schmid, J. Sesterhenn, Dynamic mode decomposition of numerical and experimental data, in: 61st Annual Meeting of the APS Division of Fluid Dynamics, Vol. 53(15), American Physical Society, San Antonio, Texas, 2008. 18
[19] K. K. Chen, J. H. Tu, C. W. Rowley, Variants of dynamic mode decomposition: boundary condition, Koopman and Fourier analyses, Nonlinear Science 22 (2012) 887–915. [20] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, J. N. Kutz, On dynamic mode decomposition: Theory and applications, Journal of Computational Dynamics 1 (2) (2014) 391–421. [21] M. R. Jovanovic, P. J. Schmid, J. W. Nichols, Low-rank and sparse dynamic mode decomposition, Center for Turbulence Research Annual Research Briefs (2012) 139–152. [22] J. N. Kutz, X. Fu, S. L. Brunton, N. B. Erichson, Multi-resolution dynamic mode decomposition for foreground/background separation and object tracking, in: 2015 IEEE International Conference on Computer Vision Workshop (ICCVW), no. INSPEC Accession Number:15790263, IEEE, Santiago, 2015, pp. 921 – 929. [23] J. N. Kutz, X. Fu, S. L. Brunton, Multi-resolution dynamic mode decomposition, SIAM Journal of Applied Dynamical Systems 15 (2016) 713–735. [24] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data–driven approximation of the Koopman operator: extending dynamic mode decomposition, Journal of Nonlinear Science 25 (6) (2015) 1307–1346. [25] B. R. Noack, W. Stankiewicz, M. Morzynski, P. Schmid, Recursive dynamic mode decomposition of a transient cylinder wake, physics.flu-dyn arXiv:1511.06876v1. [26] B. R. Noack, W. Stankiewicz, M. Morzynski, P. Schmid, Recursive dynamic mode decomposition of transient and post-transient wake flows, Journal of Fluid Mechanics 809 (2016) 843–872. [27] J. L. Proctor, S. L. Brunton, J. N. Kutz, Dynamic mode decomposition with control, SIAM Journal of Applied Dynamical Systems 15 (1) (2016) 142–161. [28] D. A. Bistrian, I. M. Navon, Randomized dynamic mode decomposition for nonintrusive reduced order modelling, International Journal for Numerical Methods in Engineering 112 (2017) 3–25. doi:10.1002/nme.5499. [29] D. A. Bistrian, I. M. Navon, Efficiency of randomized dynamic mode decomposition for reduced order modelling, International Journal of Computational Fluid Dynamics. doi: 10.1080/10618562.2018.1511049. [30] S.E. Ahmed, P.H. Dabaghian, O. San, D.A. Bistrian, I.M. Navon, Dynamic mode decomposition with core sketch, Physics of Fluids (2022). doi:10.1063/5.0095163. [31] I. Mezić, On Numerical Approximations of the Koopman Operator, Mathematics 10 (2022), 1180. doi.org/10.3390/math10071180. [32] A. Alla, J. N. Kutz, Nonlinear model order reduction via dynamic mode decomposition, SIAM Journal on Scientific Computing 39 (5) (2016) B778–B796. [33] J. M. Burgers, A mathematical model illustrating the theory of turbulence, Advances in Applied Mechanics 1 (1948) 171–199. [34] C.M. Cuesta, I. Pop, Numerical schemes for a pseudo-parabolic Burgers equation: Discontinuous data and long-time behaviour, Journal of Computational and Applied Mathematics 224 (1) (2009) 269–283. [35] J. N. Kutz, J. L. Proctor, S. L. Brunton, Koopman theory for partial differential equations, arXiv:1607.07076 (2016). 19
[36] G. Golub, C. F. van Loan, Matrix Computations, Third Edition, The Johns Hopkins University Press, 1996. [37] P. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656 (2010) 5–28. [38] A. K. Chopra, Dynamics of Structures, 4th Edition, Prentice-Hall International Series in Civil Engineering and Engineering Mechanics, 2000. [39] B. R. Noack, M. Morzynski, G. Tadmor, Reduced-Order Modelling for Flow Control, Springer, 2011. [40] G. Tissot, L. Cordier, N. Benard, B. R. Noack, Model reduction using dynamic mode decomposition, Comptes Rendus Mecanique 342 (2014) 410–416. [41] A. K. Alekseev, D. A. Bistrian, A. E. Bondarev, I. M. Navon, On linear and nonlinear aspects of dynamic mode decomposition, International Journal for Numerical Methods in Fluids 82 (2015) 348–371. doi:10.1002/fld.4221. [42] L. Cordier, B. A. E. Majd, J. Favier, Calibration of POD reduced-order models using Tikhonov regularization, International Journal for Numerical Methods in Fluids 63 (2010) 269–296. [43] O. San, T. Iliescu, A stabilized proper orthogonal decomposition reduced-order model for large scale quasigeostrophic ocean circulation, Advances in Computational Mathematics 41 (5) (2015) 1289–1319. [44] Y. Wang, I. Navon, X. Wang, Y. Cheng, 2d Burgers equations with large Reynolds number using POD/DEIM and calibration, International Journal for Numerical Methods in Fluids 82 (12) (2016) 909–931. [45] K. Narendra, K. Parthasarathy, Identification and control of dynamic systems using neural networks, IEEE Transactions on Neural Networks 1 (1) (1990) 4–27. [46] A. Juditsky, H. Hjalmarsson, A. Benveniste, B. Deylon, L. Ljung, J. Sloberg, Q. Zhang, Nonlinear black-box models in system identification: Mathematical foundations, Automatica 31 (12) (1995) 1725–1750. [47] L. Ljung, System Identification: Theory for the User, Second Edition, Prentice Hall Information and System Sciences Series, 1999. [48] O. Nelles, Nonlinear System Identification: From Classical Approaches to Neural Networks and Fuzzy Models, Engineering online library. Springer, 2001. [49] G. Liu, V. Kadirkamanathan, S. Billings, Predictive control for non-linear systems using neural networks, International Journal of Control 71 (1998) 1119–1132. [50] J. Moody, The effective number of parameters: An analysis of generalization and regularization in nonlinear learning system, Neural Information Processing Systems (NIPS 1991), Morgan Kaufmann (1991) 847–854. [51] H. Peng, T. Ozaki, V. Haggan-Ozaki, Y. Toyoda, Structured parameter optimization method for the radial basis function-based state-dependent autoregressive model, International Journal of Systems Science 33 (2002) 1087–1098. [52] J. Nocedal, S. J. Wright, Numerical Optimization, Second Edition, Springer, 2006.
20
[53] X. S. Yang, Engineering Optimization: An Introduction with Metaheuristic Applications, Wiley, USA, 2010. [54] Z. Wang, D. Xiao, F. Fang, R. Govindan, C. Pain, Y. Guo, Model identification of reduced order fluid dynamics systems using deep learning, International Journal for Numerical Methods in Fluids 86 (4) (2018) 255–268. [55] T. Tieleman, G. Hinton, Rmsprop: Divide the gradient by a running average of its recent magnitude, COURSERA: Neural Networks for Machine Learning, 2012.
21