2.1. Problem setting

Suppose that a nominal physics-based model and a dataset $ \mathbf{\mathcal{D}}={\left\{\left({\boldsymbol{x}}_{\mathrm{i}},{\boldsymbol{y}}_{\mathrm{i}},{\boldsymbol{c}}_{\mathrm{i}}\right)\right\}}_{\mathrm{i}=1}^N $ , composed of $ N $ triplets of response measurements $ {\boldsymbol{x}}_{\mathrm{i}} $ , domain variables $ {\boldsymbol{c}}_{\mathrm{i}} $ and class variables $ {\boldsymbol{y}}_{\mathrm{i}} $ are available for a given system under investigation. In structural and mechanical engineering applications, the response measurements will often be displacements, strains or accelerations, measured under operating conditions, that describe the static or dynamic performance of the system. The domain variables $ {\boldsymbol{c}}_i $ can be measurements of environmental and operational parameters, such as the location, temperature, humidity or other properties of a structure or sensor. The class variables $ {\boldsymbol{y}}_i $ describe properties of the system that are cumbersome to obtain, such as assessments of the health condition of one or more structural components performed by experts or extracted from inspection reports. It is assumed that $ \boldsymbol{y} $ is a quantity of interest to be predicted for new incoming observations of $ \boldsymbol{x} $ and $ \boldsymbol{c} $ . Our goal is to simultaneously perform reconstruction of $ \boldsymbol{x} $ and regression or classification on $ \boldsymbol{y} $ , and furthermore to utilize the observed domain and class variables to account for the impact of the domain and class influences on the measured response.

To highlight the intended application setting, three examples are presented in Figure 1 consisting of a beam, an oscillator, and population of bridges. In each example the available physics-based model fails to account for domain and class influences that are present in the measured response signals. In all three examples, the domain and class observables provide valuable information about the system that is necessary for accurately inferring a distribution over parameters of the physics-based model and reconstructing the response of the physical system. Furthermore, the class observables are significantly harder to measure and will not be available for future experiments. The oscillator and bridge examples also include unknown confounding influences, for which neither observations nor a physical description are available. The three examples are summarized below:

Figure 1. Illustrative examples of the problem setting: a) Beam, b) Oscillator, and c) One member of a population of bridges. The objective is to learn components of the measured response (bottom row) that are not explicitly included in the nominal physics-based model (top row) using observations of related quantities.

(a) Beam: In this example, the available physics-based model of a beam assumes simply supported boundary conditions and a point load acting on an unknown position. Noisy measurements of the displacement field are obtained from a set of sensors, equally spaced along the length of the beam. The domain influence is introduced as a dependence of the rotational stiffness of the right support on the temperature, and the class influence is taken as damage causing a reduction in the vertical stiffness of the right support that varies between experiments.

(b) Oscillator: Noisy displacement time-series are obtained from multiple experiments, where a mass-spring-dashpot system is deflected from the equilibrium position and released to perform a harmonic oscillation. The impact of damping on the motion is neglected in the physics-based model and must be inferred from additional observations of properties of the medium. The spring stiffness is taken to depend on the ambient temperature, and there is unknown variability in the initial displaced position. The damping, temperature and variability in initial position are taken as class, domain and unknown confounding influences, respectively.

(c) Population of bridges: A vehicle is used to excite the response of a large set of bridges belonging to a homogeneous population, with uncertain vertical stiffness of the supports and varying position of the central pier. Each bridge is monitored by a point strain gauge that yields an influence line for the moving load. The strain gauge measurements are supplemented by qualitative assessments of the condition of the deck, obtained during inspections performed by experts and considered as class variables. The variability in the vehicle velocity and the existence of deterioration in the deck have an influence on the measured strain but are deemed too complicated to model, while the position of the pier is a known domain parameter and can be included in the physics-based model. A variability in the vehicle load is considered as an unknown confounding influence.

2.2. Epistemic uncertainty in Bayesian model updating

Uncertainty in the modeling of physical systems can generally be classified as either aleatoric or epistemic. Aleatoric uncertainty is the component of the uncertainty due to the inherent randomness of a physical process that can not be reduced, while epistemic uncertainty stems from lack of knowledge regarding the physical process (Kiureghian and Ditlevsen, Reference Kiureghian and Ditlevsen2009). Epistemic uncertainty is always present to some degree in practical applications, either due to lack of knowledge or due to simplifications and approximations used to make the evaluation of the physics-based model computationally tractable. The reader is referred to Kamariotis et al., Reference Kamariotis, Vlachas, Ntertimanis, Koune, Cicirello and Chatzi2024 for an extended overview on classification and treatment of uncertainties in SHM applications.

Suppose that an analytical or numerical physics-based model of a structure, defined as a function $ f\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ is available. The response measurements $ \boldsymbol{x} $ can then be expressed as $ \boldsymbol{x}=f\left({\boldsymbol{z}}_{\mathrm{x}}\right)+\epsilon $ , where $ \epsilon $ is a realization of a random variable quantifying the discrepancy between $ f\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ and $ \boldsymbol{x} $ due to the combined influence of aleatoric and epistemic uncertainties. In the Bayesian model updating framework, the measured response of a physical system is used to update the prior knowledge, expressed in terms of a prior distribution $ p\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ over physically meaningful latent variables $ {\boldsymbol{z}}_{\mathrm{x}} $ . This is achieved by approximating the data generating process (i.e. the real-world process that generated the observations) as a combination of a deterministic physics-based model and a probabilistic model (Kennedy and O’Hagan, Reference Kennedy and O’Hagan2001), where the latter accounts for the combined influence of epistemic and aleatoric uncertainty. In this work, it is assumed that the epistemic uncertainty stems from the confounding influence of the domain and class of a structure, and our inability (e.g. due to cost or time constraints) to account for these influences in the form and parameters of the physics-based model.

2.3. The variational autoencoder

In practical applications, the available domain knowledge is often not sufficient to guarantee that the coupled probabilistic-physical model is an accurate description of the data generating process, limiting the applicability of physics-based modeling. To remedy the lack of domain knowledge, data-driven models based on machine learning techniques have emerged as an alternative to physics-based models, where the unknown physical process is learned from measurements using flexible parametrized approximations. VAE (Kingma and Welling, Reference Kingma and Welling2019, Reference Kingma and Welling2022) are a popular data-driven approach for learning a joint distribution of data and the latent variables that are assumed to have generated the data using amortized variational inference (VI) (Blei et al., Reference Blei, Kucukelbir and McAuliffe2017). In VAE, the per-datapoint posterior distribution is approximated using a parametrized family of distributions, where the optimal parameters are obtained by minimizing the Kullback–Leibler divergence (KLD) between the true and approximate posteriors. The VAE is composed of an encoder network $ {q}_{\phi}\left(\boldsymbol{z}|\boldsymbol{x}\right) $ and a decoder network $ {p}_{\theta}\left(\boldsymbol{x}|\boldsymbol{z}\right) $ , parametrized by $ \phi $ and $ \theta $ _, respectively, where $ \boldsymbol{z} $ denotes latent variables that can not be observed directly and must be inferred from measurements. The encoder is typically implemented as a feed-forward NN that maps the inputs $ \boldsymbol{x} $ to a conditional density over latent variables $ \boldsymbol{z} $ . The decoder network $ {p}_{\theta}\left(\boldsymbol{x}|\boldsymbol{z}\right) $ works in the opposite direction by approximating the density of $ \boldsymbol{x} $ conditioned on $ \boldsymbol{z} $ . The training process for VAE consists of simultaneously optimizing the parameters of the decoder that reconstructs the observations given samples of the latent variables, and the encoder that maps inputs to a posterior distribution over these latent variables. Optimization is performed by maximizing a lower bound on the marginal likelihood of the data known as the Evidence Lower BOund (ELBO), denoted as $ {\mathrm{\mathcal{L}}}_{\mathrm{VAE}} $ in Equation (2.1). Sampling $ z\sim {q}_{\phi}\left(\boldsymbol{z}|\boldsymbol{x}\right) $ and evaluating the decoder yields samples from the learned distribution of the data, which in the context of civil and mechanical structural systems can be used for downstream tasks such as remaining useful life assessment.

(2.1) $$ {\displaystyle \begin{array}{c}{\mathrm{\mathcal{L}}}_{\mathrm{VAE}}\left(\boldsymbol{\theta}, \boldsymbol{\phi}; \boldsymbol{x}\right)={\unicode{x1D53C}}_{q_{\boldsymbol{\phi}}\left(\boldsymbol{z}|\boldsymbol{x}\right)}\left[\log {p}_{\boldsymbol{\theta}}\left(\boldsymbol{x}|\boldsymbol{z}\right)\right]-{D}_{\mathrm{KL}}\left({q}_{\boldsymbol{\phi}}\left(\boldsymbol{z}|\boldsymbol{x}\right)\Big\Vert {p}_{\boldsymbol{\theta}}\left(\boldsymbol{z}\right)\right)\\ {}=\log {p}_{\boldsymbol{\theta}}\left(\boldsymbol{x}\right)-{D}_{\mathrm{KL}}\left({q}_{\boldsymbol{\phi}}\left(\boldsymbol{z}|\boldsymbol{x}\right)\Big\Vert {p}_{\boldsymbol{\theta}}\left(\boldsymbol{z}|\boldsymbol{x}\right)\right)\\ {}\le \log {p}_{\boldsymbol{\theta}}\left(\boldsymbol{x}\right)\end{array}} $$

