2.1. Encoder (inference model)

Consider an observation (measurement) dataset $ \mathcal{D}={\left\{{\mathbf{x}}^{(i)}\right\}}_{i=1}^N $ with $ N $ independent sequences of time series data. Each sequence reflects a multi-DOF time series record, defined as $ {\mathbf{x}}^{(i)}={\left\{{\mathbf{x}}_0,,,{\mathbf{x}}_1,,,\dots, {\mathbf{x}}_t,,,\dots, {\mathbf{x}}_T\right\}}^{(i)} $ , where the observation vector at time instance $ t $ , $ {\mathbf{x}}_t\in {\mathrm{\mathbb{R}}}^m $ , reflects $ m $ monitored DOFs. When the underlying physics equations are known, the observation $ {\mathbf{x}}_t $ at each time instance $ t $ can be assumed to be derived from a corresponding latent (state) variable $ {\mathbf{z}}_t $ , assumed to completely describe the embedded dynamical state. In practice, a common issue is that the latent variables are usually unobserved or only partially observed, via indirect measurements. This limitation is often tackled in prior art via the use of an encoder parameterized by a NN $ {\Psi}_{\mathrm{NN}} $ , which is employed to infer the latent variables from observation data.

In delivering such an estimate, we adopt a temporal version (Girin et al., Reference Girin, Leglaive, Bie, Diard, Hueber and Alameda-Pineda2020) of the VAE (Kingma and Welling, Reference Kingma and Welling2013), that has been implemented in existing literature (Krishnan et al., Reference Krishnan, Shalit and Sontag2017; Yildiz et al., Reference Yildiz, Heinonen and Lähdesmäki2019; Liu et al., Reference Liu, Lai, Bacsa and Chatzi2022). The encoder $ {\Psi}_{\mathrm{NN}} $ can be mathematically described as:

(1a) $$ {\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)={\Psi}_{\mathrm{NN}}\left(\left[\begin{array}{c}{\mathbf{q}}_0\\ {}{\dot{\mathbf{q}}}_0\end{array}\right]|{\mathbf{x}}_{0:{n}_t}\right)=\mathcal{N}\left(\left[\begin{array}{c}{\mu}_{{\mathbf{q}}_0}\\ {}{\mu}_{{\dot{\mathbf{q}}}_0}\end{array}\right],\left[\begin{array}{ll}\operatorname{diag}\left({\sigma}_{{\mathbf{q}}_0}^2\right)& \mathbf{0}\\ {}\mathbf{0}& \operatorname{diag}\left({\sigma}_{{\dot{\mathbf{q}}}_0}^2\right)\end{array}\right]\right), $$

where the first few observations from $ {\mathbf{x}}_0 $ to $ {\mathbf{x}}_{n_t} $ (denoted by $ {\mathbf{x}}_{0:{n}_t} $ ) are used for inferring $ {\mathbf{z}}_0 $ , that is, $ {\mathbf{z}}_0 $ is conditioned on $ {\mathbf{x}}_0 $ to $ {\mathbf{x}}_{n_t} $ ; the latent variables $ {\mathbf{z}}_t\in {\mathrm{\mathbb{R}}}^{2p} $ are assumed to have dimension of $ 2p $ , and the output of the encoder is intentionally split into $ {\mathbf{q}}_0\in {\mathrm{\mathbb{R}}}^p $ and $ {\dot{\mathbf{q}}}_0\in {\mathrm{\mathbb{R}}}^p $ that are corresponding to displacement and velocity states, respectively, that is, $ {\mathbf{z}}_0=\left[\begin{array}{c}{\mathbf{q}}_0\\ {}{\dot{\mathbf{q}}}_0\end{array}\right] $ . It is further assumed that the inferred state variable $ {\mathbf{z}}_0 $ is a stochastic one, which is in this case essential for reflecting uncertainties, and follows a normal distribution, of mean value $ \left[\begin{array}{c}{\mu}_{{\mathbf{q}}_0}\\ {}{\mu}_{{\dot{\mathbf{q}}}_0}\end{array}\right] $ and diagonal covariance matrix $ \left[\begin{array}{ll}\operatorname{diag}\left({\sigma}_{{\mathbf{q}}_0}^2\right)& \mathbf{0}\\ {}\mathbf{0}& \operatorname{diag}\left({\sigma}_{{\dot{\mathbf{q}}}_0}^2\right)\end{array}\right] $ . It should though be noted that it is common to model uncertainty in structural systems, which are subjected to random environmental influences, using a normal distribution. For most of dynamical VAEs frameworks, which are adopted in the context of modeling dynamical systems with uncertainty, the inherent uncertainties are accounted for via the use of normal distributions, as summarized in the work of Girin et al. (Reference Girin, Leglaive, Bie, Diard, Hueber and Alameda-Pineda2020).