While data-driven approaches might excel in accurately predicting the response of a physical system for a given set of input parameters when sufficient training data is available, the resulting models are typically black boxes that lack interpretability and yield no useful insights about the underlying physical process that generated the measurements. In cases where both the domain knowledge and the available data are limited, purely physics-based or data-driven approaches become infeasible, necessitating a compromise between the two extremes. Physics-enhanced machine learning (PEML) encompasses a wide range of approaches that combine machine learning with domain knowledge (Cicirello, Reference Cicirello2024; Cross et al., Reference Cross, Gibson, Jones, Pitchforth, Zhang, Rogers, Cury, Ribeiro, Ubertini and Todd2022; Haywood-AIexander et al., Reference Haywood-AIexander, Liu, Bacsa, Lai and Chatzi2024; von Rueden et al., Reference von Rueden, Mayer, Beckh, Georgiev, Giesselbach, Heese, Kirsch, Pfrommer, Pick, Ramamurthy, Walczak, Garcke, Bauckhage and Schuecker2023). In this paradigm, the available domain knowledge can be supplemented with data, resulting in more accurate and interpretable models than would be possible with either domain knowledge or data alone. PEML approaches have the potential to reduce the required amount of data, improve accuracy and generalization performance and ensure that model predictions are consistent with the known physics. Importantly, incorporating the known physics can yield interpretable representations of physically meaningful quantities, and models that are robust and explainable.

2.4. Challenges in combining physics-based and data-driven components in VAE

A straightforward approach to account for epistemic uncertainty in a data-driven manner would be to approximate the measured response $ \boldsymbol{x} $ as the sum of the physics-based model $ f\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ and a trainable NN-based function $ {g}_{\boldsymbol{\theta}}\left(\cdot \right) $ , where the latter corrects the discrepancies between the physics-based model predictions and measurements. It is assumed that the gradients of the physics-based model with respect to the inputs can be evaluated efficiently to obtain a computationally tractable optimization problem. This type of hybrid model is referred to as a residual model. Parametrizing the data driven component of the residual model as $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ is not feasible when it is required that $ {\boldsymbol{z}}_{\mathrm{x}} $ is interpretable: The resulting hybrid generative model has a posterior distribution $ {p}_{\theta}\left({\boldsymbol{z}}_{\mathrm{x}}|\boldsymbol{x}\right) $ , where the latent variables $ {\boldsymbol{z}}_{\mathrm{x}} $ are the input to a coupled physics-based and data-driven model, and thus no longer physically meaningful. Instead, the latent space can be partitioned as $ \left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)\in \boldsymbol{z} $ and the data driven component parametrized as $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ , where $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ are physically meaningless latent variables intended to capture variability in the measured response due to the influence of the domain and the class. Assuming that the remaining aleatory uncertainties (e.g., caused by measurement noise) are independent of the signal being measured and can be sufficiently modeled as independent and identically distributed (i.i.d.) samples of Gaussian white noise with standard deviation $ {\sigma}_{\mathrm{x}} $ , the response measurements can be expressed as:

(2.2) $$ \boldsymbol{x}=f\left({\boldsymbol{z}}_{\mathrm{x}}\right)+{g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)+{\boldsymbol{\epsilon}}_{\mathrm{x}}, $$

where $ {\boldsymbol{\epsilon}}_{\mathrm{x}}\sim \mathbf{\mathcal{N}}\left(\mathbf{0},{\sigma}_{\mathrm{x}}^2\boldsymbol{I}\right) $ and $ \boldsymbol{I} $ is the identity matrix. Substituting $ {\hat{\boldsymbol{x}}}_{\mathrm{p}}=f\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ and $ {\hat{\boldsymbol{x}}}_{\mathrm{d}}={g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ for clarity, the resulting generative model is defined as:

(2.3) $$ {p}_{\theta}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right):= \mathbf{\mathcal{N}}\left({\hat{\boldsymbol{x}}}_{\mathrm{p}}+{\hat{\boldsymbol{x}}}_{\mathrm{d}},{\sigma}_{\mathrm{x}}^2\boldsymbol{I}\right) $$

The observables $ \boldsymbol{c} $ and $ \boldsymbol{y} $ can be used to ensure that the latent variables $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ encode information about the domain and class of the structure by simultaneously training auxiliary tasks $ {r}_{\mathrm{c}}\left(\boldsymbol{c}|{\boldsymbol{z}}_{\mathrm{c}}\right) $ and $ {r}_{\mathrm{y}}\left(\boldsymbol{y}|{\boldsymbol{z}}_{\mathrm{y}}\right) $ . The resulting hybrid generative model $ {p}_{\boldsymbol{\theta}}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ can be coupled with variational posteriors $ {q}_{\phi_{\mathrm{x}}}\left({\boldsymbol{z}}_{\mathrm{x}}|\boldsymbol{x}\right){q}_{\phi_{\mathrm{c}}}\left({\boldsymbol{z}}_{\mathrm{c}}|\boldsymbol{x}\right){q}_{\phi_{\mathrm{y}}}\left({\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right) $ to yield an architecture similar to VAE. This is the architecture derived in Section 3, without the additional constraints. It should be noted that the assumption of an additive structure for the discrepancy term $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ and the uncertainty term $ {\boldsymbol{\epsilon}}_{\mathrm{x}} $ presented in Equation (2.2) is suitable for many physical systems and is commonly employed in hybrid models (Cross et al., Reference Cross, Gibson, Jones, Pitchforth, Zhang, Rogers, Cury, Ribeiro, Ubertini and Todd2022). We use it without loss of generality with the aim of promoting clarity and interpretability. Depending on domain knowledge regarding the problem at hand, a multiplicative or other form can also be specified. The implications of the additivity assumption are further discussed in Sections 3.2.3 and 6.2.

Neither the additive structure of the hybrid physics-based and data-driven model, nor the specified parametrization of the residual term $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ ensure that the known physics will be utilized by the model, or that $ {\boldsymbol{z}}_{\mathrm{x}} $ will be physically meaningful. Without further constraints, the model can learn combinations of arbitrary predictions from the physics-based and data-driven components $ f\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ and $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ that sum to an accurate prediction. To see why, it is sufficient to consider the form of the objective given in Equation (2.1). Both the encoder and decoder are aligned in the task of maximizing the reconstruction term $ {\unicode{x1D53C}}_{q_{\boldsymbol{\phi}}\left(\boldsymbol{z}|\boldsymbol{x}\right)}\left[\log {p}_{\boldsymbol{\theta}}\left(\boldsymbol{x}|\boldsymbol{z}\right)\right] $ , and the data-driven components of the model will account for discrepancies between the physics-based model prediction and measurements up to some level of noise. This results in an entangled representation, where the data-driven components override the known physics and the physics-grounded latent variables loose their physical meaning.

This issue is illustrated in Figure 2 using the beam case study shown in Figure 1(a). Further details of the case study are provided in Section 5.1. In this example, a VAE trained on a dataset $ \mathbf{\mathcal{D}}={\left\{\left({\boldsymbol{x}}_i,{\boldsymbol{c}}_i,{\boldsymbol{y}}_i\right)\right\}}_{i=1}^N $ is evaluated on a new set of input measurements $ \boldsymbol{x} $ , generated from the ground truth data generating process by linearly varying the position of the load $ {x}_{\mathrm{F}} $ . It can be seen that the effect of this variation on the measured response is largely captured by the data-driven component of the decoder, which overrides the known physics, despite the fact that the physics-based model includes the load position as an input parameter. The extent to which the data-driven components override the known physics can be inconsistent and hard to predict, and will depend on the neural network architectures and the physics of the problem at hand.

Figure 2. Demonstration of the data-driven component of the decoder $ {g}_{\boldsymbol{\theta}}\left({\mathbf{z}}_{\mathrm{c}},{\mathbf{z}}_{\mathrm{y}}\right) $ overriding the physics-based model $ f\left({\mathbf{z}}_{\mathrm{x}}\right) $ . The effect of varying the position of the load $ {x}_{\mathrm{F}} $ should be described by the known physics, but is instead captured by the data-driven components.

The results presented in Figure 2 are obtained under the assumption of a factorized variational posterior. Ideally, a single encoder with shared parameters $ {q}_{\phi}\left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right) $ would be used for the three subsets of the latent space $ {\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ . However, using a shared encoder leads to further degradation of the performance as noted by Ilse et al., Reference Ilse, Tomczak, Louizos, Welling, Arbel, Ben Ayed, Bruijne, Descoteaux, Lombaert and Pal2020. This example demonstrates how the interaction between the physics-based and data-driven components, which determines the resulting learned representation, depends on the capacity and flexibility of the individual machine learning components. Standard VAE offer no mechanism to control this interaction, and are therefore unable to guarantee that the physics will be utilized correctly in the presence of flexible NN-based decoder components.

3.1. Detailed description of the model

It is assumed that three generative factors, the underlying physics of the structure, the domain, and the class, contribute to the measured response $ \boldsymbol{x} $ . Conversely, the latent variables are partitioned into subsets $ \left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)\in \boldsymbol{z} $ . It is emphasized that the separation of the latent variables is only semantic and used for clarity. In practice, they can be the output of a single encoder with shared parameters. The latent variables are the input to the hybrid probabilistic decoder $ {p}_{\theta}\left(\boldsymbol{x}|\boldsymbol{z}\right) $ , wherein a NN-based function $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ accounts for discrepancies between the measured response $ \boldsymbol{x} $ and physics-based model prediction $ f\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ . To ensure that information relevant to the domain and class is encoded in the corresponding subsets $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ , we utilize two auxiliary decoders $ {r}_{\psi_{\mathrm{c}}}\left(\boldsymbol{c}|{\boldsymbol{z}}_c\right) $ and $ {r}_{\psi_{\mathrm{y}}}\left(\boldsymbol{y}|{\boldsymbol{z}}_y\right) $ . The latent variables $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ are assigned conditional prior distributions $ {p}_{\theta_{\mathrm{c}}}\left({\boldsymbol{z}}_{\mathrm{c}}|\boldsymbol{c}\right) $ and $ {p}_{\theta_{\mathrm{y}}}\left({\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{y}\right) $ respectively, while the physics-grounded latent variables $ {\boldsymbol{z}}_{\mathrm{x}} $ are assigned a distribution $ p\left({\boldsymbol{z}}_{\mathrm{x}}\right) $ based on the available prior knowledge. A schematic illustration of the architecture is provided in Figure 3(a).

Figure 3. a) Schematic diagram illustrating the components of the model and the encoder-decoder architecture, and b) Detailed structure of the dependencies in the generative and inference models.

To minimize the reconstruction error, the encoder tends to maximize the information in the posterior distribution over $ \boldsymbol{z} $ that can be used to predict $ \boldsymbol{x} $ , $ \boldsymbol{c} $ and $ \boldsymbol{y} $ , subject to the regularization imposed by the prior distribution. For $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ , this information includes features from the input signal $ \boldsymbol{x} $ that are predictive of $ \boldsymbol{c} $ and $ \boldsymbol{y} $ , but also irrelevant features that are only predictive of $ \boldsymbol{x} $ . These features can include systematic errors stemming from partial knowledge of the physics, and the influence of unknown confounding factors in the measurements. This superfluous information (Federici et al., Reference Federici, Dutta, Forré, Kushman and Akata2020), i.e. information in $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ that is not predictive of $ \boldsymbol{c} $ and $ \boldsymbol{y} $ , can enable the data-driven component of the decoder to override the known physics. Motivated by this observation we aim to simultaneously maximize the information in $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ that is predictive of $ \boldsymbol{c} $ and $ \boldsymbol{y} $ , while minimizing the information that is predictive of $ \boldsymbol{x} $ . This trade-off can be formalized in terms of the mutual information (MI), a measure of the dependence between two random variables (Cover and Thomas, Reference Cover and Thomas2006). Denoting the MI between $ \boldsymbol{x} $ and $ \boldsymbol{z} $ for an encoder parametrized by $ \boldsymbol{\phi} $ as $ {I}_{\phi}\left(\boldsymbol{x};\boldsymbol{z}\right) $ , and introducing the trade-off parameters $ {\lambda}_{\mathrm{c}},{\lambda}_{\mathrm{y}} $ , we define the following relaxed Lagrangian objectives:

(3.1) $$ {\displaystyle \begin{array}{c}{\mathrm{\mathcal{L}}}_{\mathrm{y}}\left(\boldsymbol{\phi}; {\lambda}_{\mathrm{y}}\right)={I}_{\phi}\left(\boldsymbol{y};{\boldsymbol{z}}_{\mathrm{y}}\right)-{\lambda}_{\mathrm{y}}{I}_{\phi}\left(\boldsymbol{x};{\boldsymbol{z}}_{\mathrm{y}}\right)\\ {}{\mathrm{\mathcal{L}}}_{\mathrm{c}}\left(\boldsymbol{\phi}; {\lambda}_{\mathrm{c}}\right)={I}_{\phi}\left(\boldsymbol{c};{\boldsymbol{z}}_{\mathrm{c}}\right)-{\lambda}_{\mathrm{c}}{I}_{\phi}\left(\boldsymbol{x};{\boldsymbol{z}}_{\mathrm{c}}\right)\end{array}} $$

The quantities described in Equation (3.1) are optimized indirectly through a latent bottleneck structure (Alemi et al., Reference Alemi, Fischer, Dillon and Murphy2019; Fischer, Reference Fischer2020; Moyer et al., Reference Moyer, Gao, Brekelmans, Galstyan, Ver Steeg, Bengio, Wallach, Larochelle, Grauman, Cesa-Bianchi and Garnett2018; Tishby et al., Reference Tishby, Pereira and Bialek2000) combined with adversarial training, as described in the following informal sketch. Note that in the inference model shown in Figure 3(b), the latent variables $ {\boldsymbol{z}}_{\mathrm{y}} $ do not depend directly on $ \boldsymbol{y} $ (and analogously for $ \boldsymbol{c} $ ). This is the conditional independence assumption typically used in the Information Bottleneck framework. Intuitively, the encoder is forced to distill the relevant information in $ \boldsymbol{x} $ that is necessary for reconstructing $ \boldsymbol{c} $ and $ \boldsymbol{y} $ into the latent variables $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ . This results in the maximization of the $ {I}_{\phi}\left(\boldsymbol{c};{\boldsymbol{z}}_{\mathrm{c}}\right) $ and $ {I}_{\phi}\left(\boldsymbol{y};{\boldsymbol{z}}_{\mathrm{y}}\right) $ terms in Equation (3.1) during training. The additional requirement of minimizing $ {I}_{\phi}\left(\boldsymbol{x};{\boldsymbol{z}}_{\mathrm{c}}\right) $ and $ {I}_{\phi}\left(\boldsymbol{x};{\boldsymbol{z}}_{\mathrm{y}}\right) $ can be satisfied by introducing a Gradient Reversal Layer (GRL) (Ganin and Lempitsky, Reference Ganin, Lempitsky, Bach and Blei2015) at the input of the data-driven decoder component $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ . During optimization, the gradient signal propagated backwards from $ {g}_{\boldsymbol{\theta}}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ to the encoder is scaled by $ -\lambda $ , while the forward pass remains unchanged. Therefore, the GRL can be thought of as a pseudo function $ {R}_{\lambda}\left(\boldsymbol{z}\right) $ such that $ {R}_{\lambda}\left(\boldsymbol{z}\right)=\boldsymbol{z} $ and $ \frac{\mathrm{d}{R}_{\lambda }}{\mathrm{d}\boldsymbol{z}}=-\lambda \boldsymbol{I} $ . Positive values of $ \lambda $ correspond to adversarial training. Conversely, negative values make the training “collaborative”. The absolute value of $ \lambda $ determines the strength of the adversarial or collaborative objective, with larger values corresponding to a stronger regularization effect. By turning the decoder $ {p}_{\theta}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ into an adversary, the GRL penalizes information in $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ that contributes to the reconstruction of $ \boldsymbol{x} $ , biasing the encoder towards representations that are minimally informative about $ \boldsymbol{x} $ .

The full structure of the model is shown in Figure 3(b). The variational lower bound can be obtained by considering the marginal likelihood over observed variables as shown in Equation (3.2).

(3.2) $$ {\displaystyle \begin{array}{c}\mathrm{\mathcal{L}}\left(\boldsymbol{\theta}, \phi; \boldsymbol{x},\boldsymbol{c},\boldsymbol{y}\right)={\unicode{x1D53C}}_{q_{\phi}\left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right)}\left[\log \frac{p_{\theta}\left(\boldsymbol{x},\boldsymbol{c},\boldsymbol{y},{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)}{q_{\phi}\left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right)}\right]\\ {}\hskip1.3em \le \log {p}_{\theta}\left(\boldsymbol{x},\boldsymbol{c},\boldsymbol{y}\right)\end{array}} $$

Rearranging the terms in Equation (3.2), noting that the generative model factorizes as $ p\left(\boldsymbol{x},\boldsymbol{c},\boldsymbol{y},{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)={p}_{\theta}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)p\left({\boldsymbol{z}}_{\mathrm{x}}\right){p}_{\theta_{\mathrm{c}}}\left({\boldsymbol{z}}_{\mathrm{c}}|\boldsymbol{c}\right){p}_{\theta_{\mathrm{y}}}\left({\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{y}\right)p\left(\boldsymbol{c}\right)p\left(\boldsymbol{y}\right) $ , yields the following expression for the lower bound:

(3.3) $$ {\displaystyle \begin{array}{c}\mathrm{\mathcal{L}}\left(\boldsymbol{\theta}, \boldsymbol{\phi}; \boldsymbol{x},\boldsymbol{c},\boldsymbol{y}\right)={\unicode{x1D53C}}_{q_{\phi}\left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right)}\left[\log {p}_{\theta}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)\right]\\ {}-{D}_{\mathrm{KL}}\left({q}_{\phi}\left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right)\;\left\Vert\;p\left({\boldsymbol{z}}_{\mathrm{x}}\right){p}_{\theta_{\mathrm{c}}}\right({\boldsymbol{z}}_{\mathrm{c}}|\boldsymbol{c}\left){p}_{\theta_{\mathrm{y}}}\right({\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{y}\Big)\right)\end{array}} $$

Including the auxiliary tasks and additional regularization hyperparameters commonly used in representation learning, we rewrite the loss function as:

(3.4) $$ {\displaystyle \begin{array}{c}\mathrm{\mathcal{L}}\left(\boldsymbol{\theta}, \boldsymbol{\phi}, \boldsymbol{\psi}; \boldsymbol{x},\boldsymbol{c},\boldsymbol{y}\right)={\unicode{x1D53C}}_{q_{\phi \left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right)}}\left[{\alpha}_{\mathrm{x}}\log {p}_{\theta}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right)+{\alpha}_{\mathrm{c}}\log {r}_{\psi_{\mathrm{c}}}\Big(\boldsymbol{c}|{\boldsymbol{z}}_{\mathrm{c}}\left)+{\alpha}_{\mathrm{y}}\log {r}_{\psi_{\mathrm{y}}}\right(\boldsymbol{y}|{\boldsymbol{z}}_{\mathrm{y}}\Big)\right]\\ {}-\beta {D}_{\mathrm{KL}}\left({q}_{\phi}\left({\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{x}\right)\left\Vert p\left({\boldsymbol{z}}_{\mathrm{x}}\right){p}_{\theta_{\mathrm{c}}}\right({\boldsymbol{z}}_{\mathrm{c}}|\boldsymbol{c}\left){p}_{\theta_{\mathrm{y}}}\right({\boldsymbol{z}}_{\mathrm{y}}|\boldsymbol{y}\Big)\right)\end{array}} $$