In practice, $ {\Psi}_{\mathrm{NN}} $ is comprised of a feed-forward NN (MLP) and a RNN. We assume that the displacement quantity $ {\mathbf{q}}_0 $ only depends on $ {\mathbf{x}}_0 $ , per the assumption adopted in Yildiz et al. (Reference Yildiz, Heinonen and Lähdesmäki2019):

(1b) $$ \left({\mu}_{{\mathbf{q}}_0},{\sigma}_{{\mathbf{q}}_0}^2\right)=\mathrm{MLP}\left({\mathbf{x}}_0\right). $$

The output of this MLP is a stochastic variable of mean $ {\mu}_{{\mathbf{q}}_0} $ and variance $ {\sigma}_{{\mathbf{q}}_0}^2 $ ; the velocity quantity $ {\dot{\mathbf{q}}}_0 $ is inferred from the first leading observations $ {\mathbf{x}}_{0:{n}_t} $ , thus a RNN is implemented to take $ {\mathbf{x}}_0,\hskip0.35em {\mathbf{x}}_1,\hskip0.35em \dots, \hskip0.35em {\mathbf{x}}_{n_t} $ into account:

(1c) $$ \left({\mu}_{{\dot{\mathbf{q}}}_0},{\sigma}_{{\dot{\mathbf{q}}}_0}^2\right)=\mathrm{RNN}\left({\mathbf{x}}_{0:{n}_t}\right), $$

where the output of the RNN is a stochastic variable of mean $ {\mu}_{{\dot{\mathbf{q}}}_0} $ and variance $ {\sigma}_{{\dot{\mathbf{q}}}_0}^2 $ ; $ {n}_t $ need not necessarily reflect a large number, larger $ {n}_t $ might dilute the inference of the velocity quantity; for instance, based on empirical trial, in our implementations $ {n}_t=10 $ . Once the normal distribution defined in equation (1a) is derived, one can sample $ {\mathbf{z}}_0 $ from this distribution, and use it for computing the evolution of the latent dynamics over time. We use $ {\theta}_{\mathrm{enc}} $ to denote all the parameters used in the $ {\Psi}_{\mathrm{NN}} $ , that is, all the hyperparameters involved in the formulation of the MLP and RNN architectures.

2.2. Modeling latent dynamics via Pi-Neural ODEs

There are generally two strategies in terms of how the temporal dependence between states $ \mathbf{z} $ can be modeled. The first strategy is to use a discrete-time model to describe the embedded dynamics, where the first-order Markovian property is assumed. A popular example, in the deep learning context, can be found in the deep Markov models (Krishnan et al., Reference Krishnan, Shalit and Sontag2015, Reference Krishnan, Shalit and Sontag2017; Liu et al., Reference Liu, Lai, Bacsa and Chatzi2022). An alternative lies in adopting continuous models, usually in the form of differential equations, to describe the temporal dependence embedded in the data. The neural ODEs (Chen et al., Reference Chen, Rubanova, Bettencourt and Duvenaud2018) form a recently proposed tool that parameterizes the governing differential equations by feed-forward NNs in a continuous format. A specific merit of a continuous modeling approach is that nonequidistant sequential data can be used for training the model. As the Neural ODEs effectively represent a differential equation construct, the trained model can, in turn, be used as a generative model, meaning as a model which can predict the system response given initial conditions or external excitation.

In the previous work of the authors (Lai et al., Reference Lai, Mylonas, Nagarajaiah and Chatzi2021), we introduced a Pi-Neural ODEs scheme, assuming that a system can be modeled as a superposition of a physics-based modeling term and a learning-based term, where the latter aims to capture the discrepancy between the physics-based model and the actual system. A similar scheme is further discussed in Wagg et al. (Reference Wagg, Worden, Barthorpe and Gardner2020) for application within the context of digital twinning, as the learning-based term allows for adaptation. The scheme is formally described as follows:

(2) $$ \dot{\mathbf{z}}={f}_{\theta_{\mathrm{dyn}}}\left(\mathbf{z}\right)={f}_{\mathrm{phy}}\left(\mathbf{z}\right)+{f}_{\mathrm{NN}}\left(\mathbf{z}\right), $$

where $ {f}_{\mathrm{phy}}\left(\mathbf{z}\right) $ is a physics-based model, which can be built by leveraging the best possible knowledge of the system; $ {f}_{\mathrm{NN}}\left(\mathbf{z}\right) $ is the learning-based model that is materialized as a NN function of $ \mathbf{z} $ . It is noted that the former term $ {f}_{\mathrm{phy}}\left(\mathbf{z}\right) $ is of a fixed and preassigned structure, while the latter term is adjustable during the process of training the model. The parameter vector $ {\theta}_{\mathrm{dyn}} $ , reflects the set of hyperparameters involved in the NN representation $ {f}_{\mathrm{NN}}\left(\mathbf{z}\right) $ .

In this article, we adopt this modeling scheme for use within a reduced order modeling (ROM) setting, to model the latent dynamics of a high-dimensional system. We restrict ourselves in this initiating effort to the modeling of linear or mildly nonlinear systems. The mildly nonlinear system we refer to in this article is that the system can be well approximated by the linearization of the system—the first-order Taylor expansion.

In such a case, an approximation of the dynamics can be derived through the solution of an eigenvalue problem of the structural matrices of the physics-based model (in the case of a nonlinear system, we rely on the linearized part), and is reflected in the following decoupled low-dimensional linearized form:

(3a) $$ \left[\begin{array}{c}\dot{\mathbf{q}}\\ {}\ddot{\mathbf{q}}\end{array}\right]=\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {}-\boldsymbol{\Lambda} & -\boldsymbol{\Gamma} \end{array}\right]\left[\begin{array}{c}\mathbf{q}\\ {}\dot{\mathbf{q}}\end{array}\right], $$

where,

(3b) $$ \boldsymbol{\Lambda} =\left[\begin{array}{cccc}{\omega}_1^2& & & \\ {}& {\omega}_2^2& & \\ {}& & \ddots & \\ {}& & & {\omega}_p^2\\ {}& & & \end{array}\right]\hskip0.24em \boldsymbol{\Gamma} =\left[\begin{array}{cccc}2{\xi}_1{\omega}_1& & & \\ {}& 2{\xi}_2{\omega}_2& & \\ {}& & \ddots & \\ {}& & & 2{\xi}_p{\omega}_p\\ {}& & & \end{array}\right], $$

where $ \boldsymbol{\Lambda} $ and $ \boldsymbol{\Gamma} $ are both diagonal matrices; $ {\omega}_1,\hskip0.35em {\omega}_2,\hskip0.35em \dots, \hskip0.35em {\omega}_p $ are the first $ p $ leading natural frequencies (the first $ p $ maximum frequencies in a descending order) that are retrieved from an eigenanalysis of an a priori available physics-based model; $ {\xi}_1,\hskip0.35em {\xi}_2,\hskip0.35em \dots, \hskip0.35em {\xi}_p $ are the corresponding modal damping ratios; $ \mathbf{I}\in {\mathrm{\mathbb{R}}}^{p\times p} $ denotes the identity matrix.

Our premise is that the physics-based model in equation (3a) does not fully represent the actual system, which implies that the model-derived modal parameters can be different from the parameters that describe the actual operating system as-is, or that additionally, further to the parameters, the structure of the model is lacking. The latter implies that certain mechanisms are not fully understood and are, thus, modeled inaccurately, for instance, mechanisms related to nonlinearities or damping. To account for such sources of error or discrepancies, we add a learning-based term to model the dynamics that are unaccounted for, with equation (3a) now defined as:

(4) $$ \dot{\mathbf{z}}=\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {}-\boldsymbol{\Lambda} & -\boldsymbol{\Gamma} \end{array}\right]\mathbf{z}+\left[\begin{array}{c}\mathbf{0}\\ {}\mathrm{NN}\left(\mathbf{z}\right)\end{array}\right]\hskip0.24em \mathrm{with}\hskip0.24em \mathbf{z}(0)={\mathbf{z}}_0, $$

where $ \mathbf{z}=\left[\begin{array}{c}\mathbf{q}\\ {}\dot{\mathbf{q}}\end{array}\right] $ ; NN represents a feed-forward NN that is a function of $ \mathbf{z} $ . It is noted that the structure presented in equation (4) has the potential of breaking the fully decoupled structure, which is defined by the first term. This is in fact welcomed since the hypothesis of fully decoupled damping matrices, relating to a Rayleigh viscous damping assumption (Craig and Kurdila, Reference Craig and Kurdila2006), is a known source of modeling discrepancies for real-world systems (Satake et al., Reference Satake, Suda, Arakawa, Sasaki and Tamura2003). The learning-based term $ \mathrm{NN}\left(\mathbf{z}\right) $ is thus added to account for possible sources of inconsistency and error. In this physics-informed architecture, during training, the estimated gradients are obtained as the sum of the corresponding gradients derived from the physics-based and learning-based terms. Since the gradients from the physics-based term are fixed, only the gradients of the learning-based term are to be estimated. The combined gradients are restricted in a regime that is closer to the true function’s gradients. Supplementary Appendix further elaborates on the benefit of this physics-informed architecture, which boosts the search for the governing equations close to the actual systems.