The loss function described in Equation (3.4)) includes additional regularization hyperparameters that can be used to balance the contribution of different terms. These are included for completeness, and are not used in the experiments described in Section 5. A scaling factor $ \beta >0 $ on the KLD is commonly included in the ELBO as a means of adjusting the strength of the regularization imposed by the KLD term, and to control the capacity of the probabilistic encoder. It is often beneficial to begin training with $ \beta =0 $ and gradually increase it to $ \beta =1 $ using an annealing scheme such as the one proposed by Bowman et al., Reference Bowman, Vilnis, Vinyals, Dai, Jozefowicz and Bengio2016. Annealing $ \beta $ can prevent the model from getting stuck in local minima of the KLD, and the posterior distribution from degenerating to the prior distribution. Conversely, setting $ \beta >1 $ can promote unsupervised disentanglement (Higgins et al., Reference Higgins, Matthey, Pal, Burgess, Glorot, Botvinick, Mohamed and Lerchner2016). The impact of $ \beta $ is extensively discussed in the relevant literature, provided in Section 4. Additionally, the log-likelihood function $ \log {p}_{\theta}\left(\boldsymbol{x}|{\boldsymbol{z}}_{\mathrm{x}},{\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ , and auxiliary decoders $ \log {r}_{\psi_{\mathrm{c}}}\left(\boldsymbol{c}|{\boldsymbol{z}}_{\mathrm{c}}\right) $ and $ \log {r}_{\psi_{\mathrm{y}}}\left(\boldsymbol{y}|{\boldsymbol{z}}_{\mathrm{y}}\right) $ are assigned weights $ {\alpha}_{\mathrm{x}} $ , $ {\alpha}_{\mathrm{c}} $ and $ {\alpha}_{\mathrm{y}} $ respectively to allow for balancing the relative strength of these terms (see e.g. (Ilse et al., Reference Ilse, Tomczak, Louizos, Welling, Arbel, Ben Ayed, Bruijne, Descoteaux, Lombaert and Pal2020; Joy et al., Reference Joy, Schmon, Torr, Siddharth and Rainforth2022; Sun et al., Reference Sun, Pears and Gu2022). It is important to note that using values other than unity for $ {\alpha}_{\mathrm{x}},{\alpha}_{\mathrm{c}},{\alpha}_{\mathrm{y}} $ and $ \beta $ can have a significant impact on the interpretation of the ELBO and the inferred posterior distribution. Details of the implementation can be found in Appendix A.

3.2. Discussion of the approach

The latent bottleneck structure and GRL have several important implications for inference, conditional generation, and uncertainty quantification. These are discussed here, and demonstrated through the case studies presented in Section 5.

3.3. Description of the quantitative assessment approach

Quantitatively assessing disentanglement is a challenging problem and several metrics have been proposed (Chen et al., Reference Chen, Li, Grosse, Duvenaud, Bengio, Wallach, Larochelle, Grauman, Cesa-Bianchi and Garnett2018; Higgins et al., Reference Higgins, Matthey, Pal, Burgess, Glorot, Botvinick, Mohamed and Lerchner2016; Kim and Mnih, Reference Kim, Mnih, Dy and Krause2018). This difficulty can be partially attributed to the lack of a consistent definition of disentanglement (Locatello et al., Reference Locatello, Bauer, Lucic, Rätsch, Gelly, Schölkopf and Bachem2019). In practice, the degree of disentanglement achieved by a model is often evaluated based on subjective expectations stemming from domain knowledge (Vowels et al., Reference Vowels, Camgoz and Bowden2019). Our proposed approach aims to achieve “one-way” disentanglement and is conceptually more akin to techniques that temper the influence of misspecified model components (Carmona and Nicholls, Reference Carmona, Nicholls, Chiappa and Calandra2020; Yu et al., Reference Yu, Nott and Smith2022). Intuitively, variations in generative factors that are described by the known physics should not affect the data-driven subsets of the latent variables $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ , while variations in generative factors not included in the known physics should still be preferentially captured by the physics-grounded subset of the latent variables. The degree to which this is achieved can be evaluated by comparing the amount of information captured by each subset $ {\boldsymbol{z}}_{\mathrm{x}} $ , $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ about a specified generative factor. When a subset of the latent variables is informative about a generative factor, it should be possible to train a regressor to predict the value of the generative factor from samples drawn from this subset of the latent variables. Based on this, we propose the following procedure to assess the amount of information about a given generative factor that is encoded in a subset of the latent variables for the trained model:

1. 1. Draw two sets of samples of generative factors $ {\left\{\left({\boldsymbol{s}}_{\mathrm{x}}^{(i)},{\boldsymbol{s}}_{\mathrm{c}}^{(i)},{\boldsymbol{s}}_{\mathrm{y}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{train}}} $ and $ {\left\{\left({{\boldsymbol{s}}^{\prime}}_{\mathrm{x}}^{(i)},{{\boldsymbol{s}}^{\prime}}_{\mathrm{c}}^{(i)},{{\boldsymbol{s}}^{\prime}}_{\mathrm{y}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{test}}} $ from $ {p}_{\mathrm{gt}}\left({\boldsymbol{s}}_{\mathrm{x}},{\boldsymbol{s}}_{\mathrm{c}},{\boldsymbol{s}}_{\mathrm{y}}\right) $ and generate two datasets $ D={\left\{{\boldsymbol{x}}_i\right\}}_{i=1}^{N_{\mathrm{train}}} $ and $ {D}^{\prime }={\left\{{\boldsymbol{x}}_i^{\prime}\right\}}_{i=1}^{N_{\mathrm{test}}} $ of response measurements from the ground truth generative process.

2. 2. Draw a single sample from each of the approximate posterior distributions $ {\boldsymbol{z}}_i\sim {q}_{\phi}\left({\boldsymbol{z}}_i|{\boldsymbol{x}}_i\right) $ and $ {\boldsymbol{z}}_i^{\prime}\sim {q}_{\phi}\left({\boldsymbol{z}}_i^{\prime }|{\boldsymbol{x}}_i^{\prime}\right) $ for each $ {\boldsymbol{x}}_i\in D $ and $ {\boldsymbol{x}}_i^{\prime}\in {D}^{\prime } $ respectively, using the trained model. This yields two sets of samples from the latent variables $ {\left\{\left({\boldsymbol{z}}_{\mathrm{x}}^{(i)},{\boldsymbol{z}}_{\mathrm{c}}^{(i)},{\boldsymbol{z}}_{\mathrm{y}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{train}}} $ , and $ {\left\{\left({{\boldsymbol{z}}^{\prime}}_{\mathrm{x}}^{(i)},{{\boldsymbol{z}}^{\prime}}_{\mathrm{c}}^{(i)},{{\boldsymbol{z}}^{\prime}}_{\mathrm{y}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{test}}} $ .

3. 3. Train a regressor to predict the value of each set of generative factors $ {\left\{{s}_j^{(i)}\right\}}_{i=1}^{N_{\mathrm{train}}} $ from each subset $ {\left\{{\boldsymbol{z}}_{\mathrm{x}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{train}}} $ , $ {\left\{{\boldsymbol{z}}_{\mathrm{c}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{train}}} $ and $ {\left\{{\boldsymbol{z}}_{\mathrm{y}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{train}}} $ , for $ j=1,\dots, {N}_f $ , where $ {N}_f $ is the number of generative factors. This process yields $ 3\times {N}_{\mathrm{f}} $ regressors.

4. 4. Compute the $ {R}^2 $ value between each subset $ {\left\{{{\boldsymbol{z}}^{\prime}}_{\mathrm{x}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{test}}} $ , $ {\left\{{{\boldsymbol{z}}^{\prime}}_{\mathrm{c}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{test}}} $ and $ {\left\{{{\boldsymbol{z}}^{\prime}}_{\mathrm{y}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{test}}} $ and each set of generative factors $ {\left\{{s^{\prime}}_j^{(i)}\right\}}_{i=1}^{N_{\mathrm{test}}} $ using the corresponding trained regression model.

This procedure yields $ {N}_f $ sets of pair-wise $ {R}^2 $ values $ {\left\{{R}_{{\boldsymbol{z}}_{\mathrm{x}}\to {s}_j}^2,{R}_{{\boldsymbol{z}}_{\mathrm{c}}\to {s}_j}^2,{R}_{{\boldsymbol{z}}_{\mathrm{y}}\to {s}_j}^2\right\}}_{j=1}^{N_f} $ with each of the $ {N}_f $ sets corresponding to a single generative factor. A more informative subset of the latent variables should yield a more accurate regressor than an uninformative subset, and therefore also a higher $ {R}^2 $ value. It is emphasized that the metric described here is only intended to be a surrogate quantity for the amount of information encoded in each subset of the latent variables, and not a metric of disentanglement. Furthermore, this metric requires access to the ground truth distribution and data generating process, and is therefore not generally applicable.

Case study objectives

Each case study addresses a different set of challenges. The beam case study demonstrates that the proposed approach preferentially utilizes the known physics, yielding an interpretable and parsimonious representation of the physical system. The oscillator case study highlights how the adversarial training can prevent the model from learning arbitrary components of the response in a data-driven manner, and investigates the impact of the GRL hyperparameter. Finally, the bridge case study demonstrates the feasibility of using the model for damage detection in a more complex synthetic case, and compares the performance to that of existing data-driven approaches. It is noted that the case studies are only meant as didactic examples, intended to elaborate the issues with combining physics-based and data-driven components in VAE, provide intuition about the interaction between these components, and demonstrate the behavior of the model. Therefore, emphasis is placed on clarity rather than realism.

Physics-based models

Three separate physics-based models are considered for every case study: A high-fidelity simulator, a full model, and a nominal model. The simulator is an accurate but generally computationally expensive model of the physical system, typically in the form of a finite element (FE) model, used to train the full and nominal models for each case study. The full model is a computationally efficient surrogate model of the ground truth data generating process: one or more structures with varying physical characteristics, subject to operational and environmental conditions, damage and degradation. To produce the training dataset for the full model, the simulator is evaluated on a set of generative factors $ {\left\{\left({\boldsymbol{s}}_{\mathrm{x}}^{(i)},{\boldsymbol{s}}_{\mathrm{c}}^{(i)},{\boldsymbol{s}}_{\mathrm{y}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{full}}} $ , sampled uniformly and independently from prescribed ranges of values. The ranges are chosen to provide sufficient coverage over the support of the corresponding ground truth distribution $ {p}_{\mathrm{gt}}\left({\boldsymbol{s}}_{\mathrm{x}},{\boldsymbol{s}}_{\mathrm{c}},{\boldsymbol{s}}_{\mathrm{y}}\right) $ . The full model is then obtained by fitting a NN-based surrogate to the dataset composed of $ {N}_{\mathrm{full}} $ input-output pairs $ {D}_{\mathrm{full}}={\left\{\left({\boldsymbol{s}}_{\mathrm{x}}^{(i)},{\boldsymbol{s}}_{\mathrm{c}}^{(i)},{\boldsymbol{s}}_{\mathrm{y}}^{(i)},{\boldsymbol{x}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{full}}} $ obtained from the simulator. Using a NN as the forward model for the data generating process enables the efficient visualization of the latent space and the reconstructions generated by the VAE for different inputs, simplifies the generation of test data to evaluate the performance of the VAE, and makes it possible to account for randomness in the hyperparameter initialization and data generation by averaging results over multiple runs with i.i.d. datasets. The nominal model corresponds to the available incomplete representation of the physics of the system under investigation. When an analytical expression describing the partially known physics is available, this is used as the nominal model. Alternatively, the nominal model is built by training a NN-based surrogate on a limited dataset $ {D}_{\mathrm{nom}}={\left\{\left({\boldsymbol{s}}_{\mathrm{x}}^{(i)},{\boldsymbol{x}}^{(i)}\right)\right\}}_{i=1}^{N_{\mathrm{nom}}} $ , obtained by evaluating the simulator only on the physics-based subset of the generative factors $ {\left\{{\boldsymbol{s}}_{\mathrm{x}}^{(i)}\right\}}_{i=1}^{N_{\mathrm{nom}}} $ , while $ {\boldsymbol{s}}_{\mathrm{c}} $ and $ {\boldsymbol{s}}_{\mathrm{y}} $ are set to a constant reference value corresponding to the nominal condition of the structure.

Synthetic data generation

The VAE is trained and validated on a dataset composed of $ {N}_{\mathrm{total}}={N}_{\mathrm{train}}+{N}_{\mathrm{val}} $ triplets of observables $ \mathbf{\mathcal{D}}={\left\{\left({\boldsymbol{x}}_i,{\boldsymbol{c}}_{\mathbf{i}},{\boldsymbol{y}}_{\mathbf{i}}\right)\right\}}_{i=1}^{N_{\mathrm{total}}} $ . This dataset is generated by first drawing samples of the generative factors from the ground truth distribution $ \left({\boldsymbol{s}}_{\mathrm{x}}^{(i)},{\boldsymbol{s}}_{\mathrm{c}}^{(i)},{\boldsymbol{s}}_{\mathrm{y}}^{(i)}\right)\sim {p}_{\mathrm{gt}}\left({\boldsymbol{s}}_{\mathrm{x}},{\boldsymbol{s}}_{\mathrm{c}},{\boldsymbol{s}}_{\mathrm{y}}\right) $ for $ i=1,\dots, {N}_{\mathrm{total}} $ , applying a set of deterministic transformations, and subsequently adding i.i.d. samples of zero-mean Gaussian white noise $ {\boldsymbol{\epsilon}}_{\mathrm{x}}\sim \mathbf{\mathcal{N}}\left(\mathbf{0},{\sigma}_{\mathrm{x}}^2\boldsymbol{I}\right) $ , $ {\boldsymbol{\epsilon}}_{\mathrm{c}}\sim \mathbf{\mathcal{N}}\left(\mathbf{0},{\sigma}_{\mathrm{c}}^2\boldsymbol{I}\right) $ and $ {\boldsymbol{\epsilon}}_{\mathrm{y}}\sim \mathbf{\mathcal{N}}\left(\mathbf{0},{\sigma}_{\mathrm{y}}^2\boldsymbol{I}\right) $ . Denoting the full model as $ {h}_{\mathrm{x}}\left(\cdot \right) $ , the response observables are obtained as $ {\boldsymbol{x}}_i={h}_{\mathrm{x}}\left({\boldsymbol{s}}_{\mathrm{x}}^{(i)},{\boldsymbol{s}}_{\mathrm{c}}^{(i)},{\boldsymbol{s}}_{\mathrm{y}}^{(i)}\right)+{\boldsymbol{\epsilon}}_{\mathrm{x}} $ . The domain and class observables are obtained as $ {\boldsymbol{c}}_i={h}_{\mathrm{c}}\left({\boldsymbol{s}}_{\mathrm{c}}^{(i)}\right)+{\boldsymbol{\epsilon}}_{\mathrm{c}} $ and $ {\boldsymbol{y}}_i={h}_{\mathrm{y}}\left({\boldsymbol{s}}_{\mathrm{y}}^{(i)}\right)+{\boldsymbol{\epsilon}}_{\mathrm{y}} $ respectively. This procedure is illustrated in Figure 4.

Figure 4. Illustration of the procedure used to obtain the full and nominal physics-based models (left), and to generate the datasets used in the case studies (right).

Implementation details

For all the case studies presented in this section, $ {h}_{\mathrm{c}} $ and $ {h}_{\mathrm{y}} $ are taken as the identity function for simplicity. Furthermore we use $ {N}_{\mathrm{train}}=1024 $ and $ {N}_{\mathrm{val}}=512 $ , and consider no other regularization except for the GRL, i.e. $ \beta ={\alpha}_{\mathrm{x}}={\alpha}_{\mathrm{c}}={\alpha}_{\mathrm{y}}=1.0 $ . Unless stated otherwise, the number of the domain and class latent variables are taken to be twice the number of domain and class generative factors. The intention behind this choice is to avoid biasing the model towards a disentangled representation by matching the number of latent variables to the ground truth generative factors, ensuring that any disentanglement in the learned representation is not a consequence of limited latent space capacity. To enable the formulation of problems with bounded latent variables and to ensure a stable optimization procedure, the physics-grounded latent variables are obtained through a sequence of invertible transformations applied to the encoder output, mapping samples from an unbounded base latent space to the target latent space. All physics-grounded latent variables are constrained to lie within ranges that ensure consistency with the underlying physics of each system.

Visualization

A particularly useful tool for assessing the learned representation is to “traverse” the latent space and the space of reconstructions of the VAE (Higgins et al., Reference Higgins, Matthey, Pal, Burgess, Glorot, Botvinick, Mohamed and Lerchner2016; Kim and Mnih, Reference Kim, Mnih, Dy and Krause2018). This can be achieved by generating synthetic data while interpolating over a specified generative factor, and setting the remaining generative factors to a constant reference value. The VAE is then evaluated on the generated data, yielding samples from the latent space, realizations of the reconstructed input $ \hat{\boldsymbol{x}} $ , and the mean physics-based and data-driven components $ {\hat{\boldsymbol{x}}}_{\mathrm{p}} $ and $ {\hat{\boldsymbol{x}}}_{\mathrm{d}} $ . Each generative factor is linearly interpolated within the $ {1}^{\mathrm{st}} $ and $ {99}^{\mathrm{th}} $ percentiles of the corresponding ground truth distribution.

5.1. Beam case study

5.1.1. Case study description

The case study consists of a beam with fixed length $ L=1.0 $ m and a point load with magnitude $ F=1.0 $ N acting on an uncertain position $ {x}_{\mathrm{F}} $ along the length of the beam. The material is linear elastic with uncertain Young’s modulus $ E $ , Poisson ratio $ \nu =0.3 $ , area moment of inertia $ I=2\cdot {10}^{-6} $ m $ {}^4 $ and cross-sectional area $ A=2.4\cdot {10}^{-3} $ m $ {}^2 $ . The rotational stiffness of the right-hand side support is temperature dependent, with the dependence modeled as an increase in the rotational stiffness of the support at lower temperatures. The relationship between the temperature and the support rotational stiffness is formulated as $ \log {k}_{\mathrm{r}}=8-\frac{10}{1+{e}^{-T/2}} $ . The beam is subject to variability in the vertical stiffness of the right-hand side support, e.g. due to damage or a deficiency of the support, simulated as a translational spring boundary condition with stiffness $ {k}_{\mathrm{v}} $ . This quantity can span several orders of magnitude, and therefore we parametrize the model using $ \log {k}_{\mathrm{v}} $ instead. The beam is equipped with $ {d}_{\mathrm{x}}=32 $ sensors measuring the vertical displacement, equally spaced along the length as shown in Figure 1(a).

The Young’s modulus $ E $ and the position of the point load $ {x}_{\mathrm{F}} $ are considered as uncertain latent variables, such that $ {\boldsymbol{z}}_{\mathrm{x}}=\left(E,{x}_{\mathrm{F}}\right) $ . It is assumed that the temperature is an observed domain variable such that $ \boldsymbol{c}=(T) $ , and that the vertical spring log-stiffness $ \log {k}_{\mathrm{v}} $ is taken as a class variable representing damage in the structure, such that $ \boldsymbol{y}=\left(\log {k}_{\mathrm{v}}\right) $ . Since the class variable $ \boldsymbol{y} $ represents damage in the structure it will not be quantitatively measurable. In a realistic scenario, observations of the condition of the support on a qualitative scale (e.g. from $ 0 $ representing no damage to $ 5 $ denoting a fully damaged support) might be available. In this example we simplistically consider $ \boldsymbol{y} $ as the ground truth value of $ \log {k}_{\mathrm{v}} $ with some added noise. The variable symbols, units, types, as well as the prior distributions over the physics-grounded latent variables and the ground truth distributions of the generative factors used to generate the training data are summarized in Table 1. We additionally provide a reference value which is used to produce the figures as discussed in Section 5. To ensure physical consistency and to avoid numerical issues, the Young’s modulus is truncated below a small positive value, and the load position $ {x}_{\mathrm{F}} $ is restricted to the range $ \left(0,1\right) $ .

Table 1. Summary of generative factors and the corresponding ground truth and prior distributions

A partial description of the physics is available, in the form of an analytical expression for the vertical deflection of a simply supported Euler-Bernoulli beam with a point load acting at $ {x}_{\mathrm{F}} $ :

(5.1) $$ w(x)=\left\{\begin{array}{ll}\frac{Pbx\left({L}^2-{b}^2-{x}^2\right)}{6 LEI},& 0\le x\le {x}_{\mathrm{F}}\\ {}\frac{Pbx\left({L}^2-{b}^2-{x}^2\right)}{6 LEI}+\frac{P{\left(x-{x}_{\mathrm{F}}\right)}^3}{6 EI},& {x}_{\mathrm{F}}<x\le L\end{array}\right. $$

where $ b=L-{x}_{\mathrm{F}} $ , and the non-bold $ x $ refers to the position along the beam. This nominal model represents the beam in the undamaged condition at a reference temperature, and is directly incorporated in the physics-based branch of the VAE decoder.

Following the procedure described in Section 5, the full model (trained on input-output pairs from an FE-based simulator) is used to produce synthetic data by first drawing samples of the input parameters from the ground truth distribution, and subsequently contaminating the resulting model predictions with zero-mean Gaussian white noise with standard deviation $ {\sigma}_{\mathrm{x}}=0.02 $ m. The dataset used to train the VAE is composed of $ {N}_{\mathrm{train}} $ measurements of the beam displacement $ \boldsymbol{x}={\left\{{\boldsymbol{x}}_i\right\}}_{i=1}^{N_{\mathrm{train}}} $ where each element $ {\boldsymbol{x}}_i $ is a vector of length $ {d}_{\mathrm{x}}=32 $ . The domain and class observables $ \boldsymbol{c}={\left\{{\boldsymbol{c}}_i\right\}}_{i=1}^{N_{\mathrm{train}}} $ and $ \boldsymbol{y}={\left\{{\boldsymbol{y}}_i\right\}}_{i=1}^{N_{\mathrm{train}}} $ are obtained as the ground truth values used to generate the dataset, with the addition of i.i.d. samples of Gaussian white noise with standard deviations of $ {\sigma}_c={\sigma}_y=0.02 $ _, respectively. The dimensionality of the domain and class latent variables is taken as $ {d}_{z_{\mathrm{c}}}={d}_{z_{\mathrm{y}}}=2 $ .

5.1.2. Qualitative assessment of disentanglement

After training, the disentanglement between physics-grounded and data-driven components is qualitatively assessed by examining the latent space and samples of the reconstructed response of the beam. The predicted physics-based $ {\hat{\boldsymbol{x}}}_{\mathrm{p}} $ and data-driven $ {\hat{\boldsymbol{x}}}_{\mathrm{d}} $ components, as well as the combined prediction $ \hat{\boldsymbol{x}} $ are shown in Figure 5. It can be observed that the data-driven component of the reconstruction $ {\hat{\boldsymbol{x}}}_{\mathrm{d}} $ (middle row) is invariant to changes in $ E $ and $ {x}_F $ contributing only a constant deformed shape to the total predicted response. On the other hand, the physics-based model captures variability in both the physics-grounded generative factors $ {\boldsymbol{s}}_{\mathrm{x}}=\left(E,{x}_{\mathrm{F}}\right) $ , but also domain and class generative factors $ {\boldsymbol{s}}_{\mathrm{c}}=(T) $ and $ {\boldsymbol{s}}_{\mathrm{y}}=\left(\log {k}_{\mathrm{v}}\right) $ . In contrast to the behavior of the unconstrained VAE, presented in Section 2.4, here the model preferentially utilizes the known physics. Only the variability in the measured response due to $ \log {k}_{\mathrm{v}} $ and $ T $ that can not be captured by the physics-based model is accounted for by the data-driven part of the decoder, indicating that the model can disentangle components of the response that can be attributed to the known physics from those that cannot. A key aspect of the adversarial training is the degree to which it allows interaction between the physics-based and data-driven components of the prediction. In this case study, the additional displacement of the beam due to the reduced vertical stiffness of the right-hand side support will also depend on the load position $ {x}_{\mathrm{F}} $ . More positive values of $ \lambda $ tend to prevent the model from capturing this interaction, whereas more negative values enable it but may result in the data-driven components overriding the known physics.

Figure 5. Mean prediction and $ \pm 2\sigma $ uncertainty bounds for the physics-based $ {\hat{\mathbf{x}}}_{\mathrm{p}} $ and data-driven $ {\hat{\mathbf{x}}}_{\mathrm{d}} $ components, and combined prediction $ \hat{\mathbf{x}} $ while traversing the generative factors. The input response measurements are denoted as dots in the bottom row.

To further highlight the impact of the GRL, the latent space traversals of the unconstrained model and the model trained adversarially are compared in Figure 6. Without adversarial training, the domain latent variables $ {\boldsymbol{z}}_{\mathrm{c}} $ encode the variability in the load position $ {x}_{\mathrm{F}} $ as shown in Figure 6(a), providing the data-driven decoder components with the information needed to reconstruct this component of the measured response and resulting in an entangled representation, as discussed in Section 2.4. In contrast, when $ \lambda =1/256 $ the adversarial training results in a posterior distribution over $ {\boldsymbol{z}}_{\mathrm{c}} $ that is invariant to changes in $ {x}_{\mathrm{F}} $ . Instead, the variability is captured by the corresponding physics-grounded latent variable, as shown in Figure 6(b), indicating disentanglement of the physics-grounded and domain generative factors. The results shown previously suggest that the latent bottleneck architecture and GRL regularization result in a sparse and parsimonious representation of the physical system, and can yield domain and class latent variables that are invariant to changes in the underlying physics. The influence of the GRL hyperparameter $ \lambda $ is further investigated in the oscillator case study presented below.

Figure 6. Visualizations of the VAE latent space during traversal of the generative factors $ {x}_{\mathrm{F}} $ and $ \log {k}_{\mathrm{v}} $ . Each column corresponds to variation of a single generative factor, and each row shows the marginal approximate posterior distribution of a single latent variable.

5.2. Oscillator case study

5.2.1. Case study description

This example demonstrates how the adversarial training prevents the model from compensating for all discrepancies between the physics-based model predictions and measurements. Suppose that a mass-spring-dashpot system undergoes damped harmonic motion, starting from an initial displaced position $ {x}_0 $ , with no external excitation. It is assumed that each experiment is performed under varying temperature $ T $ , which is taken as the domain variable. The temperature affects the spring stiffness through the relationship $ k(T)={k}_{\mathrm{ref}}+{\alpha}_{\mathrm{T}}\left({T}_{\mathrm{ref}}-T\right) $ , with $ {T}_{\mathrm{ref}}=20.0 $ C $ {}^o $ and $ {\alpha}_{\mathrm{T}}=0.01 $ . The reference spring stiffness at $ {T}_{\mathrm{ref}}=20.0 $ C $ {}^o $ is assumed known and equal to $ {k}_{\mathrm{ref}}=1.0 $ N/m. The mass $ m $ is considered unknown and treated as a physics-grounded latent variable to be inferred from data. The viscous damping coefficient $ {c}_{\mathrm{d}} $ is taken to define the class of the system. Finally, it is assumed that the observations are subject to an unknown confounding influence in the form of small random perturbations of the initial displacement $ {x}_0 $ . A summary of the generative factors, the prior distribution and the ground truth distribution is provided in Table 2.

Table 2. Summary of generative factors and the corresponding ground truth and prior distributions

The equation of motion describing the system can be written as:

(5.2) $$ m\frac{{\mathrm{d}}^2x(t)}{\mathrm{d}{t}^2}+{c}_{\mathrm{d}}\frac{\mathrm{d}x(t)}{\mathrm{d}t}+k(T)x(t)=0 $$

A partial description of the physics is available in the form of an analytical solution under the assumption that the initial displacement is $ {x}_0=1.0 $ m, and the initial velocity is $ {\dot{x}}_0=0.0 $ m/s for all experiments, and that there is no damping affecting the motion of the oscillator. Furthermore, it is assumed that the relationship between temperature and stiffness is not known, and the temperature effect is therefore not included in the nominal physics-based model. Under the assumptions described previously, the displacement of the oscillator at time $ t $ can be expressed as:

(5.3) $$ x(t)=\cos \left(\sqrt{\frac{k_{\mathrm{ref}}}{m}}t\right) $$

Each triplet of observations is composed of a noisy displacement time series, and noisy measurements of the viscous damping coefficient $ {c}_{\mathrm{d}} $ and temperature $ T $ , which are considered as class and domain variables respectively such that $ \boldsymbol{c}=(T) $ and $ \boldsymbol{y}=\left({c}_{\mathrm{d}}\right) $ . Training and validation datasets are generated by drawing samples from the ground truth distribution and generating the oscillator displacement time-series using the full model, trained on input-output pairs simulated using the equation of motion shown in Equation (5.2). Each of the measured time-series is a vector of $ 64 $ measurements, equally spaced within a time interval $ t\in \left[0,10\right] $ s. The synthetic response measurements are subsequently contaminated with i.i.d. realizations of Gaussian white noise with standard deviation $ {\sigma}_{\mathrm{x}}=0.01 $ m. The standard deviation of the measurement uncertainty of the domain and class observables are taken as $ {\sigma}_{\mathrm{c}}=0.01 $ and $ {\sigma}_{\mathrm{y}}=0.01 $ respectively. To ensure that the model has sufficient capacity to learn the unknown confounding influence if the adversarial training were not present, the dimensionality of the latent space is specified to be significantly larger than the number of ground truth generative factors. The domain and class latent space dimensions are taken as $ {d}_{z_{\mathrm{c}}}={d}_{z_{\mathrm{y}}}=4 $ .

5.2.2. Model behavior in the presence of unknown confounders

The proposed model with no adversarial training ( $ \lambda =-1 $ ) is trained and evaluated on the synthetic example. The reconstructed response obtained from a traversal of the initial displacement $ {x}_0 $ , shown in Figure 7a, highlights another issue that occurs when combining physics-based and data-driven components in VAE: Although the variability in $ {x}_0 $ can not be accounted for by the physics-based model, and there is no information in the domain or class variables regarding the value of $ {x}_0 $ , the lack of regularization results in a model that is free to capture the components of the measured displacement stemming from the variability in the initial displacement $ {x}_0 $ . Although in this case the effect is benign, in more complex physical systems it can result in the model learning unknown confounding influences in a non-interpretable black-box manner. When the GRL regularization is utilized (Figure 7b), the data-driven encoder is unable to capture the variability in $ {x}_0 $ , depriving the data-driven decoder from the information needed to reconstruct this component of the input measurements.

Figure 7. Physics-based model prediction $ {\hat{\mathbf{x}}}_{\mathrm{p}} $ , data-driven model prediction $ {\hat{\mathbf{x}}}_{\mathrm{d}} $ , and combined prediction $ \hat{\mathbf{x}} $ for varying initial displacement $ {x}_0 $ . With $ \lambda =-1.0 $ (top) the data-driven components in the VAE are free to account for the variability in the initial position. For $ \lambda =1/128 $ (bottom) the model does not learn this component of the response.

The results shown in Figure 7 demonstrate how the unconstrained VAE will compensate for discrepancies between the physics-based model prediction and the measurements caused by unknown confounding influences. The reason why this can be detrimental for the learning task is illustrated in Figure 8a. The unconstrained VAE accounts for the learned confounding influence of $ {x}_0 $ , which can lead to underestimation of the uncertainty over the latent variables and predictions. In contrast, when the model is trained with the adversarial objective, the encoder is prevented from learning a representation of $ {x}_0 $ . The uncertainty stemming from the partial knowledge of the physics, including the unknown influence of the viscous damping and the variability in $ {x}_0 $ , is more accurately accounted for in the reconstructed response, as shown in (Figure 8b). The additional uncertainty can also be attributed to the fact that the mass-spring-dashpot system does not satisfy the additivity assumption described in Section 2.4.

Figure 8. Physics-based model prediction $ {\hat{\mathbf{x}}}_{\mathrm{p}} $ , data-driven model prediction $ {\hat{\mathbf{x}}}_{\mathrm{d}} $ , and combined prediction $ \hat{\mathbf{x}} $ for varying viscous damping coefficient $ {c}_{\mathrm{d}} $ . The data-driven decoder components are prevented from fully accounting for the discrepancies between the physics-based model and measurements, resulting in wider uncertainty bounds for the proposed model.

5.2.3. Quantitative assessment of disentanglement

The trade-off between invariance of the domain and class latent variables to non-domain or class influences and prediction accuracy can be adapted by tuning the GRL hyperparameter $ \lambda $ . A parameter study is performed to assess the impact of different choices for $ \lambda $ on the learned representation. The model is trained for varying values of $ \lambda =\left\{-1,-1/10,-1/100,-1/\mathrm{1000,0,1}/1000,1/\mathrm{100,1}/10,1\right\} $ , and the metric described in Section 3.3 is computed for each trained model using linear regression. The training and testing datasets $ D $ and $ {D}^{\prime } $ are composed of $ 2048 $ samples each. To account for the impact of randomness in the synthetic dataset, neural network parameter initialization, and training procedure, the results for each value of $ \lambda $ are averaged over multiple runs. The values of the metric for each generative factor and subset of the latent space, as a function of $ \lambda $ , are shown in Figure 9.

Figure 9. $ {R}^2 $ value per subset of the latent variables and generative factor as a function of $ \lambda $ , averaged over $ 6 $ runs. The shaded intervals correspond to two standard deviations.

It can be seen that the sign of $ \lambda $ determines the nature of the training procedure, with positive values resulting in adversarial training. For negative values of $ \lambda $ the training becomes collaborative, in the sense that the encoder attempts to find approximate posterior distributions over $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ that are jointly informative about their respective modality as well as the response measurements $ \boldsymbol{x} $ . The magnitude of $ \lambda $ determines the strength of the adversarial or collaborative training. To aid in the interpretation of the results, the behavior of the model is classified into four regimes. When $ \lambda $ approaches $ -1 $ from above, the training is strongly collaborative, and the tasks of minimizing the error in the reconstruction of $ \boldsymbol{x} $ and the prediction of the domain and class variables $ \boldsymbol{c} $ and $ \boldsymbol{y} $ are jointly prioritized. This is reflected by the relatively high scores obtained by the subsets $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ for the generative factors $ m $ and $ {x}_0 $ . For $ \lambda \to {0}^{-} $ , the auxiliary tasks are prioritized over the main task, and the amount of information about $ m $ and $ {x}_0 $ that is encoded in $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ is limited. In this regime the behavior can be characterized as weakly collaborative. Conversely, for small positive values of $ \lambda $ the training becomes weakly adversarial. In this regime the encoder will seek latent codes over $ {\boldsymbol{z}}_{\mathrm{c}} $ and $ {\boldsymbol{z}}_{\mathrm{y}} $ that are uninformative about $ \boldsymbol{x} $ . Further increasing the GRL coefficient such that $ \lambda \to 1 $ yields a strongly adversarial model, and any information that can be used to reconstruct the domain and class variables is heavily weighted against the potential improvement in the reconstruction of $ \boldsymbol{x} $ . In this case the encoder fails to capture the variation in any of the generative factors.

5.3. Bridge case study

5.3.1. Case study description

The final synthetic case study utilizes the two-span bridge benchmark presented in Tatsis and Chatzi, Reference Tatsis and Chatzi2019, illustrated in Figure 1(c). Members of a homogeneous population (Bull et al., Reference Bull, Gardner, Gosliga, Rogers, Dervilis, Cross, Papatheou, Maguire, Campos and Worden2021) of bridges are subjected to controlled loading tests, where a vehicle with known mass and moving at a constant velocity is used to excite the bridge response. The response is obtained as a strain influence line, expressed in parts per thousand (‰), measured by a point strain gauge placed at a distance of $ 5.625 $ m from the start of the bridge, and at a height of $ 0.1 $ m from the bottom of the cross section. Each time-series is composed of $ 64 $ measurements, equally spaced in time $ t\in \left[1,21\right] $ s, where $ t=0 $ s is the moment the vehicle enters the bridge.

The behavior of each bridge is partially determined by the unknown vertical stiffnesses of the supports $ {k}_{\mathrm{v},1} $ , $ {k}_{\mathrm{v},2} $ and $ {k}_{\mathrm{v},3} $ , which are taken to vary between different bridges due to variability in the design, construction and soil conditions. The boundary conditions are known to be symmetric such that $ {k}_{\mathrm{v},1}={k}_{\mathrm{v},3} $ . The base-10 logarithms of the vertical stiffnesses are considered as physics-grounded latent variables. In the horizontal direction, only the left support has a large stiffness, while the rest are unconstrained. It is assumed that the position of the central pier can vary between members of the population by up to $ \pm 1.0 $ m from $ L/2 $ . Furthermore, fluctuations from the prescribed reference vehicle velocity $ {v}_{\mathrm{ref}}=1 $ m/s were observed during the tests that can not be accounted for in the nominal physics-based model. These fluctuations are modeled as a multiplicative term $ {\delta}_{\mathrm{v}} $ such that $ v={\delta}_{\mathrm{v}}\cdot {v}_{\mathrm{ref}} $ . Noisy measurements of the vehicle velocity, and the known pier offsets $ {\delta}_{\mathrm{s}} $ are included as domain variables. It is assumed that the bridge decks are prone to deterioration in a region around the supports. During inspections performed by experts, each bridge is assigned scores $ \boldsymbol{y}=\left({y}_1,{y}_2\right) $ , quantifying the deterioration of the structure near the left and middle supports respectively. Each $ {y}_i $ takes values between zero and one, where zero represents pristine condition and unity corresponds to severe damage of the cross-section at that position. These scores are considered as class observables. Finally, a small variability is considered in the vehicle load such that $ F={\delta}_{\mathrm{F}}\cdot {F}_{\mathrm{ref}} $ , where $ {F}_{\mathrm{ref}}=100 $ kN is the reference load. This variability is caused by deviation in the transverse position of the vehicle, and is considered as an unknown confounding influence. The quantities involved in the case study along with their prior and ground truth distributions are summarized in Table 3.

Table 3. Summary of physics-based, class and domain variables for the two-span bridge case study

Training data for the full and nominal physical models is generated using the FE model of Tatsis and Chatzi, Reference Tatsis and Chatzi2019 as a simulator. The FE model is composed of quadrilateral isoparametric plane stress elements with $ 9 $ nodes each, arranged in a $ 200\times 6 $ grid. The length, width and height are assumed constant and equal to $ L=25.0 $ m, $ w=0.1 $ m and $ h=0.6 $ m for all of the bridges. The material is linear elastic, with Young’s modulus $ E=200 $ GPa, density $ \rho =7850 $ kg/m $ {}^3 $ and Poisson’s ratio $ \nu =0.3 $ . All supports are modeled as linear springs in the vertical direction. The equations of motion are integrated from $ {t}_0=0 $ seconds to $ {t}_1=25 $ seconds with a timestep of $ \mathrm{d}t=0.00025 $ seconds, using an implicit Newmark scheme with parameters $ \gamma =1/2 $ and $ \beta =1/6 $ . The deterioration is modeled as a reduction of the cross section width, ranging from $ 0\% $ (for $ {y}_i=0.0 $ ) to $ 90\% $ (for $ {y}_i=1.0 $ ) in the region spanning $ L/20 $ around the corresponding support. It is noted that the deterioration is intentionally exaggerated to a large extent to ensure that it can be observed in the strain influence lines. The damaged regions are illustrated in Figure 1(c).

The nominal physics-based model is assumed to be a simplistic but representative model for the behavior of the bridges in the population under study in their nominal condition, obtained through the procedure described in Section 5. In addition to the physics-grounded latent variables, the nominal model also includes the parameter $ {\delta}_{\mathrm{s}} $ as an input, describing the offset of the pier relative to the center of the bridge in the longitudinal direction. Given the vertical stiffness of the abutments and support and the offset of the central pier, the nominal model returns a time-series of strains.

The response, domain and class observables are contaminated with i.i.d. samples of Gaussian white noise with standard deviations $ {\sigma}_{\mathrm{x}}={\sigma}_{\mathrm{c}}={\sigma}_{\mathrm{y}}={10}^{-4} $ . It is important to note that this case study is not intended to be a realistic representation of system identification and SHM for bridges, since it circumvents several important practical difficulties such as ensuring the consistency and alignment of data collected over long time-scales from a large number of structures. Furthermore, it is assumed for simplicity that the degradation condition of the bridges does not change significantly within the amount of time required to obtain the dataset.

5.3.2. Qualitative assessment of disentanglement

The model is trained with $ \lambda =1/1024 $ and $ {d}_{z_{\mathrm{c}}}={d}_{z_{\mathrm{y}}}=4 $ . The predictions generated by the model while traversing each of the generative factors are shown in Figure 10. It can be seen that the data-driven component of the decoder is prevented from capturing variability in the reconstructed response when varying $ {\log}_{10}\;{k}_{\mathrm{v},1} $ , $ {\log}_{10}\;{k}_{\mathrm{v},2} $ and $ {\delta}_{\mathrm{F}} $ , but is able to contribute to the components caused by the variation of the domain and class generative factors. Furthermore, the figure illustrates that the unknown confounder $ {\delta}_{\mathrm{F}} $ can be partially accounted for by the physics-based model. This is in contrast to the oscillator example (Section 5.2) where the influence of the unknown confounder could not be accounted for by the known physics.

Figure 10. Mean prediction and $ \pm 2\sigma $ uncertainty bounds for the physics-based $ {\hat{\mathbf{x}}}_{\mathrm{p}} $ and data-driven $ {\hat{\mathbf{x}}}_{\mathrm{d}} $ components, and combined prediction $ \hat{\mathbf{x}} $ while traversing the generative factors. The input response measurements are denoted as dots in the bottom row.

The previous conclusions are further supported by the traversal of the latent space, shown in Figure 11, which indicates that the domain and class subsets of the latent variables encode information that enables the auxiliary decoders to predict the domain and class labels, and the response decoder to correct the physics-based model prediction. It can be seen that the influence of the unknown confounder $ {\delta}_{\mathrm{F}} $ is partly captured as variability in the physics-based subset $ {\boldsymbol{z}}_{\mathrm{x}} $ , indicating that model form uncertainty is compensated by inferring an “effective” value of the physics-grounded latent variables. Figure 11 also suggests that the domain and class latent variables only capture variability in the corresponding generative factors, whereas the physics-grounded latent variables are always active, providing additional evidence for the claim that the adversarial objective induces a disentangled representation while prioritizing the use of the known physics. Importantly, Figures 10 and 11 illustrate that the adversarial training can feasibly constrain $ {g}_{\theta}\left({\boldsymbol{z}}_{\mathrm{c}},{\boldsymbol{z}}_{\mathrm{y}}\right) $ , such that it only contributes to the prediction when justified by additional domain and class observables.

Figure 11. Visualization of the VAE latent space during traversal of the generative factors. Each column corresponds to variation of a single generative factor, and each row shows the marginal approximate posterior distribution of a single latent variable.

5.3.3. Application to damage identification

As discussed in Section 3, the model is trained in a fully supervised manner to simultaneously reconstruct the domain and class variables from the input measurements, making it possible to handle tasks such as damage detection, where predicting the class labels $ \boldsymbol{y} $ from input measurements $ \boldsymbol{x} $ is of interest. Given a trained model, the condition labels of a similar bridge can be predicted from response measurements. The performance is evaluated in two different cases, illustrated in Figure 12, referred to as “interpolation” and “extrapolation,” respectively. For each case, the space of physics-grounded generative factors is subdivided into four quarters. In the interpolation case, the model is trained on $ {N}_{\mathrm{train}}=1024 $ samples from three quarters and evaluated on $ {N}_{\mathrm{test}}=512 $ samples from the fourth. In the extrapolation case, the model is trained on data from a single quarter and evaluated on the remaining three, using the same train and test set sizes. All other generative factors are sampled from the ground truth distributions presented in Table 3. To obtain a more comprehensive evaluation, each of the two cases is divided into four sub-cases, over which the results are averaged.

Figure 12. Samples of physics-grounded generative factors used for creating the synthetic training set (blue) and test set (orange). Two cases are constructed in order to evaluate performance in interpolation (top) and extrapolation (bottom).

The proposed disentangled physics-informed variational autoencoder (DPIVAE), using two different hyperparameter settings denoted as DPIVAE-A and DPIVAE-B, is compared with linear regression (LIN), Gaussian process regression (GPR) and a multi-layer perceptron (MLP). For DPIVAE-A the GRL is not utilized, i.e., $ \lambda =-1 $ , and separate encoders are used for each subset of the latent variables. For DPIVAE-B the GRL hyperparameter is taken as $ \lambda =1/1024 $ . The GPR is implemented with a radial basis function kernel and additive Gaussian white noise. The MLP is formulated with two hidden layers, each with a width of $ 64 $ units and a rectified linear unit (ReLU) activation function. Results in terms of interpolation and extrapolation performance of each model is quantified in terms of the R $ {}^2 $ and mean squared error (MSE), shown in Table 4. These results are intended to highlight that the performance of the different models is comparable, and that the proposed model can feasibly be used to predict the class variables in a complex high-dimensional case study. Manual hyperparameter tuning is performed for the models involved in the comparison.

Table 4. Mean and standard deviation of $ {R}^2 $ and $ \mathrm{MSE} $ for the task of predicting $ \mathbf{y} $ , averaged over $ 6 $ runs

The proposed approach using adversarial training doesn’t result in any improvement over existing approaches for the specific interpolation and extrapolation tasks in the present case study. This can be attributed to the adversarial training, which forces the encoder to weigh any information that is relevant to the prediction of $ \boldsymbol{y} $ against the potential improvement it provides towards the reconstruction of $ \boldsymbol{x} $ . It can be seen that the GPR outperforms all other models while only using a fraction of the parameters, possibly due to the smoothness of the input influence line measurements. The results also indicate that the model performs better when using a single encoder combined with adversarial training in both interpolation and extrapolation. Despite this negative result, we speculate that the disentangled representation induced by the architecture, the invariance of the class latent variables to unknown confounding influences, and the incorporation of the known physics, might be beneficial in certain tasks. Future work will aim to investigate the factors that affect the performance of the proposed approach, and the conditions under which it can provide a benefit in class prediction tasks. This analysis indicates that the proposed model performs on par with other commonly used data-driven models, but with the added benefit of ensuring the proper use of the known physics and the additional interpretability of the physics-grounded latent space.