The Physics-informed Neural ODE equation (4) governs the evolution of the dynamics. The dynamics of $ \mathbf{z}(t) $ can be solved by numerically integrating $ \mathbf{z}(t)={\int}_{t_0}^t{f}_{\theta_{\mathrm{dyn}}}\left(\mathbf{z}\right) dt $ from $ {t}_0 $ to $ t $ given initial conditions $ {\mathbf{z}}_0 $ , with the estimate of the latent state vector $ \mathbf{z}(t) $ at each time $ t $ offered as:

(5) $$ \mathbf{z}(t)=\mathrm{ODESOLVE}\left({f}_{\theta_{\mathrm{dyn}}},{\mathbf{z}}_0,{t}_0,t\right), $$

where ODESOLVE reflects the chosen numerical integration scheme, with Runge–Kutta methods comprising a typical example of such solvers. The dynamics of the latent state $ \mathbf{z} $ , with realization of $ {\mathbf{z}}_0,\hskip0.35em {\mathbf{z}}_1,\hskip0.35em \dots, \hskip0.35em {\mathbf{z}}_t,\hskip0.35em \dots, \hskip0.35em {\mathbf{z}}_T $ (where $ {\mathbf{z}}_t=\left[\begin{array}{c}{\mathbf{q}}_t\\ {}{\dot{\mathbf{q}}}_t\end{array}\right] $ ), are thus computed at each time step, and can be subsequently fed into the decoder model to reconstruct the full field response, as described in what follows.

2.3. Decoder

In the case of a linear dynamical system, the full-order response $ {\mathbf{x}}_t^{\mathrm{full}}\in {\mathrm{\mathbb{R}}}^g $ comprises a modal representation of $ {\mathbf{x}}_t^{\mathrm{full}}\approx \hskip0.35em {\Phi}_{\hskip0.35em }{\mathbf{q}}_t\hskip0.24em \left({\Phi}_p\in {\mathrm{\mathbb{R}}}^{g\times p};{\mathbf{q}}_t\in {\mathrm{\mathbb{R}}}^p;p\hskip0.35em \le \hskip0.35em g\right) $ , where $ {\Phi}_p $ is the truncated eigenvector matrix, that is, the leading $ p $ columns of full-order eigenvector matrix $ \Phi $ (corresponding to the largest $ p $ eigenvalues).

As illustrated in Figure 1, an estimate of the evolution of the latent state over time $ {\mathbf{z}}_0,\hskip0.35em {\mathbf{z}}_1,\hskip0.35em \dots, \hskip0.35em {\mathbf{z}}_T $ can be obtained by solving the Pi-Neural ODEs via equation (5). It is noted that, within the structural dynamics context, important measurable quantities such as accelerations $ \ddot{\mathbf{q}} $ can further be computed on the basis of the governing equation (4): $ \ddot{\mathbf{q}}=\left[-\boldsymbol{\Lambda} \hskip0.5em -\boldsymbol{\Gamma} \right]\mathbf{z}+\mathrm{NN}\left(\mathbf{z}\right) $ . Thus, beyond the latent states $ \mathbf{q} $ , $ \dot{\mathbf{q}} $ , we can derive further response quantities of interest, such as the acceleration $ \ddot{\mathbf{q}} $ .

Each response quantity can be respectively emitted to the corresponding full-order response vector (involving all structural DOFs) via the decoder $ {\Phi}_p $ ( $ {\mathrm{\mathbb{R}}}^p\to {\mathrm{\mathbb{R}}}^g $ ):

(6a) $$ {\displaystyle \begin{array}{ll}& \mathrm{displacement}:\hskip0.36em {\mathbf{x}}_t^{\mathrm{full}}={\Phi}_p\left({\mathbf{q}}_t\right),\\ {}& \mathrm{velocity}:\hskip0.36em {\dot{\mathbf{x}}}_t^{\mathrm{full}}={\Phi}_p\left({\dot{\mathbf{q}}}_t\right),\\ {}& \mathrm{acceleration}:\hskip0.36em {\ddot{\mathbf{x}}}_t^{\mathrm{full}}={\Phi}_p\left({\ddot{\mathbf{q}}}_t\right),\hskip0.36em \left(t=0,1,\dots, T\right),\\ {}& \end{array}} $$

where $ {\mathbf{x}}_t^{\mathrm{full}} $ , $ {\dot{\mathbf{x}}}_t^{\mathrm{full}} $ , and $ {\ddot{\mathbf{x}}}_t^{\mathrm{full}} $ denote the reconstructed full-order displacement, velocity, and acceleration, respectively. It is noted that further response quantities of interest, such as potentially strains, can be inferred due to availability of a FEM model.

We can only measure a limited of DOFs, $ {\mathbf{x}}_t\in {\mathrm{\mathbb{R}}}^m $ (we use $ {\mathbf{x}}_t $ to denote measured quantities while $ {\hat{\mathbf{x}}}_t $ denoting the corresponding estimated quantities), via use of appropriate sensors, which form a subset of the full response vector:

(6b) $$ {\hat{\mathbf{x}}}_t=\mathbf{E}\left[\begin{array}{c}{\mathbf{x}}_t^{\mathrm{full}}\\ {}{\dot{\mathbf{x}}}_t^{\mathrm{full}}\\ {}{\ddot{\mathbf{x}}}_t^{\mathrm{full}}\end{array}\right], $$

where $ \mathbf{E}\in {\mathrm{\mathbb{R}}}^{m\times 3g} $ is a selection matrix (each row is a one-hot row vector), selecting corresponding monitored quantities; $ {\hat{\mathbf{x}}}_t $ can represent an extended set of the estimated observations, which can correspond to displacement, velocity, acceleration, or further computable response quantities (such as strains). Since we only consider mild nonlinearity, we rely on the observability of the linearized part of the system, where classical observability theory (Kalman, Reference Kalman1960) can be applied to analyze the observability—estimating the full state vector from limited measurements.

The architecture of the proposed framework essentially comprises a sequential version of the VAE, exploiting the presence of an underlying low-dimensional latent representation in the observed dynamics. In the original VAE, the decoder is parameterized by a NN without regularization, which flexibly fits the training data, without necessarily embodying a physical connotation. From an engineering perspective, however, it would be beneficial if the decoder is bestowed with a direct linkage to physical DOFs. One way to achieve this is to seed the modal shape information, computed from physics-based models, which carries within it the spatial information of how each element/node in $ \mathbf{x} $ is interconnected. Therefore, we forcibly implement eigenmodes $ {\Phi}_p $ as the decoder for emitting the latent variables to the observation space. $ {\Phi}_p=\left[{\phi}_1,{\phi}_2,.\dots, {\phi}_p\right] $ , where each column represents a single eigenmode, can be derived from the structural matrices of the physics-based full order model, for example, a FE model. It is noted that $ {\Phi}_p $ is assumed to be time-invariant, thus reflecting an invariant encoding of the spatial relationship between structural DOFs. However, the residual term $ \mathrm{NN}\left(\mathbf{z}\right) $ in the Pi-Neural ODE in equation (4) adaptively accounts for discrepancies that stem from mild nonlinearities, which would also violate the assumption of invariance. We remind that, in Section 2.2, a decoupled structure is adopted as a prior model to encourage the model to mimic the process of a modal decomposition–reconstruction.

It is worth mentioning that, in this framework, the encoder process can be viewed as the transformation from full-order physical coordinates to modal coordinates $ {\Psi}_{\mathrm{NN}}:\mathbf{x}\to \mathbf{z} $ . In real scenarios that involve weakly nonlinear systems, this can be thought of as a “modal-like” coordinate as the learning term $ \mathrm{NN}\left(\mathbf{z}\right) $ can violate the decoupled structure, while the decoder is viewed as the operator which enables the transformation from the modal coordinates’ space to the measured physical coordinates ( $ {\Phi}_p:\mathbf{z}\to \mathbf{x} $ ).

2.4. Loss function

For the purpose of training the suggested Neural Modal ODE models, which capitalize on the availability of physics information and data, we calculate the measurement prediction error. The model delivers an estimate $ {\hat{\mathbf{x}}}_{0:T} $ of the measured response quantities $ {\mathbf{x}}_{0:T} $ , which in turn allows to minimize the error between the predicted and actual observations, to train the model. The training of the encoder, decoder, and latent dynamic models are performed simultaneously, and the loss function of the framework is given as:

(7) $$ \mathcal{L}(\boldsymbol{\theta}; \mathbf{x})=\mathcal{L}\{\mathrm{D}\mathrm{E}\mathrm{C}\mathrm{O}\mathrm{D}\mathrm{E}\mathrm{R}[\mathrm{O}\mathrm{D}\mathrm{E}\mathrm{S}\mathrm{O}\mathrm{L}\mathrm{V}\mathrm{E}({f}_{{\boldsymbol{\theta}}_{\mathrm{dyn}}},{\Psi}_{\mathrm{NN}}({\mathbf{x}}_{0:{n}_t}),{t}_0,T)]\}, $$

where $ \theta ={\theta}_{\mathrm{enc}}\cup {\theta}_{\mathrm{dyn}} $ are all the parameters involved in the deep learning model; $ {\mathbf{x}}_{0:{n}_t} $ is the first $ {\mathbf{x}}_0 $ to $ {\mathbf{x}}_{n_t} $ data fed into the encoder $ {\Psi}_{\mathrm{NN}} $ ; $ {\mathbf{x}}_{0:T} $ is the whole sequence of the data set used for the decoder; the notation DECODER denotes the process given in equations (6a) and (6b).

In the VAE formulation (Kingma and Welling, Reference Kingma and Welling2013), the loss function $ \mathrm{\mathcal{L}} $ is used to maximize a variational lower bound of the data log-likelihood $ \log p\left(\mathbf{x}\right) $ ; here x is short for $ {\mathbf{x}}_{0:T} $ . Using the variational principle with the inference model $ {\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right) $ , which is only used to infer the initial condition of $ {\mathbf{z}}_0 $ , the evidence lower bound (ELBO) of the data log-likelihood, which is the loss function, is given as follows:

(8) $$ \mathrm{\mathcal{L}}\left(\theta; \mathbf{x}\right)=\sum \limits_{t=0}^T\left\{{\unicode{x1D53C}}_{\Psi_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)}\left[\log p\left({\mathbf{x}}_t|{\mathbf{z}}_t\right)\right]-{\unicode{x1D53C}}_{\Psi_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)}\left[\mathrm{KL}\left({\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)\Big\Vert p\left({\mathbf{z}}_0\right)\right)\right]\right\}, $$

where KL stands for the Kullback–Leibler divergence; a statistical measure that evaluates the closeness of two probability distributions $ {p}_1 $ and $ {p}_2 $ , defined as $ \mathrm{KL}\left({p}_1\left(\mathbf{z}\right)\Big\Vert {p}_2\left(\mathbf{z}\right)\right):= \int {p}_1\left(\mathbf{z}\right)\log \frac{p_1\left(\mathbf{z}\right)}{p_2\left(\mathbf{z}\right)}d\mathbf{z} $ . In the loss function, the first term $ {\sum}_{t=0}^T{\unicode{x1D53C}}_{\Psi_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)}\left[\log p\left({\mathbf{x}}_t|{\mathbf{z}}_t\right)\right] $ evaluates the reconstruction accuracy: $ {\mathbf{z}}_0 $ is sampled from the distribution given in equation (1a), and with this given initial condition, one can compute the predicted $ {\hat{\mathbf{x}}}_t\sim \mathcal{N}\left({\hat{\mu}}_t,{\hat{\varSigma}}_t\right)\hskip0.24em \left(t=0,1,\dots, T\right) $ via the latent dynamics model in equation (5) followed by the decoder. Thus, this term can be computed as $ {\sum}_{t=0}^T\log p\left({\mathbf{x}}_t\right) $ given $ {\mathbf{z}}_0\sim {\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right) $ , and $ \log p\left({\mathbf{x}}_t\right) $ has an analytical form when $ p\left({\mathbf{x}}_t\right) $ follows a normal distribution:

(9) $$ \log p\left({\mathbf{x}}_t\right)=-\frac{1}{2}\left[\log |{\hat{\varSigma}}_t|+{\left({\mathbf{x}}_t-{\hat{\mu}}_t\right)}^T{\hat{\varSigma}}_t^{-1}\left({\mathbf{x}}_t-{\hat{\mu}}_t\right)+{d}_{\mathbf{x}}\log \left(2\pi \right)\right], $$

which is the log-likelihood, and the training of the model is expected to maximize this likelihood given the actual observation data $ {\mathbf{x}}_t $ ; $ {d}_{\mathbf{x}} $ is the dimension of $ {\mathbf{x}}_t $ .

The second term $ -{\sum}_{t=0}^T{\unicode{x1D53C}}_{\Psi_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)}\left[\mathrm{KL}\left({\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)\Big\Vert p\left({\mathbf{z}}_0\right)\right)\right] $ evaluates the closeness of the inferred initial condition with a prior distribution $ p\left({\mathbf{z}}_0\right) $ . In practice, $ p\left({\mathbf{z}}_0\right) $ can be assumed as a normal distribution $ \mathcal{N}\left(\mathbf{0},\mathbf{I}\right) $ if no further prior knowledge is given. The KL terms acts as a penalty term when the inferred initial value is distant from the prior distribution. This term can be alternatively computed as $ -{\sum}_{t=0}^T\mathrm{KL}\left({\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)\Big\Vert p\left({\mathbf{z}}_0\right)\right) $ given $ {\mathbf{z}}_0\sim {\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right) $ . $ \mathrm{KL}\left({p}_1\left(\mathbf{z}\right)\Big\Vert {p}_2\left(\mathbf{z}\right)\right) $ is described by an analytical formula when both $ {p}_1\left(\mathbf{z}\right) $ and $ {p}_2\left(\mathbf{z}\right) $ are normal distributions and $ {p}_2\left(\mathbf{z}\right)\sim \mathcal{N}\left(\mathbf{0},\mathbf{I}\right) $ :

(10) $$ \mathrm{KL}\left({\Psi}_{\mathrm{NN}}\left({\mathbf{z}}_0|{\mathbf{x}}_{0:{n}_t}\right)\Big\Vert p\left({\mathbf{z}}_0\right)\right)=-\log \mid \operatorname{diag}\left({\sigma}_{{\mathbf{z}}_0}\right)\mid +\frac{{\left\Vert {\sigma}_{{\mathbf{z}}_0}\right\Vert}^2+{\left\Vert {\mu}_{{\mathbf{z}}_0}\right\Vert}^2}{2}-\frac{d_{\mathbf{z}}}{2}, $$

in which, $ {\sigma}_{{\mathbf{z}}_{\mathbf{0}}}=\left[\begin{array}{c}{\sigma}_{{\mathbf{q}}_{\mathbf{0}}}\\ {}{\sigma}_{{\dot{\mathbf{q}}}_{\mathbf{0}}}\end{array}\right] $ ; $ {\mu}_{{\mathbf{z}}_{\mathbf{0}}}=\left[\begin{array}{c}{\mu}_{{\mathbf{q}}_{\mathbf{0}}}\\ {}{\mu}_{{\dot{\mathbf{q}}}_{\mathbf{0}}}\end{array}\right] $ ; $ \mid \cdot \mid $ is the determinant of a matrix; $ \left\Vert \cdot \right\Vert $ is the modulus of a vector; $ {d}_{\mathbf{z}} $ is the dimension of $ \mathbf{z} $ .

2.5. Prediction of learned dynamics

The completion of the training process results in the definition of the hyperparameter sets $ {\boldsymbol{\theta}}_{\mathrm{enc}} $ and $ {\boldsymbol{\theta}}_{\mathrm{dyn}} $ . This delivers an encoder $ {\Psi}_{\mathrm{NN}} $ together with a learned dynamic model $ \dot{\mathbf{z}}={f}_{\theta_{\mathrm{dyn}}}\left(\mathbf{z}\right) $ , which retains the structure of differential equations. Equations (5) and (6a) can be used for predicting the dynamics given an initial state $ {\mathbf{z}}_0 $ . $ {\mathbf{z}}_0 $ can be either be inferred from the observation dataset via the learned encoder $ {\Psi}_{\mathrm{NN}} $ , or—when using the derived model as a generative model—the modeler can assign other specific values for the initial condition $ {\mathbf{z}}_0 $ .

For those readers that are interested in reusing the developed algorithms, a demonstrative implementation in Python, reproducing all steps from Sections 2.1 to 2.5, will be made available at: https://github.com/zlaidyn/Neural-Modal-ODE-Demo, including both linear and nonlinear cases of a demonstrative example introduced in the next section.

4.1. Experimental setup and data description

As shown in Figure 4a, this model bridge consists of one 6-m continuous beam, 2 towers, and 16 cables. The beam and towers are made of aluminum alloy, and additional metal weights are attached onto the beam and towers, ensuring that the scaled model’s dynamic properties closely approximate those of the real cable-stayed bridge.

Figure 4. Scale model cable-stayed bridge: (a) in situ photo; (b) diagram of the finite element model (unit: mm). The eight deployed accelerometers are labeled as A1, A2,…, A8, with arrows indicating the sensing directions; the coordinate system is defined by the direction of “X–Y” in the diagram.

Cable-stayed bridges are known for exhibiting geometric nonlinearities. Generally speaking, the nonlinear effects in cable-stayed bridges include: (a) cable sag effects: cables sag because of their self-weight, resulting in variation of their axial stiffness; (b) P-delta effects: the horizontal components of cable forces bend the vertically compressed bridge pylon, introducing additional bending moments. The main girder of a cable-stayed bridge also suffers from the P-delta effects, where the bending girder is compressed by the horizontal components of cable forces; (c) large displacement effects: the displacement of the girder can be large, as the main girder of a cable-stayed bridge is mainly supported by flexible cables; thus, the small deformation assumption and linear beam theory do not apply in this scenario.

In this study, the model cable-stayed bridge exhibits nonlinearity in terms of the P-delta and large displacement effects, but these two effects are mild owing to the relatively small dimension of the scaled model. Further, the cable sag effects are negligible, as the steel cables are light. As a result, this scaled cable-stayed bridge model manifests mild nonlinearity and can be well approximated by the proposed scheme.

To measure the dynamic response of the bridge model, as highlighted in Figure 4, eight MEMS (Micro-Electro-Mechanical System) accelerometers—labeled as A1–A8—are deployed on the structure, and a wired connection is used to collect the acceleration data to a digital data acquisition system. Acceleration measurements are collected at a sampling rate of 100 Hz, while the collected raw data is low-pass filtered at 30 Hz, as the dominant power in the spectrum of the raw signal lies below 30 Hz.

A “pull-and-release” action was used to excite the bridge model. A 1 kg iron weight was hung on node 19 with a wire. When the bridge model and the weight were both stationary, the wire was abruptly cut, inducing a damped free vibration of the bridge model, and those bridge responses were recorded with the accelerometers. It is worth mentioning that the weight was hanged at the exact lateral center of the beam, so the out-of-plane vibration such as torsion was supposed to be negligible.

Five repeated tests were performed. Four of these tests were used for training, with the remaining test serving as a testing dataset.

4.2. Finite element modeling

A two-dimensional (2-D) finite element model (FEM) of the scaled bridge has been developed, in a MATLAB environment (MATLAB, 2019), which serves as the physics-based model to be adopted within the proposed deep learning framework. We consider a 2-D model as the expected motion and the deployed sensors lie within a plane. The dimension, boundary conditions, node number, coordinate system, and sensor position of the FEM are displayed in Figure 4b. Each node corresponds to three DOFs: horizontal ( $ x $ ), vertical ( $ y $ ), and rotation. The notation “010” in Figure 4b signifies that vertical movement is restricted, while the horizontal and rotational movement are free (nodes 1, 4, 7, 23, 26, 29 are of this case); “111” signifies that all the three possible DOFs are restricted (nodes 30 and 43 are of this case). The beam and towers are simulated using the Euler–Bernoulli beam element, and the cables are modeled with the tension-only truss element. The total number of DOFs of the FEM model is 153, after applying the boundary conditions.

Eigenanalysis is performed on the FEM model of this bridge, and the first four mode shapes ( $ {\phi}_1 $ to $ {\phi}_4 $ ) and corresponding frequencies ( $ \frac{\omega_1}{2\pi } $ to $ \frac{\omega_4}{2\pi } $ ) are shown in Figure 5, which are a horizontal drifting mode (1.6387 Hz), followed by three vertical bending modes (3.4529, 6.3667, and 11.2516 Hz).

Figure 5. The first four mode shapes (denoted by red lines) derived from the eigenvalue analysis of the FEM model.

4.3. Model implementation

In this example, for modeling the latent dynamics via equation (4), we adopt the first 10 modes to construct the latent dynamics, that is, $ p=10 $ ; $ \boldsymbol{\Lambda} =\operatorname{diag}\left({\omega}_1^2,{\omega}_2^2,\dots, {\omega}_{10}^2\right) $ and $ \boldsymbol{\Gamma} =\operatorname{diag}\left(2{\xi}_1{\omega}_1,2{\xi}_2{\omega}_2,\dots, 2{\xi}_{10}{\omega}_{10}\right) $ . The decoder $ {\Phi}_p=\left[{\phi}_1,{\phi}_2,\dots, {\phi}_{10}\right]\in {\mathrm{\mathbb{R}}}^{153\times 10} $ , mapping the lower dimensional latent variables back to the full order of 153; $ {\phi}_1 $ to $ {\phi}_{10} $ are the first 10 mode shapes.

To train the model, the channels A1, and A3–A8 are used, while it is noted that the channel A2 is left out (considered as “unmeasured”) to be used for evaluating the performance of reconstruction, that is, the model uses the sensor data at a few DOFs to reconstruct a full-order response.

The data set includes multiple repeated free-vibration cases of the bridge, introduced by cutting a string that hangs a 1 kg mass on node 19. The whole data set is divided into batches for training the model, and the number of time steps for each batch is equally 500. Thus, for each batch, the initial conditions are different, which is beneficial for training the encoder of the model. In addition, we normalize the measured acceleration across from A1 to A8, so that the maximum amplitude is 1.0, which is unitless. The details of the involved NNs are listed in Table 3.

Table 3. Implementation details for the experimental study.