# Model-Order-Reduction Approach for Structural Health Monitoring of Large
Deployed Structures with Localized Operational Excitations111This work was
supported by the Office of Naval Research [N00014-17-1-2077]; and the Army
Research Office [W911NF1910098]

Mohamed Aziz BHOURI222Present address: 3401 Walnut St, Wing A, Office 536,
Philadelphia, PA 19104, USA

Preprint submitted to Computers & Structures<EMAIL_ADDRESS>Department of
Mechanical Engineering, Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, MA 02139 USA

###### Abstract

We present a simulation-based classification approach for large deployed
structures with localized operational excitations. The method extends the two-
level Port-Reduced Reduced-Basis Component (PR-RBC) technique to provide
faster solution estimation to the hyperbolic partial differential equation of
time-domain elastodynamics with a moving load. Time-domain correlation
function-based features are built in order to train classifiers such as
artificial neural networks and perform damage detection. The method is tested
on a bridge example with a moving vehicle (playing the role of a digital twin)
in order to detect cracks’ existence. Such problem has $45$ parameters and
shows the merits of the two-level PR-RBC approach and of the correlation
function-based features in the context of operational excitations, other
nuisance parameters and added noise. The quality of the classification task is
enhanced by the sufficiently large synthetic training dataset and the accuracy
of the numerical solutions, reaching test classification errors below $0.1\%$
for disjoint training set of size $7\times 10^{3}$ and test set of size
$3\times 10^{3}$.

###### keywords:

Structural health monitoring , simulation-based classification , model order
reduction , domain decomposition , parametrized partial differential equations
, neural networks

††journal: Computers & Structures

##  1 Introduction

Structural Health Monitoring (SHM) has received a great deal of attention in
the civil engineering and machine learning communities. The ultimate goal of
SHM is to automatically identify damage before failure occurs for a given
system. One alternative to solve this task is the model-based approach [1], in
which an inverse problem is solved to determine all parameters of a model from
which the state of damage is inferred. The data-based approach is another
predominant method in current research and consists of two stages. First, in
an offline stage of a classification task, and before the structure is put
into service, a set of training data is collected from a rich representation
of possible healthy and unhealthy states of interest. This training dataset
can be constructed (i) by performing physical experiments [2], or (ii) by
performing synthetic experiments using a mathematical best-knowledge model [3,
4, 5]. Machine learning algorithms are then applied to the training dataset to
obtain a classifier which assigns data to the relevant diagnostic class label.
Second, in the online stage of the classification task, and during the normal
operation of the structure, the classifier is used to monitor the structure,
and to map measured data in the field to a best prediction of the
corresponding state of damage. Simulation-Based Classification (SBC) refers to
the particular choices of a data-based approach and use of synthetic
experiments.

SBC for SHM of large deployed mechanical structures, such as offshore
platforms and bridges, is still an active research field for several
challenges associated with such task. First, such systems are subject to
operational excitations [6, 7] and the most difficult loads to simulate are
those with the shortest time span. By consequence, harmonic analyses fail to
faithfully capture the response of these mechanical structures and time-domain
characterization is needed. Indeed, modal analyses are typically not very good
for local inhomogeneities since eigenfunctions do not well represent local
forces whose frequency spectra are broadly spread over certain intervals.
Second, numerical simulations of large structures are often challenging due to
the considerable amount of memory and computation capacity needed, which is
more hampering for systems with localized excitations. Such problems involve
different scales, and thus numerical methods need to be accurate enough to
well approximate the source terms, but also not be computationally
prohibitive, which is challenging due to the large geometric domain of the
structure. In this context, these forces are applied on regions which are
relatively small compared to the size of the global domain. Finally, the
response of the structure is inherently affected by probabilistic nuisance
parameters that should be considered in the training task. Nuisance parameter
refers to any parameter that affects the behavior of the structure but does
not influence its state of damage. As a consequence, numerical models should
consider the inherent uncertainty in the value of the parameters representing
material properties, geometric domains and input forces. These deviations may
be introduced in the manufacturing process or caused subsequently by
operational and environmental conditions. The list of possible variations —
present perhaps in both undamaged and damaged states — is thus relatively
long. Moreover, ambient operational excitation will certainly deviate from
ideal in terms of spatial structure and temporal signature. Therefore, it is
of interest to consider not just active systems with probabilistic nuisance
parameters in which one can provide forced and controlled input excitations,
but also passive systems which rely on ambient loading as naturally arises in
operation of the deployed system [2, 8]. The goal of this work is to develop a
SBC approach for such problems.

The SHM literature proposes a variety of output-only methods for passive
systems. In particular, Operational Modal Analysis (OMA) identifies modal
properties of structures from ambient vibration data [9, 7, 10, 11]. Numerous
modal-based damage detection techniques have been developed. These approaches
can rely on: structure’s natural frequencies; modal damping; modal shape of
displacement or its curvature; modal strain energy or modal flexibility.
Representative examples of OMA techniques include peak-picking [12], Frequency
Domain Decomposition [13, 14, 15], and time-domain Decomposition [16]. All OMA
techniques have been shown to be considerably sensitive to noise contamination
[17, 18, 19]. Moreover, the assumption of stationary random signals for
excitations (typically white noise) made by most frequency domain techniques
may not always be suited to real world applications. For instance, vehicle-
induced excitations on damaged bridges can be highly non-stationary.
Furthermore, the mode shapes cannot be identified precisely for complex
structures. Although techniques based on modal shape of displacement or its
curvature showed superior damage sensitivity compared to other OMA methods,
such techniques require many sensors to recover higher modes and their
performance heavily depends on the number of modes considered [2]. In
addition, computing curvatures from vibration data inherently introduces
additional errors due to the application of the finite difference
approximation method, which is further amplified for high-frequency noise.
Existing frequency domain-based classification approaches consider synthetic
datasets with size that goes up to the order of $10^{4}$ sample points [20]
but as stated earlier, frequency analyses fail to capture the structure’s
response to operational and localized excitations. One of the alternatives or
improvement to OMA techniques is an efficient utilization of machine learning
algorithms for damage detection. Recently, different vision-based methods for
detecting concrete cracks have been developed [21, 22, 23, 24]. Most of these
techniques rely on deep architectures of neural networks such as convolutional
networks [25, 26] and encoder-decoder networks [27]. The SBC approach
presented in this work is an alternative to efficiently incorporate machine
learning techniques into SHM for large mechanical structures under localized
operational excitations [28].

Appropriate choice of features is absolutely crucial for classification.
Features which are sensitive to the anticipated damage but relatively
insensitive to nuisance variables and measurement noise greatly simplify the
classification task and ultimately improve the robustness and hence
performance of the deployed classifier. Within those considerations,
structural damage detection methods using two-point time-domain correlation
functions of vibration response under stochastic excitation have been
developed [29, 30, 31]. The integrand of the correlation function clearly
depends on the wave speed since it counts for two displacements with a shift
in time. Therefore, it is expected to be sensitive to anticipated damage such
as loss of stiffness or crack existence, since these instances of damage
considerably affect the wave speed. Multiple damage detection methods have
been developed based on the correlation functions such as the correlation
function amplitude method (CCFA) [31] — also named correlation function
amplitude vector method (CorV) [32] — and the inner product vector method
(IPV) one [30]. The effectiveness of these methods has been shown for white
noise and steady random excitations with specific frequency spectrum by visual
inspection of the damage index. The latter is defined as the difference
between the considered vector (CCFA or IPV) of the intact and damaged
structures. Moreover, these methods were tested on relatively simple
structures: the CCFA method [31] was applied to a structure with $12$ degrees
of freedom, while the IPV one was applied in to a clamped beam with $8$
degrees of freedom [32] and to a cantilever beam with $25$ eight-node
quadrilateral shell elements [30]. These resolutions were justified by the
necessity of obtaining rapid calculation of the vibration responses. In
addition, these techniques were only tested on structures with loss of
stiffness as an instance of damage. No topological changes, such as the
existence of crack, were considered to further investigate their efficiency.
Finally, these methods rely on the choice of a reference sensor so as to
define the corresponding correlation-based vector to monitor in order to
detect the damage. Similar damage detection approaches based on correlation
functions have been built. The AMV and CZV methods [29] were developed as an
attempt to solve the reference sensor problem. However, similarly to the CCFA
and IPV methods, these techniques were only tested on relatively simple
problems (for instance a structure with only $8$ degrees of freedom), did not
account for nuisance parameters and also rely on visual inspection of the
damage index without generalization to a SBC task.

One of the main challenges associated with the application of SBC is the
construction of a sufficiently large dataset required to train the classifier.
Indeed, a prerequisite for good classifier performance is a sufficiently rich
description of undamaged and damaged states. Absent such a complete
description (a) a classifier will certainly not be able to discriminate
between different states of damage (not represented in the training dataset)
and (b) the method may not be able to identify features which can discriminate
between undamaged and damaged states. In particular, relying on a classical
Finite Element (FE) approach to construct the dataset will be prohibitively
costly. Parametric Model Order reduction (pMOR) is a mathematical and
computational field of study which aims to reduce the computational cost
associated with the estimation of the solution to a parameterized mathematical
model. Within our context, the latter consists in the time-domain
elastodynamics partial differential equation (PDE). pMOR approaches
intrinsically account for the probabilistic behavior of the parameters
governing the system. Therefore, the corresponding numerical model considers
the inherent uncertainty in the value of the parameters. Existing pMOR
techniques for time-dependent problems include Proper Orthogonal Decomposition
(POD) approaches [33, 34, 35], Greedy methods [36, 37], hybrid approaches
combining POD (in time) and Greedy procedures (in parameter space) [38], and
space-time approaches [39].

The dimensionality of the parameter domain which describes any sufficiently
rich set of possible systems — note that the parametrization characterizes
material properties, geometry, boundary conditions, and also topology —
precludes application of the classical pMOR approaches due to the well-known
curse of dimensionality. Within the context of time-domain elastodynamics of
large geometric domains with localized excitations, the recently developed
two-level PR-RBC method [40] is of great interest. The latter takes advantage
of domain decomposition techniques [41, 42, 43, 44, 45] and the frequency-time
duality [46, 47] to construct the reduced bases. Using domain decomposition
techniques define components which correspond to subdomains that form the
global domain when they are assembled. Hence, the two-level PR-RBC method
solves the curse of dimensionality issue by reducing the effective
dimensionality of the parameter spaces considered in the variational problems
to approximate. Therefore, it provides sufficiently accurate approximations
with relatively low dimensional approximation spaces, which considerably
reduces the computational cost of constructing the training dataset in the
context of SBC. Moreover, the two-level PR-RBC approach addresses not a
particular global system but a family of systems that can be built using a
predefined library of archetype components. Such global systems are formed by
connecting compatible replica of the archetype components. Hence, it offers
great flexibility in topology and geometry. As a consequence, the method is
particularly well-suited to SBC in which the dimension of the model
parametrization is perforce very large, and where topological damages need to
me modeled.

Our goal here is to develop an efficient simulation-based classification
approach for large deployed structures with localized operational excitations
by considering a numerical example of a bridge with a moving vehicle. For such
example, the parameters are all chosen based on literature review of actual
deployed bridges. In order to overcome the previously listed challenges
associated with such task, we extend the two-level PR-RBC method [40] to
account for moving loads. Moreover, using the two-point correlation functions,
we build time-domain-based features which are not only sensitive to the damage
— already guaranteed by the definition of the correlation function — but also
relatively insensitive to nuisance parameters and measurement noise. These
properties are obtained by using normalization and exploring different
strategies for the choice of the time shift; the sensors; and the direction of
the displacement considered for the correlation function. As for the
classifiers, distinct state-of-the-art learning algorithms were tested
including Artificial Neural Network (ANN) and one-vs-all Support Vector
Machine (ova-SVM). The bridge example will be considered to demonstrate (a)
the characterization of localized operational excitation in terms of nuisance
parameters, (b) the merits of time-domain-based correlation function features
in the context of localized operational excitation, other nuisance parameters
and added noise for damage localization and (c) the importance of the two-
level PR-RBC approach. The latter allows flexibility in topology and geometry,
permits the construction of dataset with sufficiently large size and
reasonable computation time, and guarantees the accuracy of the numerical
solution. These last two criteria need to be satisfied in order to obtain good
classification results. The methodology presented in this work can also be
applied to a variety of other deployed mechanical systems such as airframes,
railroad tracks, ships, …Our specific contributions can be summarized in the
following points:

  1. 1.

The PR-RBC-based classification permits the removal of reference sensor choice
dilemma since one can actually consider all the possible combinations of
correlation functions that can be built from a network of sensors and imbed
them into the feature that will be used for the classification. Classifiers
take care of learning from these features even if they are high dimensional.

  2. 2.

Instead of visual inspection of the damage index, the PR-RBC-based method
appeals to classifiers to discriminate between the different damage cases and
thus provides an automated method for damage detection without human
inspection. Moreover, most of the existing damage identification methods based
on correlation function rely on the difference between the damage index of the
undamaged and damaged state of the structure, which makes their deployment on
real word systems a hard task. Supervised classification solves this issue
since for any state of the structure, a feature that only depends on its
current state is considered, and no reference state is needed.

  3. 3.

The PR-RBC-based classification relies on the model order reduction technique
to construct a sufficiently large dataset with a reasonable computational cost
and a rich description of undamaged and damaged states by reducing the degrees
of freedom from the order of $10^{3}-10^{4}$ to the order of $10-10^{2}$. This
allows good classification performance since a classifier will be able to
discriminate between different states of damage (represented in the rich
training dataset) and to identify features which can discriminate between
undamaged and damaged states.

  4. 4.

The PR-RBC-based approach addresses full elasticity, which is carried out
sufficiently fast thanks to the use of model order reduction, and in
particular the two-level PR-RBC method. The method can be applied for 2-d and
3-d elasticity problems. For the current work, only the 2-d case is
considered. As a consequence, the system is modeled as a continuum and the
elastodynamics PDE governing its behavior is solved. This provides a more
faithful characterization of the structure than modeling it as an assembly of
rigid beams as it is the case for most of SHM systems detailed in the
literature. The rigid beam model fails to capture effects related to cracks
and other important “health” issues such as local parameters variation but was
chosen in most of the previous works due to the cheap computational cost
related to solving the corresponding mathematical model with few degrees of
freedom. Such local topological damages and other local parameters’ variation
can be considered in the elasticity model considered in this work, and the
two-level PR-RBC approach inherently permits to account for this probabilistic
nature.

  5. 5.

The PR-RBC-based classification enables including nuisance parameters which
are inherently present in actual deployed systems such as the temporal and
spatial signature of loads, damping coefficients and materials properties
which are all stochastic quantities and only determined within a certain
interval of confidence. The references cited above do not consider nuisance
parameters since they rely on visual inspection of the damage index and not on
SBC approach. The model order reduction technique and the classification task
inherently permit to account for the probabilistic nature of the nuisance
parameters.

  6. 6.

Due to its great flexibility in topology and geometry, the PR-RBC method can
accommodate for more realistic damage instances such as crack existence and is
thus particularly well-suited to SHM. In such context, the computational
framework should be able to switch in/out damage instances efficiently,
including topological changes like cracks, and should be able to treat local
variations including local excitations with support very small compared to the
size of the structure dimension which is the case for bridges in instance with
moving vehicles. The PR-RBC approach also offers flexibility with respect to
re-use of components either for variations on a given structure or even for
new but similar structures built from the same library of archetype components
considered for the two-level PR-RBC method.

  7. 7.

Finally, relying on the correlation function-based features, the bridge
example considered in this work shows that only a sparse configuration of
sensors can be required to reach sufficiently small classification errors for
damage existence identification.

This paper is organized as follows. In Section 2, we present an overview of
the recently developed two-level PR-RBC method and how to extend it in order
to apply it for elastodyanmics PDE with moving loads. In Section 3, we define
the Simulation-Based Classification task and the correlation function-based
features considered in this work. In order to demonstrate the capability and
assess the performance of the proposed technique, a numerical example is
presented in Section 4. Finally, in section 5 we summarize our key results,
discuss the limitations of the proposed approach, and carve out potential
directions for future investigation.

##  2 Two-Level PR-RBC Method

In this section we present a summary of the two-level PR-RBC method and its
extension to systems with moving loads. The two-level PR-RBC method can be
applied to any linear time-domain PDE which admits an affine representation of
the parameter. The latter can be recovered by means of empirical quadrature
procedure (EQP) [48] or empirical interpolation method EIM [49, 36] if needed.
Within our context of SBC, we restrict ourselves to the PDE of linear
elastodynamics in this work.

###  2.1 Time-Domain Equation

Let $\Omega\subset\mathbb{R}^{d}$, $d=2$ be a bounded domain, with boundary
$\partial\Omega$. The proposed method can be naturally applied to the
3D-elastodynamics case, but for purposes of presentation we limit ourselves to
the 2D case in this work. The boundary $\partial\Omega$ is assumed to be
partitioned into $\Gamma^{D}$ and $\Gamma^{N}$, such that Dirichlet boundary
conditions are imposed on $\Gamma^{D}$, while natural boundary conditions are
satisfied on $\Gamma^{N}$. Without loss of generality, the Dirichlet boundary
conditions are assumed to be homogeneous, and the non-essential boundary
conditions are assumed to be of Neumann type; more details on treating non-
homogeneous Dirichlet boundary conditions can be found in [50]. The simulation
time interval is noted $[0,T_{\rm final}]$, $T_{\rm final}>0$. Let
${X}\equiv\\{v\in[H^{1}({\Omega})]^{d}\ |\ v|_{{\Gamma}^{D}}=0\\}$ be the
Hilbert space of admissible real-valued functions; ${X}$ is imbued with inner
product $(w,v)_{{X}}\equiv\int_{{\Omega}}\nabla w\cdot\nabla v+wv\,dV$ and
induced norm $\|w\|_{{X}}\equiv\sqrt{(w,w)_{{X}}}$. The problem
parameterization is denoted ${\mu}\in{\mathcal{P}}$, where
${\mathcal{P}}\subset\mathbb{R}^{{n}_{P}}$ is a suitable compact set.

Let
$m(\cdot,\cdot;\cdot);c(\cdot,\cdot;\cdot);a(\cdot,\cdot;\cdot):\big{(}H^{1}({\Omega})\big{)}^{d}\times\big{(}H^{1}({\Omega})\big{)}^{d}\times\mathcal{P}\rightarrow\mathbb{R}$
be the bilinear forms corresponding to the mass, damping and stiffness terms
of the 2D-elastodynamics equation respectively. These bilinear forms are
defined as follows:

$m(w,v;\mu)\equiv\rho\;\int_{\Omega}w\cdot v\ dx\ ,$ (1)
$a(w,v;\mu)\equiv\frac{\nu\;E}{(1+\nu)\;(1-2\nu)}\;\int_{\Omega}\frac{\partial
w_{i}}{\partial x_{j}}\frac{\partial{v_{k}}}{\partial
x_{l}}\delta_{ik}\delta_{jl}\
dx+\frac{E}{2(1+\nu)}\;\int_{\Omega}\frac{\partial w_{i}}{\partial
x_{j}}\frac{\partial{v_{k}}}{\partial
x_{l}}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\ dx\ ,$ (2)
$c(w,v;\mu)\equiv\alpha_{\rm Ray}\;m(w,v;\mu)+\beta_{\rm Ray}\;a(w,v;\mu)\ ,$
(3)

where in equation (2) we use the convention of summation over repeated
indices, $\nu$ denotes the Poisson ratio, $\rho>0$ is the material density,
$\alpha_{\rm Ray}>0$ and $\beta_{\rm Ray}>0$ are the Rayleigh damping
coefficients, and $E>0$ is the Young’s modulus. Here
$f(\cdot,\cdot;\cdot):\big{(}H^{1}({\Omega})\big{)}^{d}\times[0,T_{\rm
final}]\times\mathcal{P}\rightarrow\mathbb{R}$ refers to linear form
corresponding to the Neumann boundary conditions imposed on $\Gamma^{N}$. In
order to realize an efficient offline-online decomposition, the bilinear and
linear forms are assumed to have an affine dependence on the parameter $\mu$.
In this work, we are interested in moving loads with time-independent shapes.
Hence, we recover the affine dependence of the linear form
$f(\cdot,\cdot,\mu_{f})$ on the parameters $\mu_{f}$ by means of EIM [49, 36],
where $\mu_{f}$ denote the parameters within $\mu$ that govern
$f(\cdot,\cdot,\mu)$. For the EIM training, we consider the location of the
load as a parameter instead of its speed. Such parameter spans the whole
domain of possible locations of the moving load. As a consequence, the speed
of the moving load, noted $V$, is simply replaced by the location of the load,
noted $l$, within the parameters $\mu_{f}$ to form the parameters
$\tilde{\mu}_{f}$ that we consider when performing the EIM.

The variational formulation of the elastodynamics equation then reads as
follows: Find $u(t\in[0,T_{\rm final}];\mu)$ such that $\forall t\in[0,T_{\rm
final}]$ ,

$m\Big{(}\frac{\partial^{2}u(t;\mu)}{\partial
t^{2}},v;\mu\Big{)}+c\Big{(}\frac{\partial u(t;\mu)}{\partial
t},v;\mu\Big{)}+a\Big{(}u(t;\mu),v;\mu\Big{)}=f(v,t;\mu)\ ,\forall v\in X\
,\forall t\in[0,T_{\rm final}]\ ,$ (4) $u(t;\mu)=0\ ,{\rm on}\ \Gamma^{D}\
,\forall t\in[0,T_{\rm final}]\ ,$ (5)

and

$u(t=0;\mu)=0\ ;\frac{\partial u}{\partial t}\Big{(}t=0;\mu\Big{)}=0\ .$ (6)

Extension to non-zero initial conditions can also be considered.

###  2.2 Frequency-Domain Equation

Let $\hat{X}\equiv\\{v\ |\ v=w+\textrm{i}y\ ,w,y\in X\\}$ be the Hilbert space
of admissible complex-valued functions; ${\hat{X}}$ is imbued with inner
product $(w,v)_{\hat{X}}\equiv\int_{{\Omega}}\nabla w\cdot\overline{\nabla
v}+w\bar{v}\,dV$ and induced norm
$\|w\|_{\hat{X}}\equiv\sqrt{(w,w)_{\hat{X}}}$, where $\overline{\ \cdot\ }$
refers to the complex conjugate operator. The complex problem parameterization
is denoted $\tilde{\mu}=(\mu,\omega)\in\tilde{\mathcal{P}}$, which corresponds
to the real-valued problem parameter concatenated with the angular frequency
$\omega$ as an additional parameter; here
$\tilde{\mathcal{P}}\in\mathbb{R}^{{n}_{P}+1}$ refers to the augmented compact
parameter set.

We assume that we have non-zero damping such that
$\int\limits_{0}^{\infty}|u(t;\mu)|dt<\infty$, and write
$u(t;\mu)=\Re\\{\hat{u}\big{[}\tilde{\mu}=(\mu,\omega)\big{]}e^{\textrm{i}\omega
t}\\}$, where $\Re$ refers to real part. It follows that
$\hat{u}(\tilde{\mu})$ satisfies the variational formulation of the Helmholtz
equation: Find $\hat{u}(\tilde{\mu})\in\hat{X}$ such that:

$\hat{a}(\hat{u}(\tilde{\mu}),v;\tilde{\mu})=\hat{f}(v;\tilde{\mu})\ ,\forall
v\in\hat{X}\ ,$ (7)

where
$\hat{a}(\cdot,\cdot;\cdot):\hat{X}\times\hat{X}\times\tilde{\mathcal{P}}\rightarrow\mathbb{C}$
is a sesquilinear form given by

$\hat{a}(\cdot,\cdot;\cdot)=-\omega^{2}\;\hat{m}(\cdot,\cdot;\cdot)+\textrm{i}\omega\;\hat{c}(\cdot,\cdot;\cdot)+\hat{k}(\cdot,\cdot;\cdot)\
.$ (8)

Here $\hat{m}(\cdot,\cdot;\cdot)$, $\hat{c}(\cdot,\cdot;\cdot)$ and
$\hat{k}(\cdot,\cdot;\cdot)$ are defined as

$\hat{m}(w,v;\tilde{\mu})=\rho\;\int_{\Omega}w\cdot\bar{v}\ dx\ ,$ (9)
$\hat{k}(w,v;\tilde{\mu})=\frac{\nu\;E}{(1+\nu)\;(1-2\nu)}\;\int_{\Omega}\frac{\partial
w_{i}}{\partial x_{j}}\frac{\partial\overline{v_{k}}}{\partial
x_{l}}\delta_{ik}\delta_{jl}\
dx+\frac{E}{2(1+\nu)}\;\int_{\Omega}\frac{\partial w_{i}}{\partial
x_{j}}\frac{\partial\overline{v_{k}}}{\partial
x_{l}}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\ dx\ ,$ (10)
$\hat{c}(w,v;\tilde{\mu})=\alpha_{\rm
Ray}\;\hat{m}(w,v;\tilde{\mu})+\beta_{\rm Ray}\;\hat{a}(w,v;\tilde{\mu})\ ,$
(11)

where in equation (10) we use the convention of summation over repeated
indices and
$\hat{f}(\cdot;\cdot):\hat{X}\times\tilde{\mathcal{P}}\rightarrow\mathbb{R}$
is a continuous anti-linear form associated to the Laplace transform of the
$f(\cdot,\cdot;\cdot)$ term.

Since $m(\cdot,\cdot;\mu)$, $c(\cdot,\cdot;\mu)$, $a(\cdot,\cdot;\mu)$ and
$f(\cdot,\cdot;\mu)$ are assumed to have an affine dependence on the parameter
$\mu$, it follows that $\hat{a}(\cdot,\cdot;\tilde{\mu})$ and
$\hat{f}(\cdot;\tilde{\mu})$ also have an affine dependence on the parameter
$\tilde{\mu}$.

###  2.3 Finite Element Approximation

####  2.3.1 Finite Element Discretization

In order to approximate equations (4) and (7), we consider a suitably refined
finite element (FE) Galerkin approximation: a triangulation
${\mathcal{T}}^{h}$ for domain $\Omega$; associated conforming FE
approximation spaces $X_{h}^{0}\subset\big{(}H^{1}(\Omega)\big{)}^{d}$ for
real-valued functions and $\hat{X}_{h}^{0}$ for complex-valued functions, both
of dimension $\mathcal{N}^{0}_{h}$. Let $X_{h}\equiv X\cap X_{h}^{0}$ and
$\hat{X}_{h}\equiv\hat{X}\cap\hat{X}_{h}^{0}$, and let
$\\{\varphi_{j}\\}_{j=1,\ldots,\mathcal{N}_{h}}$ denote the associated (real)
standard FE nodal basis; here $\mathcal{N}_{h}$ is the dimension of $X_{h}$
and $\hat{X}_{h}$.

##### Time-Domain equation

The FE approximation $u_{h}(t\in[0,T_{\rm final}];\mu)$ to $u(t\in[0,T_{\rm
final}];\mu)$ can be obtained by projecting equation (4) on $X_{h}$:
$u_{h}(t;\mu)\in X_{h}\ ,\forall t\in[0,T_{\rm final}]$ satisfies

$m\Big{(}\frac{\partial^{2}u_{h}(t;\mu)}{\partial
t^{2}},v;\mu\Big{)}+c\Big{(}\frac{\partial u_{h}(t;\mu)}{\partial
t},v;\mu\Big{)}+a\Big{(}u_{h}(t;\mu),v;\mu\Big{)}=f(v,t;\mu)\ ,\forall v\in
X_{h}\ ,\forall t\in[0,T_{\rm final}]\ ,$ (12) $u_{h}(t=0;\mu)=0\
;\frac{\partial u_{h}}{\partial t}\Big{(}t=0;\mu\Big{)}=0\ .$ (13)

##### Frequency-Domain equation

Similarly to the time-domain equation, the FE approximation
$\hat{u}_{h}(\tilde{\mu})$ to $\hat{u}(\tilde{\mu})$ can be obtained by
projecting equation (7) on $X_{h}$: $\hat{u}_{h}(\tilde{\mu})\in\hat{X}_{h}$
satisfies

$\hat{a}(\hat{u}_{h}(\tilde{\mu}),v;\tilde{\mu})=\hat{f}(v;\tilde{\mu})\
,\forall v\in\hat{X}_{h}\ .$ (14)

We now present the discrete equations. Let
$\underline{\hat{}\mathcal{A}_{h}}\in\mathbb{R}^{\mathcal{N}_{h}\times\mathcal{N}_{h}}$
and $\underline{\hat{f}_{h}}\in\mathbb{R}^{\mathcal{N}_{h}}$ be the frequency-
domain FE matrix and vector defined as

$\Big{(}\underline{\hat{}\mathcal{A}_{h}}(\tilde{\mu})\Big{)}_{qq^{\prime}}\equiv\hat{a}(\varphi_{q^{\prime}},\varphi_{q};\tilde{\mu})\
,1\leq q,q^{\prime}\leq\mathcal{N}_{h}\ ,$ (15)
$\Big{(}\underline{\hat{f}_{h}}(\tilde{\mu})\Big{)}_{q}\equiv\hat{f}(\varphi_{q};\tilde{\mu})\
,1\leq q\leq\mathcal{N}_{h}\ $ (16)

The FE basis vector for
$\hat{u}_{h}(\tilde{\mu})\in\mathbb{R}^{\mathcal{N}_{h}}$ is then given by

$\underline{\hat{}\mathcal{A}_{h}}(\tilde{\mu})\;\underline{\hat{u}_{h}}(\tilde{\mu})=\underline{\hat{f}_{h}}(\tilde{\mu})\
.$ (17)

Note that $\underline{\hat{}\mathcal{A}_{h}}(\tilde{\mu})$ is typically large
but very sparse.

####  2.3.2 Finite Element – Finite Difference Discretization

In order to solve equation (12), a finite-difference discretization scheme for
time marching with $N_{t}$ time steps is considered. Let $\Delta t=T_{\rm
final}/N_{t}$ ; $t^{j}\equiv j\;\Delta t,0\leq j\leq N_{t}$. Note that the
proposed method can be applied to any finite-difference scheme. A particular
scheme is selected in this work for sake of clarity: the unconditionally
stable Newmark-$\beta$ scheme with $\beta_{t}=\frac{1}{4}$ and
$\gamma_{t}=\frac{1}{2}$ (such that the average constant acceleration scheme,
or mid-point rule, is obtained) [51, 52].

Let $u_{h,\Delta t}^{j}(\mu)$, $0\leq j\leq N_{t}$ denote the finite element –
finite difference solution at time step $t^{j}$, and $\dot{u}_{h,\Delta
t}^{j}(\mu)$ and $\ddot{u}_{h,\Delta t}^{j}(\mu)$ the corresponding first and
second derivatives in time respectively. Since $u_{h,\Delta t}^{0}(\mu)=0$ and
$\dot{u}_{h,\Delta t}^{0}(\mu)=0$, then $\ddot{u}_{h,\Delta t}^{0}$ is
determined as the solution to

$m\Big{(}\ddot{u}_{h,\Delta t}^{0}(\mu),v;\mu\Big{)}=f(v,t=0;\mu)\ ,\forall
v\in X_{h}\ .$ (18)

The fields $\ddot{u}_{h,\Delta t}^{j}(\mu)$, $\dot{u}_{h,\Delta t}^{j}(\mu)$
and $u_{h,\Delta t}^{j}(\mu)$, $1\leq j\leq N_{t}$, are then determined as the
solutions to the following equations, respectively:

$m\Big{(}\ddot{u}_{h,\Delta t}^{j}(\mu),v;\mu\Big{)}+\Delta t\gamma_{t}\
c\Big{(}\ddot{u}_{h,\Delta t}^{j}(\mu),v;\mu\Big{)}+\Delta t^{2}\beta_{t}\
a\Big{(}\ddot{u}_{h,\Delta t}^{j}(\mu),v;\mu\Big{)}=f(v,t^{j};\mu)\\\
-c\Big{(}\dot{u}_{h,\Delta t}^{j-1}(\mu)+\Delta t(1-\gamma_{t})\
\ddot{u}_{h,\Delta t}^{j-1}(\mu),v;\mu\Big{)}-a\Big{(}u_{h,\Delta
t}^{j-1}(\mu)+\Delta t\ \dot{u}_{h,\Delta t}^{j-1}(\mu)+\Delta
t^{2}(1-\beta_{t})\ \ddot{u}_{h,\Delta t}^{j-1}(\mu),v;\mu\Big{)}\ ,\forall
v\in X_{h}\ ,$ (19) $\dot{u}_{h,\Delta t}^{j}(\mu)=\dot{u}_{h,\Delta
t}^{j-1}(\mu)+\Delta t\ \Big{[}(1-\gamma_{t})\ \ddot{u}_{h,\Delta
t}^{j-1}(\mu)+\gamma_{t}\ \ddot{u}_{h,\Delta t}^{j}(\mu)\Big{]}$ (20)
$u_{h,\Delta t}^{j}(\mu)=u_{h,\Delta t}^{j-1}(\mu)+\Delta t\ \dot{u}_{h,\Delta
t}^{j-1}(\mu)+\Delta t^{2}\ \Big{[}\Big{(}\frac{1}{2}-\beta_{t}\Big{)}\
\ddot{u}_{h,\Delta t}^{j-1}(\mu)+\beta_{t}\ \ddot{u}_{h,\Delta
t}^{j}(\mu)\Big{]}\ .$ (21)

We note that the scheme is implicit.

We next present the discrete equations in matrix form. Let
$\underline{\mathcal{M}_{h}}\in\mathbb{R}^{\mathcal{N}_{h}\times\mathcal{N}_{h}}$,
$\underline{\mathcal{C}_{h}}\in\mathbb{R}^{\mathcal{N}_{h}\times\mathcal{N}_{h}}$
and
$\underline{\mathcal{A}_{h}}\in\mathbb{R}^{\mathcal{N}_{h}\times\mathcal{N}_{h}}$
be the mass, damping and stiffness FE matrices, respectively:

$\Big{(}\underline{\mathcal{M}_{h}}(\mu)\Big{)}_{qq^{\prime}}\equiv
m(\varphi_{q^{\prime}},\varphi_{q};\mu)\ ,1\leq
q,q^{\prime}\leq\mathcal{N}_{h}\ ,$ (22)
$\Big{(}\underline{\mathcal{C}_{h}}(\mu)\Big{)}_{qq^{\prime}}\equiv
c(\varphi_{q^{\prime}},\varphi_{q};\mu)\ ,1\leq
q,q^{\prime}\leq\mathcal{N}_{h}\ ,$ (23)
$\Big{(}\underline{\mathcal{A}_{h}}(\mu)\Big{)}_{qq^{\prime}}\equiv
a(\varphi_{q^{\prime}},\varphi_{q};\mu)\ ,1\leq
q,q^{\prime}\leq\mathcal{N}_{h}\ .$ (24)

Furthermore, let $\underline{f_{h}^{j}}\in\mathbb{R}^{\mathcal{N}_{h}}$,
$0\leq j\leq N_{t}$, be the FE vectors corresponding to the linear form of the
time-domain variational formulation at time instance $t^{j}$, and let
$\underline{\mathcal{T}_{h}}\in\mathbb{R}^{\mathcal{N}_{h}\times\mathcal{N}_{h}}$
be the time marching matrix:

$\Big{(}\underline{f_{h}^{j}}(\mu)\Big{)}_{q}\equiv
f(\varphi_{q},t_{j};\tilde{\mu})\ ,1\leq q\leq\mathcal{N}_{h}\ ;$ (25)
$\underline{\mathcal{T}_{h}}(\mu)=\underline{\mathcal{M}_{h}}(\mu)+\Delta
t\;\gamma_{t}\;\underline{\mathcal{C}_{h}}(\mu)+\Delta
t^{2}\;\beta_{t}\;\underline{\mathcal{A}_{h}}(\mu)\ .$ (26)

Finally, if $\underline{u_{h,\Delta
t}^{j}}(\mu)\in\mathbb{R}^{\mathcal{N}_{h}}$, $\underline{\dot{u}_{h,\Delta
t}^{j}}(\mu)\in\mathbb{R}^{\mathcal{N}_{h}}$ and
$\underline{\ddot{u}_{h,\Delta t}^{j}}(\mu)\in\mathbb{R}^{\mathcal{N}_{h}}$,
$0\leq j\leq N_{t}$, denote the FE basis vectors for ${u_{h,\Delta
t}^{j}(\mu)}$, ${\dot{u}_{h,\Delta t}^{j}(\mu)}$ and ${\ddot{u}_{h,\Delta
t}^{j}(\mu)}$ respectively, we initialize $\underline{u_{h,\Delta
t}^{0}}(\mu)=\underline{0}$, $\underline{\dot{u}_{h,\Delta
t}^{0}}(\mu)=\underline{0}$, and $\underline{\dot{u}_{h,\Delta t}^{0}}(\mu)$
solution of

$\underline{\mathcal{M}_{h}}(\mu)\;\underline{\ddot{u}_{h,\Delta
t}^{0}}(\mu)=\underline{f_{h}^{0}}(\mu)\ ;$ (27)

we then solve for $\underline{u_{h,\Delta t}^{j}}(\mu)$,
$\underline{\dot{u}_{h,\Delta t}^{j}}(\mu)$, and
$\underline{\ddot{u}_{h,\Delta t}^{j}}(\mu)$, for $1\leq j\leq N_{t}$, from

$\underline{\mathcal{T}_{h}}(\mu)\;\underline{\ddot{u}_{h,\Delta
t}^{j}}(\mu)=\Big{[}\underline{f_{h}^{j}}(\mu)-\underline{\mathcal{C}_{h}}(\mu)\;\Big{(}\underline{\dot{u}_{h,\Delta
t}^{j-1}}(\mu)+\Delta t(1-\gamma_{t})\ \underline{\ddot{u}_{h,\Delta
t}^{j-1}}(\mu)\Big{)}-\underline{\mathcal{A}_{h}}(\mu)\;\Big{(}\underline{u_{h,\Delta
t}^{j-1}}(\mu)+\Delta t\ \underline{\dot{u}_{h,\Delta t}^{j-1}}(\mu)+\Delta
t^{2}(1-\beta_{t})\ \underline{\ddot{u}_{h,\Delta
t}^{j-1}}(\mu)\Big{)}\Big{]}\ ,$ (28) $\underline{\dot{u}_{h,\Delta
t}^{j}}(\mu)=\underline{\dot{u}_{h,\Delta t}^{j-1}}(\mu)+\Delta t\
\Big{[}(1-\gamma_{t})\ \underline{\ddot{u}_{h,\Delta
t}^{j-1}}(\mu)+\gamma_{t}\ \underline{\ddot{u}_{h,\Delta t}^{j}}(\mu)\Big{]}\
,$ (29) $\underline{u_{h,\Delta t}^{j}}(\mu)=\underline{u_{h,\Delta
t}^{j-1}}(\mu)+\Delta t\ \underline{\dot{u}_{h,\Delta t}^{j-1}}(\mu)+\Delta
t^{2}\ \Big{[}\Big{(}\frac{1}{2}-\beta_{t}\Big{)}\underline{\ddot{u}_{h,\Delta
t}^{j-1}}(\mu)+\beta_{t}\ \underline{\ddot{u}_{h,\Delta t}^{j}}(\mu)\Big{]}\
.$ (30)

This completes the FE “truth” discretization.

###  2.4 Two-Level PR-RBC Method: Overview

In this section, we provide a summary of the two-level PR-RBC method and its
extension to systems with moving loads. The two-level PR-RBC method is a
domain decomposition technique in which the global system is decomposed into
smaller components, which are referred to as instantiated components. This
decomposition creates an ensemble of parameterized instantiated components
which can be mapped to an ensemble of parameterized archetype components;
multiple instantiated components of the global system can correspond to the
same archetype component. Moreover, the domain decomposition creates an
ensemble of ports, defined as the intersection of the closures of each two
adjacent instantiated components, with the latter forming a parameterized bi-
component system. These ports can also be mapped to an ensemble of reference
ports associated with archetype bi-component systems. For simplicity, the
ports are presumed to be mutually disjoint, such that the reference port is
associated to two local ports. In the two-level PR-RBC method, the port and
component parameterizations correspond to the frequency-domain equation, and
the corresponding sesquilinear and anti-linear forms are assumed to have an
affine dependence on the parameters. By consequence, any global system built
from the library of the archetype components will be governed by a frequency-
domain PDE whose sesquilinear and anti-linear forms have an affine dependence
on the system parameter. The latter is intrinsically related to the parameters
considered for the different instantiated components forming the global
system.

####  2.4.1 Offline Stage

The two-level PR-RBC offline stage corresponds to the construction of reduced
bases to approximate the solution to the frequency-domain equation (14) within
the parameterized archetype components and over reference ports. Hence, it is
informed by the library of archetype components and reference ports, and is
independent of any subsequent (feasible) system assembly. For every reference
port, a reduced port space is built in order to approximate the solution on
the reference port joining each compatible pair of archetype components [45,
40].

The PR-RBC offline stage also includes the construction of reduced bubble
spaces for port mode liftings, and a reduced bubble space for each archetype
component with non-zero linear form [53]. These reduced bases approximate the
solution inside the archetype components domains (zero on the ports). Within
our context, the archetype component bubble spaces reduce to the bubble spaces
for inhomogeneity associated with non-zero source terms. In this work, all
reduced bubble spaces (for port mode liftings and for inhomogeneity) are
constructed by Proper Orthogonal Decomposition (POD). Finally, in the context
of an offline-online decomposition, all parameter-independent sesquilinear and
anti-linear forms needed for the PR-RBC online stage (detailed in Section
2.4.2) are computed and stored once in an offline stage.

We emphasize the important role of components. In general, the components
distribute the parameter domain: we reduce a large problem with many global
parameters to many small problems each with just a few (local) parameters.
Components also permit consideration of very large systems: even in the PR-RBC
offline stage, we are required to solve FE problems over at most pairs of
components — never the full system. Also, components provide geometry and
topology parametric variation. And finally, components permit us to more
easily justify the PR-RBC Offline investment: we may amortize the offline
effort not only over many queries for any particular global system, but over
all possible global systems in our family.

As detailed in the development of the two-level PR-RBC method [40], the
frequency-domain equation will be used in order to build the PR-RBC reduced
spaces. Since $f(t;\mu)\rightarrow u(t;\mu)$ is a linear time invariant
system, we can simply consider the frequency-domain representation of the load
to be independent of the frequency $\omega$ when forming the PR-RBC reduced
spaces in the offline stage. The exact time-dependency of the load $f(t;\mu)$
will be used within the online stage when solving the time-domain equation
after projecting it on the final reduced space. In this work, we are
interested in moving loads with time-independent shapes. By consequence, in
the offline stage, the anti-linear form $\hat{f}(\cdot;\tilde{\mu})$
introduced in (7) is actually taken as the linear form corresponding to the
actual load, but at one given location within its domain of existence. Hence,
the parameters governing $\hat{f}(\cdot;\cdot)$ are taken as the parameters
$\tilde{\mu}_{f}$ introduced in Section 2.1. The parameters $\tilde{\mu}_{f}$
correspond to the parameters $\mu_{f}$, with the load location $l$ replacing
the load speed $V$. For instance, if we consider a moving load with a Gaussian
shape, then the load amplitude and its spatial width can be taken as the
parameters forming $\mu_{f}$, along with the load speed. In the offline stage,
the load amplitude and its spatial width are kept as parameters governing
$\hat{f}(\cdot;\cdot)$, along with the load location instead of its speed such
that $\tilde{\mu}_{f}$ correspond to the load amplitude, its spatial width and
its spatial location.

####  2.4.2 Online Stage

For the online stage, we consider a global system characterized by a global
parameter $\mu$ and defined as an assembly of instantiated archetype
components. The online stage has two levels.

##### Level 1

The first-level reduction consists of evaluating the PR-RBC solutions to the
frequency-domain PDE for the global system at well-selected frequencies using
the reduced bases constructed at the offline stage. To this end, we consider a
sufficiently rich angular frequency $\omega$ set $\Xi_{\omega}$, of size
$n_{\omega}$. We also consdier a set $\Xi_{l}$ of random locations of the load
within the domain of existence of the moving load, as we did for the EIM
training to approximate $f(\cdot,\cdot;\cdot)$ (see Section 2.1) and for the
offline stage training (see Section 2.4.1). $\Xi_{l}$ has of the same size
$n_{\omega}$ as $\Xi_{\omega}$. We then define the online-train (o-t) dataset
$\Xi_{\rm
o-t}=\\{\tilde{\mu}\equiv(\mu,l,\omega);\omega\in\Xi_{\omega},l\in\Xi_{l}\\}$.
Hence $\Xi_{\rm o-t}$ can be expressed as: $\Xi_{\rm
o-t}=\\{\tilde{\mu}_{j},1\leq j\leq n_{\omega}\\}$, where
$\tilde{\mu}_{j}=(\mu,l_{j},\omega_{j})$ for $\omega_{j}\in\Xi_{\omega}$ and
$l_{j}\in\Xi_{l}$. Let $\sigma_{t}^{\rm ref}$ refer to a characteristic time
of $f(\cdot,\cdot;\cdot)$; then $\Xi_{\omega}$ is chosen as

$\Xi_{\omega}=\\{0,d\omega,\ldots,\omega_{\rm max}\\}\ ,\
d\omega=\frac{1}{\underline{c_{\omega}}\;{\sigma_{t}^{\rm ref}}}\ ,\
\omega_{\rm max}=\frac{\overline{c_{\omega}}}{{\sigma_{t}^{\rm ref}}}\ ,$ (31)

so that $n_{\omega}=\overline{c_{\omega}}\;\underline{c_{\omega}}+1$. Note
that for the PR-RBC offline stage described in Section 2.4.1, the training is
performed over a frequency set that at least contains $\Xi_{\omega}$.
Nonetheless, the PR-RBC approach is still general in the sense that we can
consider any assembly of components forming a feasible global system. For the
numerical example considered in this work, the PR-RBC offline stage is
performed over a frequency set that exactly matches $\Xi_{\omega}$.

Hence, for each $\tilde{\mu}_{j}\in\Xi_{\rm o-t}$, using the pre-computed and
stored sesquilinear and anti-linear form evaluations, the reduced bubble
functions for port mode liftings and for inhomogeneity can be computed at very
small computational cost compared to a full FE evaluation. Invoking again the
pre-computed and stored sesquilinear and anti-linear forms evaluations, the
sparse system Schur complement can be formed and solved. We opt for a Petrov-
Galerkin projection to construct the reduced system Schur complement, such
that our test space is spanned by the lifted port modes, and not by the
“harmonic” functions which form the (statically condensed) trial space. These
solutions define the PR-RBC approximations. A first crucial point related to
efficiency is the relatively low dimension of the PR-RBC space. A second
crucial point related to efficiency is the sparsity of the PR-RBC basis: the
support of a given basis function does not exceed two instantiated components
for lifted port modes, and is further restricted to just one instantiated
component for reduced bubble spaces. The level 1 reduction gives $n_{\omega}$
PR-RBC approximations, noted
$\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu}_{j})$ for
$\tilde{\mu}_{j}\in\Xi_{\rm o-t}$, to $n_{\omega}$ FE solutions
$\hat{u}_{h}(\tilde{\mu}_{j})$ for $\tilde{\mu}_{j}\in\Xi_{\rm o-t}$.

##### Level 2

The second-level reduction consists of building a “final” reduced basis from
the PR-RBC approximations
$\\{\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu}_{j})\
,\tilde{\mu}_{j}\in\Xi_{\rm o-t}\\}$ computed in Level 1 by performing a
Strong Greedy procedure to identify a reduced space $X_{\rm RB}$ of size
$\mathcal{N}^{\rm RB}$. These inexpensive PR-RBC approximations are considered
as “truth” solutions. Let $X_{\rm RB\,i}$ denote the reduced space of size
$i$, constructed prior to the $i$-th iteration of the Strong Greedy algorithm;
let $\hat{U}_{h,\mathbf{D},i}(\tilde{\mu}_{j})\in X_{\rm RB\,i}$ denote the RB
approximation to $\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu}_{j})$ obtained
using the reduced space $X_{\rm RB\,i}$. Since the PR-RBC approximations are
considered as high-fidelity solutions within the greedy algorithm, the next
snapshot we add to the reduced space $X_{\rm RB\,i}$ among the PR-RBC
solutions $\\{\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu}_{j})\
,\tilde{\mu}_{j}\in\Xi_{\rm o-t}\\}$ is selected based on the norm of the
error:
$\big{|}\big{|}\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu})-\hat{U}_{h,\mathbf{D},i}(\tilde{\mu})\big{|}\big{|}_{H^{1}(\Omega)}$,
for $\tilde{\mu}\in\Xi_{\rm o-t}$. An efficient computation of the errors
$\big{|}\big{|}\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu})-\hat{U}_{h,\mathbf{D},i}(\tilde{\mu})\big{|}\big{|}_{H^{1}(\Omega)}$,
for $\tilde{\mu}\in\Xi_{\rm o-t}$, such that only the updated quantities
depending on the new snapshot are evaluated at every iteration of the Strong
Greedy, is detailed in [50] . The Strong Greedy algorithm is stopped at the
$i$-th iteration if the relative error:

$\frac{\max\limits_{\tilde{\mu}\in\Xi_{\rm
o-t}}\big{|}\big{|}\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu})-\hat{U}_{h,\mathbf{D},i}(\tilde{\mu})\big{|}\big{|}_{H^{1}(\Omega)}}{\max\limits_{\tilde{\mu}\in\Xi_{\rm
o-t}}\big{|}\big{|}\hat{}\mathbb{U}_{h,\mathbf{D}}(\tilde{\mu})-\hat{U}_{h,\mathbf{D},1}(\tilde{\mu})\big{|}\big{|}_{H^{1}(\Omega)}}$
(32)

is below a certain threshold $\epsilon$ or if $i=\min(n_{\omega},M)$, where
$M$ is a prefixed maximum size allowed for $X_{\rm RB}$ [40].

The time-domain elastodynamics PDE is then projected and solved within the
“final” reduced space $X_{\rm RB}$ in standard fashion. The exact time
signature of the linear form is only used in the time marching performed using
the reduced space $X_{\rm RB}$. Note that the PR-RBC offline stage is
conducted prior to execution of the two-level procedure. Then, in the online
stage, for the given parameter value $\mu$, both levels of reduction are
invoked, and hence both levels must be computationally fast. In contrast, the
PR-RBC offline stage is run only once independently of the number of parameter
values considered in the evaluation of the time-domain solution, hence, we do
not give as much importance to the computational cost of the offline stage.

Similarly to the FE approximation $u_{h}(\cdot,\mu)$, the two-level reduced
basis approximation, denoted $U_{h,\mathbf{D},N}(t\in[0,T_{\rm final}];\mu)$,
is obtained by projecting equation (4) on $X_{\rm RB}$:
$U_{h,\mathbf{D},N}(t;\mu)\in X_{\rm RB}\ ,\forall t\in[0,T_{\rm final}]$ ,
such that

$m\Big{(}\frac{\partial^{2}U_{h,\mathbf{D},N}(t;\mu)}{\partial
t^{2}},v;\mu\Big{)}+c\Big{(}\frac{\partial U_{h,\mathbf{D},N}(t;\mu)}{\partial
t},v;\mu\Big{)}+a\Big{(}U_{h,\mathbf{D},N}(t;\mu),v;\mu\Big{)}=f(v,t;\mu)\
,\forall v\in X_{\rm RB}\ ,\forall t\in[0,T_{\rm final}]\ ,$ (33)
$U_{h,\mathbf{D},N}(t=0;\mu)=0\ ;\frac{\partial U_{h,\mathbf{D},N}}{\partial
t}\Big{(}t=0;\mu\Big{)}=0\ .$ (34)

We now incorporate a finite-difference scheme with the same time-
discretization notations introduced in Section 2.3.2. Let
$U_{h,\mathbf{D},N,\Delta t}^{j}(\mu)$ for $0\leq j\leq N_{t}$ denote the two-
level reduced basis-finite difference solution at time step $t^{j}$, and
$\dot{U}_{h,\mathbf{D},N,\Delta t}^{j}$ and $\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j}$ the corresponding first and second derivatives in time respectively.
We initialize $U_{h,\mathbf{D},N,\Delta t}^{0}(\mu)=0$,
$\dot{U}_{h,\mathbf{D},N,\Delta t}^{0}(\mu)=0$, and
$\ddot{U}_{h,\mathbf{D},N,\Delta t}^{0}$ solution of

$m\Big{(}\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{0}(\mu),v;\mu\Big{)}=f(v,t=0;\mu)\ ,\forall v\in X_{\rm RB}\ ;$ (35)

we then solve for $\ddot{U}_{h,\mathbf{D},N,\Delta t}^{j}(\mu)$,
$\dot{U}_{h,\mathbf{D},N,\Delta t}^{j}(\mu)$, and $U_{h,\mathbf{D},N,\Delta
t}^{j}(\mu)$, for $1\leq j\leq N_{t}$, from

$m\Big{(}\ddot{U}_{h,\mathbf{D},N,\Delta t}^{j}(\mu),v;\mu\Big{)}+\Delta
t\gamma_{t}\ c\Big{(}\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j}(\mu),v;\mu\Big{)}+\Delta t^{2}\beta_{t}\
a\Big{(}\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j}(\mu),v;\mu\Big{)}=f(v,t^{j};\mu)\\\
-c\Big{(}\dot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}(\mu)+\Delta t(1-\gamma_{t})\
\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}(\mu),v;\mu\Big{)}-a\Big{(}U_{h,\mathbf{D},N,\Delta
t}^{j-1}(\mu)+\Delta t\ \dot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}(\mu)+\Delta
t^{2}(1-\beta_{t})\ \ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}(\mu),v;\mu\Big{)}\ ,\forall v\in X_{\rm RB}\ ,$ (36)
$\dot{U}_{h,\mathbf{D},N,\Delta t}^{j}(\mu)=\dot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}(\mu)+\Delta t\ \Big{[}(1-\gamma_{t})\ \ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}(\mu)+\gamma_{t}\ \ddot{U}_{h,\mathbf{D},N,\Delta t}^{j}(\mu)\Big{]}\
,$ (37) $U_{h,\mathbf{D},N,\Delta t}^{j}(\mu)=U_{h,\mathbf{D},N,\Delta
t}^{j-1}(\mu)+\Delta t\ \dot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}(\mu)+\Delta
t^{2}\ \Big{[}\Big{(}\frac{1}{2}-\beta_{t}\Big{)}\
\ddot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}(\mu)+\beta_{t}\
\ddot{U}_{h,\mathbf{D},N,\Delta t}^{j}(\mu)\Big{]}\ ,$ (38)

respectively.

We now provide the matrix equations. Let $\underline{\mathcal{M}_{\rm
RB}}\in\mathbb{R}^{\mathcal{N}^{\rm RB}\times\mathcal{N}^{\rm RB}}$,
$\underline{\mathcal{C}_{\rm RB}}\in\mathbb{R}^{\mathcal{N}^{\rm
RB}\times\mathcal{N}^{\rm RB}}$ and $\underline{\mathcal{A}_{\rm
RB}}\in\mathbb{R}^{\mathcal{N}^{\rm RB}\times\mathcal{N}^{\rm RB}}$ be the
mass, damping and stiffness two-level reduced basis matrices, respectively,

$\underline{\mathcal{M}_{\rm RB}}(\mu)=\underline{Z_{\rm
RB}}(\mu)^{H}\;\underline{\mathcal{M}_{h}}(\mu)\;\underline{Z_{\rm RB}}(\mu),$
(39) $\underline{\mathcal{C}_{\rm RB}}(\mu)=\underline{Z_{\rm
RB}}(\mu)^{H}\;\underline{\mathcal{C}_{h}}(\mu)\;\underline{Z_{\rm RB}}(\mu),$
(40) $\underline{\mathcal{A}_{\rm RB}}(\mu)=\underline{Z_{\rm
RB}}(\mu)^{H}\;\underline{\mathcal{A}_{h}}(\mu)\;\underline{Z_{\rm RB}}(\mu),$
(41)

where $\cdot^{H}$ denotes the Hermitian transpose operator and
$\underline{Z_{\rm RB}}\in\mathbb{R}^{\mathcal{N}_{h}\times\mathcal{N}^{\rm
RB}}$ is the FE representation of $X_{\rm RB}$. The column $j$ of
$\underline{Z_{\rm RB}}$ corresponds to the coefficients of the RB basis
function $j$ from $X_{\rm RB}$ as represented by the FE nodal basis
$\\{\varphi_{i}\\}_{i=1,\ldots,\mathcal{N}_{h}}$. Similarly, let
$\underline{f_{\rm RB}^{j}}\in\mathbb{R}^{\mathcal{N}^{\rm RB}}$, $0\leq j\leq
N_{t}$, be the two-level reduced basis vectors corresponding to the linear
form of the time-domain variational formulation at time instance $t^{j}$, and
$\underline{\mathcal{T}_{\rm RB}}\in\mathbb{R}^{\mathcal{N}^{\rm
RB}\times\mathcal{N}^{\rm RB}}$ the time marching matrix,

$\underline{f_{\rm RB}^{j}}(\mu)=\underline{Z_{\rm
RB}}(\mu)^{\dagger}\;\underline{f_{h}^{j}}(\mu)$ (42)
$\underline{\mathcal{T}_{\rm RB}}(\mu)=\underline{\mathcal{M}_{\rm
RB}}(\mu)+\Delta t\;\gamma_{t}\;\underline{\mathcal{C}_{\rm RB}}(\mu)+\Delta
t^{2}\;\beta_{t}\;\underline{\mathcal{A}_{\rm RB}}(\mu)\ ,$ (43)

respectively.

Let $\underline{U_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)\in\mathbb{R}^{\mathcal{N}^{\rm RB}}$,
$\underline{\dot{U}_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)\in\mathbb{R}^{\mathcal{N}^{\rm RB}}$ and
$\underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)\in\mathbb{R}^{\mathcal{N}^{\rm RB}}$, $0\leq j\leq N_{t}$, denote
the two-level reduced basis vectors of ${U_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)$, ${\dot{U}_{h,\mathbf{D},N,\Delta t}^{j}}(\mu)$ and
${\ddot{U}_{h,\mathbf{D},N,\Delta t}^{j}}(\mu)$ respectively. We initialize
$\underline{U_{h,\mathbf{D},N,\Delta t}^{0}}(\mu)=\underline{0}$,
$\underline{\dot{U}_{h,\mathbf{D},N,\Delta t}^{0}}(\mu)=\underline{0}$, and
$\underline{\ddot{U}_{h,\mathbf{D},N,\Delta t}^{0}}(\mu)$ solution of

$\underline{\mathcal{M}_{\rm
RB}}(\mu)\;\underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{0}}(\mu)=\underline{f_{\rm RB}^{0}}(\mu)\ ;$ (44)

we then proceed to time-march as

$\underline{\mathcal{T}_{\rm
RB}}(\mu)\;\underline{\ddot{U}_{h,\mathbf{D},N,\Delta t}^{j}}(\mu)=\\\
\Big{[}\underline{f_{\rm RB}^{j}}(\mu)-\underline{\mathcal{C}_{\rm
RB}}(\mu)\;\Big{(}\underline{\dot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}}(\mu)+\Delta t(1-\gamma_{t})\
\underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}}(\mu)\Big{)}-\underline{\mathcal{A}_{\rm
RB}}(\mu)\;\Big{(}\underline{U_{h,\mathbf{D},N,\Delta t}^{j-1}}(\mu)+\Delta t\
\underline{\dot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}}(\mu)+\Delta
t^{2}(1-\beta_{t})\ \underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}}(\mu)\Big{)}\Big{]}\ ,$ (45)
$\underline{\dot{U}_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)=\underline{\dot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}}(\mu)+\Delta
t\ \Big{[}(1-\gamma_{t})\ \underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}}(\mu)+\gamma_{t}\ \underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)\Big{]}\ ,$ (46) $\underline{U_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)=\underline{U_{h,\mathbf{D},N,\Delta t}^{j-1}}(\mu)+\Delta t\
\underline{\dot{U}_{h,\mathbf{D},N,\Delta t}^{j-1}}(\mu)+\Delta t^{2}\
\Big{[}\Big{(}\frac{1}{2}-\beta_{t}\Big{)}\underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j-1}}(\mu)+\beta_{t}\ \underline{\ddot{U}_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)\Big{]}\ ,$ (47)

for $1\leq j\leq N_{t}$.

Since $X_{\rm RB}$ is a complex-valued reduced basis, the actual two-step
reduced basis approximation of the finite element - finite difference solution
$u_{h,\Delta t}^{j}(\mu)$ is given by:

$U_{h,\mathbf{D},N,\Delta
t}^{j}(\mu)=\sum\limits_{k=1}^{\mathcal{N}_{h}}\Big{(}\Re\Big{[}\underline{Z_{\rm
RB}}(\mu)\;\underline{U_{h,\mathbf{D},N,\Delta
t}^{j}}(\mu)\Big{]}\Big{)}_{k}\varphi_{k}\ ,1\leq j\leq N_{t}\ .$ (48)

If a quantity of interest $\underline{q}\in\mathbb{R}^{N_{q}}$ is considered,
then the corresponding FE output matrix
$\underline{\mathcal{Q}_{h}}\in\mathbb{R}^{N_{q}\times\mathcal{N}_{h}}$ is
multiplied by $\underline{Z_{\rm RB}}(\mu)$ once (offline) to obtain

$\underline{\mathcal{Q}_{\rm
RB}}(\mu)=\underline{\mathcal{Q}_{h}}\;\underline{Z_{\rm RB}}(\mu),$ (49)

and the two-level reduced basis approximation of $\underline{q}$ at time step
$t^{j}$ is simply given by:

$\underline{q^{j}_{\rm RB}}=\Re\Big{[}\underline{\mathcal{Q}_{\rm
RB}}(\mu)\;\underline{U_{h,\mathbf{D},N,\Delta t}^{j}}(\mu)\Big{]}\ ;1\leq
j\leq N_{t}.$ (50)

##  3 Simulation-Based Classification

In this section, a SBC approach is presented. The method takes advantage of
features that are built based on the time-domain correlation functions. The
two-level PR-RBC method detailed in section 2 is used to obtain the temporal
response of the system needed to compute the two-point correlation functions
and by consequence to inexpensively form quasi-exhaustive synthetic training
datasets of the features for the classification task. In order to obtain good
classification results, a sufficiently large dataset should be constructed,
which will be computed with reasonable computation time thanks the two-level
PR-RBC method.

First, the global parameter, defined for a given system, can be regarded as
follows:

$\mu=({\mu_{\rm damage}},{\mu_{\rm nuisance}},{\mu_{\rm
control}})\in\mathcal{P}\ ,\mathcal{P}=\mathcal{P}_{\rm
damage}\times\mathcal{P}_{\rm nuisance}\times\mathcal{P}_{\rm control}\;$ (51)

and $\mu$ follows the probability density function
$\rho_{\mu}\equiv(\rho_{{\mu_{\rm damage}}}\,\rho_{{\mu_{\rm
nuisance}}}\,\rho_{{\mu_{\rm control}}})$. ${\mu_{\rm damage}}$ refers to the
parameters that determine the damage status of the structure. This could be
the existence of a crack, or the Young modulus for a subdomain to $\Omega$ to
model the loss of stiffness as a damage instance. Note that due to its great
flexibility in topology and geometry, the two-level PR-RBC method can
naturally accommodate for crack existence and is thus particularly well-suited
to SHM. ${\mu_{\rm control}}$ refers to the parameters that the operator can
control like the location of the sensors used to monitor the structure or any
parameter relative to the operation conditions that can be chosen (e.g. load
parameters, engine parameters …) and ${\mu_{\rm nuisance}}$ englobes any other
parameter that is neither controlled by the operator nor does it define the
damage status. This includes for instance material properties that are
determined with uncertainty — such as damping coefficients for elastodynamic
problems for example — or deviations relative to the parameters used in the
model — like the excitation characteristics for instance.

Let $C^{d}$ denote the set containing the indices of the components that can
present a damage. For $c\in C^{d}$, the locations of $n_{\rm sensors}(c)$
sensors (noted as $x_{i}^{\rm sensor}(c)\ ,1\leq i\leq n_{\rm sensors}(c)$)
are specified and the PR-RBC-based approximation of the outputs for a
parameter $\mu$ is defined as:

$\underline{q^{c,j}_{\rm
RB}}(\mu)\equiv\begin{bmatrix}U_{h,\mathbf{D},N;1}(x_{1}^{\rm
sensor}(c),j\;t^{j},\mu)\\\ \vdots\\\ U_{h,\mathbf{D},N;1}(x_{n_{\rm
sensors}(c)}^{\rm sensor}(c),t^{j},\mu)\\\ U_{h,\mathbf{D},N;2}(x_{1}^{\rm
sensor}(c),t^{j},\mu)\\\ \vdots\\\ U_{h,\mathbf{D},N;2}(x_{n_{\rm
sensors}(c)}^{\rm sensor}(c),t^{j},\mu)\end{bmatrix}\ ,1\leq j\leq N_{t}\
,\Delta t=\frac{T_{\rm final}}{N_{t}}\ ,c\in C^{d}\ .$ (52)

where $U_{h,\mathbf{D},N}(\cdot,\cdot,\mu)$ refers to the two-level PR-RBC
approximation of the displacement for a parameter $\mu$, and in this case
$n_{\rm outputs}(c)=2\times n_{\rm sensors}(c)$.

###  3.1 Correlation Functions-Based Features

Appropriate choice of features is absolutely crucial for classification.
Features which are sensitive to the anticipated damage but relatively
insensitive to nuisance parameters and measurement noise greatly simplify the
classification task and ultimately improve the robustness and hence
performance of the deployed classifier. Within those considerations,
structural damage detection methods using correlation functions of vibration
response under stochastic excitation have been developed [29, 30, 31]. The
normalized two-point correlation function $C_{i,j,k,l}(\tau,\mu)$is defined
as:

$C_{i,j,k,l}(\tau,\mu,c)\equiv{\small\frac{1/T_{\rm
final}}{\max\limits_{\begin{subarray}{c}1\leq p\leq n_{\rm sensors}(c)\\\
0\leq t\leq T_{\rm final}\end{subarray}}u_{k}(x_{\rm
sensor}^{p}(c),t,\mu)\max\limits_{\begin{subarray}{c}1\leq p\leq n_{\rm
sensors}(c)\\\ 0\leq t\leq T_{\rm final}\end{subarray}}u_{l}(x_{\rm
sensor}^{p}(c),t,\mu)}}\times\int_{0}^{T_{sim}-\tau}u_{k}(x_{\rm
sensor}^{i}(c),t,\mu)u_{l}(x_{\rm sensor}^{j}(c),t+\tau,\mu)dt\\\ ,\ 1\leq
i,j\leq n_{\rm sensors}(c),1\leq k,l\leq d,0\leq\tau\leq T_{\rm final}$ (53)

where $T_{\rm final}$ denotes the simulation time and $u_{k}(x,\cdot,\cdot)$
refers to the displacement in the $k-$th direction at the point $x$ and
$x_{\rm sensor}^{i}(c)$ gives the location of the $i-$th sensor belonging to
the component $c$. Based on the choice of its integrand, the two-point
correlation function clearly depends on the wave speed since it counts for two
displacements with a shift in time. Therefore, it is expected to be sensitive
to anticipated damage such as loss of stiffness or crack existence, since
these instances of damage considerably affect the wave speed. Using
normalization and exploring different strategies for the choice of the time
shift; the sensors; and the direction of the displacement considered for the
correlation function, the goal is to build features which are not only
sensitive to the damage — already guaranteed by the definition of the
correlation function — but also relatively insensitive to nuisance variables
and measurement noise. Features based on the correlation function amplitude
(CCFA) method [31], on the inner product vector (IPV) method [30] and on the
AMV method [29] were considered. The best results were obtained using the CCFA
and the IPV-based features. Hence, in this work, a particular attention is
given to the IPV-based features for sake of clarity and conciseness. For $c\in
C^{d}$ and a global parameter $\mu$, we define the IPV-based feature
$\mathcal{F}^{\rm IPV}_{c}(\mu)$ as:

$(\mu\in\mathcal{P},c\in C^{d})\rightarrow\mathcal{F}^{\rm
IPV}_{c}(\mu)\equiv\Big{[}\big{(}C_{i,j,k,k}(\tau=0,\mu,c)\big{)}_{1\leq
i,j\leq n_{\rm sensors}(c),1\leq k\leq d}\Big{]}\ ,$ (54)

but for which no reference sensor needs to be chosen as done in the work of
[30]. Therefore, for the IPV feature, ${n_{\rm features}}(c)=d\times n_{\rm
sensors}(c)^{2}$. Features containing the two-point correlation function only
for one direction can also be considered by defining for instance the
following ${\rm IPV_{x}}$ feature, noted $\mathcal{F}_{c}^{\rm IPV_{x}}(\mu)$
for which $k=1$:

$(\mu\in\mathcal{P},c\in C^{d})\rightarrow\mathcal{F}_{c}^{\rm
IPV_{x}}(\mu)\equiv\Big{[}\big{(}C_{i,j,1,1}(\tau=0,\mu,c)\big{)}_{1\leq
i,j\leq n_{\rm sensors}(c)}]\ ,$ (55)

such that for the ${\rm IPV_{x}}$ feature, ${n_{\rm features}}(c)=n_{\rm
sensors}(c)^{2}$.

###  3.2 Classification Task

In order to perform the classification task, binary damage classes are defined
for every component in $c\in C^{d}$ by considering the mapping $D_{c}$ defined
as:

${\mu_{\rm damage}}\in\mathcal{P}_{\rm
damage}\xrightarrow{\;\,D_{c}\,\>}l\in\\{1,2\\}\ ,c\in C^{d}\ ,$ (56)

where by convention, $l=1$ corresponds to no damage case for component $c$ and
$l=2$ to the damaged case. Let $F_{c}$ denote the mapping of the outputs
$\\{\underline{q^{c,j}_{\rm RB}}(\mu)\\}_{1\leq j\leq N_{t}}$ to the feature
$\mathcal{F}_{c}(\mu)$ for a given system parameter $\mu$ as follows:

$\\{\underline{q^{c,j}_{\rm RB}}(\mu)\\}_{1\leq j\leq
N_{t}}\xrightarrow{\;\,F_{c}\;\,}\mathcal{F}_{c}(\mu)\in\mathbb{R}^{{n_{\rm
features}}(c)}\ ,c\in C^{d}\ ,$ (57)

where in our case the feature $\mathcal{F}_{c}(\mu)$ will either be
$\mathcal{F}_{c}^{\rm IPV}(\mu)$ or $\mathcal{F}_{c}^{\rm IPV_{x}}(\mu)$.

Finally, for $c\in C^{d}$, let $Q^{\rm RB}_{c}$ denote the two-level PR-RBC-
based computation process of the outputs $\\{\underline{q^{c,j}_{\rm
RB}}(\mu)\\}_{1\leq j\leq N_{t}}$ for a given system parameter $\mu$ such that
:

$\mu\xrightarrow{\;\,Q_{c}^{\rm RB}\;\,}\\{\underline{q^{c,j}_{\rm
RB}}(\mu)\in\mathbb{R}^{n_{\rm outputs}(c)}\\}_{1\leq j\leq N_{t}}\ ,c\in
C^{d}\ .$ (58)

Then, the considered classifier $C_{c}$ is applied:

$\mathcal{F}_{c}(\mu)\in\mathbb{R}^{{n_{\rm
features}}(c)}\xrightarrow{\;\,{C_{c}}\,\;}{l}\in\\{1,2\\}\ ,$ (59)

and the classifier that minimizes the expected misclassification is selected:

${C_{c}}^{*}\equiv\arg\min_{C_{c}}\;\mathbb{E}_{\mu}\Big{[}L\Big{(}D_{c}({\mu_{\rm
damage}}),{C_{c}}(\,F_{c}(Q^{\rm RB}_{c}(\mu))\,)\Big{)}\Big{]}\ ,$ (60)

where

$L(k,k^{\prime})={\scriptsize\begin{cases}0\quad k=k^{\prime}\\\ 1\quad k\neq
k^{\prime}\end{cases}}\quad k,k^{\prime}\in\mathbb{N}^{q},\ {\rm for\ some}\
q\in\mathbb{N}\ .$ (61)

As for the classifiers, Artificial Neural Network (ANN) [54] and Support
Vector Machine (SVM) [55, 56] are considered. Off-the-shelf Matlab
implementations of One-vs-all SVM and ANN are considered. More precisely,
fitcsvm Matlab function (using Sequential Minimal Optimization to solve the
dual problem) is used for binary Soft-SVM with 10-fold cross validation and
the maximum of the score (likelihood that a label comes from a particular
class) is consider as the classification criteria. Matlab train function is
called for the ANN with $1$ hidden layer containing $10$ nodes and tan-sigmoid
activation function. The output layer has a softmax function and the loss
function considered is cross-entropy. A scaled cross-conjugate gradient back-
propagation method is used with $80\%$ of the dataset for training and $20\%$
for validation. The training stops if validation error does not decrease after
6 iterations. Finally, the maximum number of epochs is fixed at $1000$.

Once a deployed structure is monitored, the classifiers would be tested on
features computed from experimental data and which can be expressed as:

$\mathcal{F}^{\rm exp}_{c}=\mathcal{F}_{c}(q_{c})+\eta+\epsilon_{b}$ (62)

where $\eta$ is the measurement noise and $\epsilon_{b}$ the bias introduced
by relying on the mathematical model to represent the deployed structure
dynamics. Hence deviations from the numerical feature in the test dataset
should be considered such that for some system parmaeter $\mu\in\mathcal{P}$,
$\mathcal{F}^{\rm exp}_{c}$ is close to $\mathcal{F}_{c}(Q^{\rm
RB}_{c}(\mu))$. Therefore, classifiers are trained on features computed using
noiseless data, and tested on features based on data perturbed with noise
$\eta$, such that for every sensor located at $x^{\rm sensor}_{i}(c)$ , $c\in
C^{d}$, $1\leq i\leq n_{\rm sensors}(c)$, the displacement
$U_{h,\mathbf{D},N;k}(x^{\rm sensor}_{i}(c),j\;\Delta t,\mu)$ for $1\leq k\leq
d$ , $1\leq j\leq N_{t}$ is replaced by:

$U_{h,\mathbf{D},N;k}(x^{\rm sensor}_{i}(c),t^{j},\mu)+\eta\ ,$

where $\eta$ is randomly sampled following the Gaussian probability density
function

$\mathcal{N}(0,\sigma\times\max\limits_{1\leq j\leq
N_{t}}\\{U_{h,\mathbf{D},N;k}(x^{\rm sensor}(c)_{i},t^{j},\mu)\\})\ ,$

where $\sigma$ is a predefined noise factor.

Then, a feature train-test dataset $\Xi^{\text{t-t}}$ of size
$n^{\textrm{t-t}}$ is considered, such that the classifier is trained on the
noiseless subset $\Xi^{\text{t-t}}_{\phi}$ of size $\phi\times
n^{\textrm{t-t}}$ and tested on the set $\Xi^{\text{t-t}}_{1-\phi}$ containing
the remaining points of $\Xi^{\text{t-t}}$ which are not picked in
$\Xi^{\text{t-t}}_{\phi}$. (Hence, $\Xi^{\text{t-t}}_{1-\phi}$ is of size
$(1-\phi)\times n^{\textrm{t-t}}$). Moreover, the data within the set
$\Xi^{\text{t-t}}_{1-\phi}$ is perturbed with noise $\eta$ using the noise
factor $\sigma$ as explained above. Finally, $n_{\rm part}$ different random
partitions for each $\Xi^{\text{t-t}}$ are considered and the
misclassification error averaged over the different random partitions is
reported.

The classification approach described above is conduct for every component in
$C^{d}$ and is summarized in the Train-Test Synthetic Supervised Learning
algorithm (Algorithm 1). In steps 10 and 11, the noise is added to the test
dataset as explained above.

Algorithm 1 T-T-Learning algorithm

1:Create t(rain)-t(est) parameter sample:
$\Xi^{\text{t-t}}\equiv\\{\mu^{\text{t-t}}_{j}\\}_{j=1,\ldots,n^{\textrm{t-t}}},\;\mu^{\text{t-t}}_{j}\text{
i.i.d.}\sim\rho_{\mu}$

2:for $p:=1{\rm\ to\ }n_{\rm part}$ do

3: Consider train-test partition of
$\Xi^{\text{t-t}}\equiv\Xi^{\text{t-t}}_{\phi}\bigcup\Xi^{\text{t-t}}_{1-\phi}$

4: for $c{\rm\ in\ }C^{d}$ do

5: Train Classifier $C_{c}$ on $\Xi^{\text{t-t}}_{\phi}$ without noise

6: for $\mu\in\Xi^{\text{t-t}}_{1-\phi}$ do

7: for $i:=1{\rm\ to\ }n_{\rm sensors}(c)$ do

8: for $k:=1{\rm\ to\ }d$ do

9: for $j:=1{\rm\ to\ }N_{t}$ do

10: $\eta\sim\mathcal{N}(0,\sigma\times\max\limits_{1\leq j\leq
N_{t}}\\{U_{h,\mathbf{D},N;k}(x^{\rm sensor}_{i}(c),t^{j},\mu)\\})$

11: $U_{h,\mathbf{D},N;k}(x^{\rm sensor}_{i}(c),t^{j},\mu)\leftarrow
U_{h,\mathbf{D},N;k}(x^{\rm sensor}_{i}(c),t^{j},\mu)+\eta$

12: end for

13: end for

14: end for

15: end for

16: Test Classifier $C_{c}$ on $\Xi^{\text{t-t}}_{1-\phi}$ perturbed by noise
$\eta$

17:
$\text{Err}_{\eta;c}(p)\leftarrow\frac{1}{(1-\phi)n^{\textrm{t-t}}}\displaystyle\sum_{\mu\in\Xi^{\text{t-t}}_{1-\phi}}L\Big{(}D_{c}(\mu),C_{c}(\,\mathcal{F}_{c}(\,Q_{c}^{\rm
RB}(\mu)\,))\Big{)}$

18: end for

19: $\text{Err}_{\eta;c}\leftarrow\frac{1}{n_{\rm
part}}\displaystyle\sum_{p=1}^{n_{\rm part}}\text{Err}_{\eta;c}(p)$

20:end for

21:return $\\{\text{Err}_{\eta;c}\\}_{c\in C^{d}}$

Afterwards, the classification task to estimate the binary state of the whole
structure is conducted, such that the structure is considered as damaged if
there is at least one damaged component. The corresponding binary damage
classes are defined by considering the mapping $D_{b}$ defined as:

${\mu_{\rm damage}}\in\mathcal{P}_{\rm
damage}\xrightarrow{\;\,D_{b}\,\>}l\in\\{1,2\\}\ ,$ (63)

such that:

$D_{b}({\mu_{\rm damage}})=\begin{cases}2\ {\rm if\ there\ is\ any}\ c\in
C^{d}\ {\rm s.t.}\ D_{c}({\mu_{\rm damage}})=2\ ,\\\ 1\ {\rm if\ for\ any}\
c\in C^{d}\ ,\ D_{c}({\mu_{\rm damage}})=1\ .\end{cases}$ (64)

No additional classification training is required for this task since the
binary structure state, noted $C_{b}(\mu)$, can be inferred using the
classifiers $C_{c}$ for $c\in C^{d}$ as follows:

$C_{b}(\mu)=\begin{cases}2\ {\rm if\ there\ is\ any}\ c\in C^{d}\ {\rm s.t.}\
{C_{c}}(\,\mathcal{F}_{c}(Q^{\rm RB}_{c}(\mu))\,)=2\ ,\\\ 1\ {\rm if\ for\
any}\ c\in C^{d}\ ,\ {C_{c}}(\,\mathcal{F}_{c}(Q^{\rm RB}_{c}(\mu))\,)=1\
.\end{cases}$ (65)

Then, the expected misclassification for the binary structure state defined
as:

$\text{Err}_{\eta;s}\equiv\mathbb{E}_{\mu}\Big{[}L\Big{(}D_{b}({\mu_{\rm
damage}}),C_{b}(\mu)\Big{)}\Big{]}\ ,$ (66)

can be estimated. Finally, the classification task to precisely determine the
state of all the components that can be damaged is considered. Thus, the goal
is to discriminate between the $2^{{\rm card}(C^{d})}$ possible values that
the vector $[D_{c}({\mu_{\rm damage}})]_{c\in C^{d}}$ can take, where ${\rm
card}(C^{d})$ refers to the size of the set $C^{d}$. Hence this classification
task tends to precisely determine the exact state of the whole structure based
on its possible local damages. The classes for such task can be defined by
considering the mapping $D$ defined as:

${\mu_{\rm damage}}\in\mathcal{P}_{\rm
damage}\xrightarrow{\;\,D\,\>}[D_{c}({\mu_{\rm damage}})]_{c\in C^{d}}\ .$
(67)

No additional classification training is required for this task since the
states of all the components $c\in C^{d}$, noted $C(\mu)$, can be inferred
using the classifiers $C_{c}$ for $c\in C^{d}$ as follows:

$C(\mu)=[{C_{c}}(\,\mathcal{F}_{c}(Q^{\rm RB}_{c}(\mu))\,)]_{c\in C^{d}}\ .$
(68)

Then, the expected misclassification for the structure state defined as:

$\text{Err}_{\eta;a}\equiv\mathbb{E}_{\mu}\Big{[}L\Big{(}D({\mu_{\rm
damage}}),C(\mu)\Big{)}\Big{]}\ ,$ (69)

can be estimated.

To obtain a classifier of sufficient quality, features should be sufficiently
discriminating between damage instances and sufficiently insensitive to small
perturbations induced by nuisance parameters. The train-test sample size
$n^{\textrm{t-t}}$ should also be sufficiently large such that damage cases
can be distinguished in the many-query context. Finally, the numerical
approximation, corresponding to the two-level PR-RBC solution
$U_{h,\mathbf{D},N}(\cdot,\cdot)$ in this work, should be sufficiently
accurate and close to the high-fidelity approximation, corresponding to the
finite element solution $u_{h}(\cdot,\cdot)$ in this work, and thus features
do not also significantly deviate from those obtained using the high-fidelity
approximation. All these criteria will be verified in the next numerical
example.

##  4 Bridge Example

In this section, the two-level PR-RBC-based SBC approach is applied to a
bridge with a moving $2$-axle vehicle problem. In the model, the vehicle-
bridge interaction is considered to be the most significant at the connection
between the bridge’s decks, since this region generally presents gaps and
metallic connectors and thus the vehicle’s passing in such region induces
vibrations that propagate through the bridge. Therefore, the vehicle-bridge
interaction is neglected at other regions than the bridge’s decks connection.

###  4.1 Problem Definition

####  4.1.1 Archetype Components and Bi-component Systems

To simulate the bridge model described above, the archetype components defined
in figure 1 are considered and the associated reference ports are shown in
figure 2. The archetype components $1$ has a rectangular geometry of dimension
$\frac{3}{2}L\times H$ and contains a homogeneous Dirichlet boundary and
homogeneous Neumann boundaries. The archetype component number $2$ has a T
shape with thickness equal to $H$ as detailed in figure 1 and also contains a
homogenous Dirichlet boundary and homogeneous Neumann boundaries. The
archetype components $3$ and $4$ have a rectangular geometry of dimension
$L\times H$ and do not contain any Dirichlet boundary. Component $3$ has only
homogeneous Neumann boundary conditions, while component $4$ will model the
connection between the decks and thus has a non-homogeneous Neumann term
corresponding to traction,

$\sigma(u)\cdot
n=\begin{bmatrix}F\;e^{-\frac{[x_{1}-(-L/2+V\;t)]^{2}}{\sigma_{x}^{2}}}\;\mathbbm{1}_{\\{x_{2}=H\\}}\;\mathbbm{1}_{\\{-d_{1}-4\;\sigma_{x}\leq-L/2+V\;t\leq
d_{2}+4\;\sigma_{x}\\}}\\\ -c_{\rm
friction}\;F\;e^{-\frac{[x_{1}-(-L/2+V\;t)]^{2}}{\sigma_{x}^{2}}}\;\mathbbm{1}_{\\{x_{2}=H\\}}\;\mathbbm{1}_{\\{-d_{1}-4\;\sigma_{x}\leq-L/2+V\;t\leq
d_{2}+4\;\sigma_{x}\\}}\end{bmatrix}\ ,\forall
x_{1}\in\Big{[}-\frac{L}{2},\frac{L}{2}\Big{]}\ ,x_{2}\in\\{0,H\\}\ ,$ (70)

applied around the midpoint of its upper boundary. Here $\sigma(u)$ denotes
the stress tensor, $F$ the load amplitude, $V$ the vehicle speed, $d_{1}$ and
$d_{2}$ the limits of the vehicle-bridge interaction, $\sigma_{x}$ the load
spatial width, $c_{\rm friction}$ the fraction coefficient, $\mathbbm{1}_{C}$
the 2D-function equal to $1$ if the condition $C$ is satisfied and $0$
otherwise, and $n$ the outer normal of the geometric domain. This particular
spatial dependence of vehicle-bridge interaction is chosen based on the
results obtained in the works of [57, 58], which show profiles close to moving
Gaussians.

Finally, the archetype component number $5$ has a rectangular geometry of
dimension $L\times H$ with non-homogeneous Neumann boundary condition similar
to the one applied for the archetype component number $4$ (70), but archetype
component $6$ also contains a crack whose dimensions are detailed in figure 1.
This archetype component $6$ will be used along with the archetype component
$4$ to model the existence or not of damage.

Figure 1: Library of archetype components (considered for construction of
elastodynamic bridges with moving vehicle) Figure 2: Library of reference
ports (archetype bi-components) constructed from the archetype components
defined in figure 1

Each of archetype components number $1$, $2$ and $3$ has $3$ parameters which
consist of the Young’s modulus and the two Rayleigh damping coefficients as
introduced in equations (2) and (3). Hence, the corresponding frequency-domain
variational problem has $4$ parameters — including the angular frequency. Each
of archetype components number $4$ and $5$ has $6$ additional parameters
defining the non-homogeneous Neumann boundary condition detailed in equation
(70): the load amplitude $F$, the load center defined by the vehicle speed
$V$, the load spatial width $\sigma_{x}$, the friction coefficient $c_{\rm
friction}$, and the two parameters $d_{1}$ and $d_{2}$ defining the interval
length around the midpoint of the upper boundary such that the vehicle-bridge
interaction is nonzero if the load center is within the spatial interval
$[-d_{1}-4\;\sigma_{x},d_{2}+4\;\sigma_{x}]$. The frequency-domain variational
problem for archetype components $4$ and $5$ has $10$ parameters in total. The
reference ports parameters and the boundary conditions considered for the
associated bi-component problems follow naturally from the parameters and the
boundary conditions defined for the archetype components. The variational
problem for archetype component-pairwise training has either $7$ parameters,
for reference ports not involving the archetype component $4$ or $5$, or $13$
parameters for the reference ports involving the archetype component $4$ or
$5$. These local parameters associated with the frequency-domain variational
problems define the parameters spaces considered in building the PR-RBC
reduced bases that are used at Level 1 reduction.

$L$ and $H$ are taken equal to $5\ {\rm m}$ and $1\ {\rm m}$ respectively. The
bridge’s material is chosen as reinforced concrete since it is one of the most
used materials in building bridges. Thus, the Rayleigh damping coefficients
$\alpha_{\rm Ray}$ and $\beta_{\rm Ray}$ follow the uniform distributions over
the intervals

$[0.566\ {\rm s^{-1}},4.311\ {\rm s^{-1}}]\ {\rm and}\ [0.009\ {\rm s},0.021\
{\rm s}]$

respectively, as established in the work of [59] for highway bridges in
concrete. The Young’s modulus also follows a uniform distributions over

$[29\ {\rm GPa},37\ {\rm GPa}]$

as determined in the work of [60] for reinforced concrete. The spatial width
of the vehicle-bridge interaction $\sigma_{x}$ and its amplitude $F$ also
follow uniform distributions over the intervals

$[2\ {\rm cm},4\ {\rm cm}]\ {\rm and}\ [10^{6}\ {\rm Pa},2\times 10^{6}\ {\rm
Pa}]$

respectively, based on the results obtained in the works of [57, 58]. Given
the work of [57] (and also using https://www.engineeringtoolbox.com), the
friction coefficient between a typical car tire and asphalt can be considered
following a uniform distribution over

$[0.5,0.7]\ .$

The distance along which the load is applied, and thus the span of the region
defining the connection between the bridge’s decks, is determined by the two
parameters $d_{1}$ and $d_{2}$. Based on the work of [61] that gives the
details of the possible connection options for prefabricated bridge elements
in reinforced concrete, $d_{1}$ and $d_{2}$ follow a uniform distribution over

$[d_{\rm min}=10\ {\rm cm},d_{\rm max}=15\ {\rm cm}]\ .$

Moreover, the geometry of the crack for component $6$ was fixed such that it
has the maximum allowable dimensions fixed by Federal Highway Administration
of the U.S. Department of Transportation [62, 63]. Finally. for reinforced
concrete, $\rho=2400\ {\rm kg.m^{-3}}$ and $\nu=0.15$ [64] (using
https://www.engineeringtoolbox.com and https://www.concrete.org/).

Since the system’s excitation is a moving vehicle, to determine the
characteristic time of the load $\sigma_{t}^{\rm ref}$ (see equation (31)) and
thus the angular frequency set $\Xi_{\omega}$, the vehicle speed needs to be
specified. Vehicles with speed following a uniform distribution between
$V_{\rm min}=15\ {\rm km.h^{-1}}$ and $V_{\rm max}=50\ {\rm km.h^{-1}}$ are
considered. We choose $\sigma_{t}^{\rm ref}=\frac{d_{\rm min}}{V_{\rm max}}$
and take

$d\omega=\frac{V_{\rm min}}{2\;d_{\rm max}}\ ,\ \omega_{\rm
max}=10\;\frac{V_{\rm max}}{2\;d_{\rm min}}\ ,$ (71)

such that $\underline{c_{\omega}}=10$, $\overline{c_{\omega}}=5$ and
$n_{\omega}=51$ (see equation (31)).

The sizes of the different reduced bases formed at Level 1 reduction are
chosen based on the decrease of the eigenvalues of the transfer eigenvalue
problem and the decrease of the POD modes for the reduced bubble space for
inhomogeneity, the reduced port space, and the reduced space for port mode
lifting. Table 1 gathers the sizes of the different reduced bases and the
computation time to run the offline stage needed by Level 1 reduction. All
simulations considered in this work were run on a 4-core laptop (with a 3.5
GHz Intel CPU and 16 GB RAM). We also provide the PR-RBC offline cost, though
this cost is amortized over the many online Level 1 - Level 2 queries.

Size of training set $\Xi_{\rm o-t}$ |  $51$
---|---
Size of port spaces |  $3$ for reference ports $1$ and $4$ $4$ for reference ports $2$ and $3$
Size of bubble spaces for port mode lifting |  $2$
Size of bubble space for inhomogeneity for archetype components $4$ and $5$ |  $1$
Computation time to run PR-RBC offline stage |  $64\;s$
Table 1: PR-RBC reduced bases sizes for elastodynamics bridge

####  4.1.2 Global System

We consider the global system consisting of $n_{\rm comp}=23$ instantiated
components presented in figure 3. The mapping of each instantiated component
of the global system to the corresponding archetype component is given in
table 2. Since each of components number $8$ and $16$ can have a crack, each
of them is either mapped to the archetype component $4$ or $5$.

Figure 3: Note-to-scale representation of the global system for the bridge: the Gaussian curves on the top boundary indicate the vehicle’s load applied on the bridge Instantiated component number in global system |  Archetype component number
---|---
$1,23$ |  $1$
$3,5,7,9,11,13,15,17,19,21$ |  $3$
$2,6,10,14,18,22$ |  $2$
$4,12,20$ |  $4$
$8,16$ |  $4$ or $5$
Table 2: Instantiated components to archetypes mapping

Note that in addition to the parameters defined for every component, two more
parameters for the global system have to be considered: the vehicle speed $V$
and the distance between the two axles noted $d_{a}$ as illustrated in figure
3. $d_{a}$ follows a Gaussian distirubtion with mean equal to
$\overline{d_{a}}=3\ {\rm m}$ and a standard deviation equal to
$\sigma_{d_{a}}=0.5\ {\rm m}$ based on generic dimensions of 2-axle vehicles,
while $V$ is sampled from a uniform probability density function define over
$[V_{\rm min},V_{\rm max}]$. The time-domain variational problem for the
global system has a total of $45$ parameters: the bridge’s two Rayleigh
damping coefficients; the Young’s modulus for each component; the vehicle-
bridge’s interaction width, the load amplitude and the friction coefficient
for each axle ($6$ parameters); the vehicle speed $V$; the distance between
the two axles $d_{a}$; the parameters $d_{1}$ and $d_{2}$ defining the spatial
region along which the load is applied ($10$ parameters); and the two
parameters defining the possible existence of crack within components $8$ and
$16$. As mentioned above, the frequency-domain variational problem has either
$4$ or $10$ parameters for the archetype components, and $7$ or $13$
parameters for the archetype bi-components defining the reference ports. Those
sizes are lower than the size of the parameter space of the global domain
(equal to $45$) and thus it shows how the two-level PR-RBC method reduces the
effective dimensionality of the parameter spaces considered in the variational
problems. This reduction of the size of the parameter spaces is even more
enhanced as larger global domains with more instantiated components are
considered.

For the time-domain problem approximation, the simulation time $T_{\rm final}$
is fixed based on the vehicle speed $V$ such that it travels across the entire
bridge during $T_{\rm final}$. The size of the reduced space $X_{\rm RB}$
constructed by Strong Greedy approach (see Paragraph Level 2 of Section 2.4.2)
is fixed by imposing $\epsilon=10^{-5}$ as a threshold for the relative error
(32). As expected for the second order Newmark-$\beta$ scheme
($\beta_{t}=\frac{1}{4}$, $\gamma_{t}=\frac{1}{2}$), the order of convergence
in time-discretization is equal to $p=2$. The number of time steps $N_{t}$,
and equivalently the step-size $\Delta t$, are fixed based on the second order
convergence of the normalized Richardson’s extrapolation-based error indicator
given by:

$\epsilon_{\Delta t}\equiv\frac{1}{\max\limits_{1\leq j\leq
N_{t}}\Big{|}\Big{|}U_{h,\mathbf{D},N,\Delta
t}^{j}(\mu)\Big{|}\Big{|}_{H^{1}(\Omega)}}\;\frac{\max\limits_{1\leq j\leq
N_{t}/2}\Big{|}\Big{|}U_{h,\mathbf{D},N,\Delta
t}^{2\;j}(\mu)-U_{h,\mathbf{D},N,2\Delta
t}^{j}(\mu)\Big{|}\Big{|}_{H^{1}(\Omega)}}{2^{p}-1}\ .$ (72)

The number of time steps $N_{t}$ is fixed such that we impose
$\epsilon_{\Delta t}\leq 10^{-4}$.

For purposes of presentation, we consider four global parameters $\mu_{{\rm
example},i}$, $i=1,\ldots,4$ corresponding to the following cases:

  1. 1.

Case 1: $\mu_{{\rm example},1}$ is such that all parameters are equal to their
average value (average of their probability density functions) and we don’t
have any crack (components $8$ and $16$ are both mapped to archetype component
$4$).

  2. 2.

Case 2: $\mu_{{\rm example},2}$ is such that all parameters are equal to their
average value and component $8$ has a crack (component $8$ is thus mapped to
the archetype component $5$, while component $16$ is mapped to the archetype
component $4$).

  3. 3.

Case 3: $\mu_{{\rm example},3}$ is such that all parameters, except $d_{a}$,
are equal to their minimum allowable values which consist in the lower bound
of the intervals used for the uniform probability density functions and
$d_{a}=\overline{d_{a}}-4\times\sigma_{d_{a}}$. Component $16$ has a crack
(component $8$ is thus mapped to the archetype component $4$, while component
$16$ is mapped to the archetype component $5$).

  4. 4.

Case 4: $\mu_{{\rm example},4}$ is such that all parameters, except $d_{a}$,
are equal to their maximum allowable values which consist in the upper bound
of the intervals used for the uniform probability density functions and
$d_{a}=\overline{d_{a}}+4\times\sigma_{d_{a}}$. We also consider 2 cracks
(component $8$ and $16$ are both mapped to the archetype component $5$).

We also consider $10$ randomly sampled global parameters $\mu_{{\rm
rand}\,i}$, $i=1,\ldots,10$ and note

$\Xi_{\rm o}\equiv\\{\mu_{{\rm example},i},i=1,\ldots,10\\}\cup\\{\mu_{{\rm
rand}\,j},j=1,\ldots,10\\}\ .$ (73)

The size of $X_{\rm RB}$ obtained by imposing $\epsilon=10^{-5}$ is on average
equal to $N=30$ for the parameters $\Xi_{\rm o}$, while imposing
$\epsilon_{\Delta t}\leq 10^{-4}$ requires $N_{t}=10^{4}$ based on the
convergence of $\epsilon_{\Delta t}$ with $N_{t}$. Table 3 gathers the
computation time to estimate $U_{h,\mathbf{D},N,\Delta t}^{j}(\mu)$ , $1\leq
j\leq N_{t}$, averaged over the parameters $\Xi_{\rm o}$. The size of the full
$\mathbb{P}_{2}$ FE approximation space is $\mathcal{N}^{h}=17532\pm 36$ and
computing one full FE simulation with $N_{t}=10^{4}$ takes $215\;s$ on average
(without performing the Richardson extrapolation). For confirmation purposes,
we consider the time-domain relative error between the two-level PR-RBC
solution and the FE approximation for a global parameter $\mu$,

$\frac{||U_{h,\mathbf{D},N,\Delta t}^{j}(\mu)-u_{h,\Delta
t}^{j}(\mu)||_{H^{1}(\Omega)}}{\max\limits_{1\leq j\leq N_{t}}||u_{h,\Delta
t}^{j}(\mu)||_{H^{1}(\Omega)}}\ ,1\leq j\leq N_{t}\ .$

Figure 4 gives the evolution of the relative error for the global parameters
$\mu_{{\rm example},i}$, $i=1,\ldots,4$ and we can verify that it is well
below $1\%$, confirming the sufficiently refined time discretization thanks to
the criterion $\epsilon_{\Delta t}<10^{-4}$, and the sufficiently rich reduced
space $X_{\rm RB}$ thanks to the strong greedy criterion imposed with
$\epsilon=10^{-5}$.

PR-RBC online stage called $n_{\omega}=51$ times |  $2.92\;s$
---|---
Strong greedy |  $0.56\;s$
Time marching |  $0.25\;s$
Total computation time to estimate $U_{h,\mathbf{D},N,\Delta t}^{j}(\mu)$ , $1\leq j\leq N_{t}$ |  $3.73\;s$
Table 3: Computation time of the two-level reduced basis method for
elastodynamics bridge Figure 4: Time-domain relative error for the
elastodynamics bridge and the global parameters $\mu_{{\rm example},i}$,
$i=1,\ldots,4$

In conclusion, the two-level reduction approach has a computation cost $58$
times lower than the FE simulation for this bridge example. In the next
section, we show the usefulness of such computational speedup in the context
of Simulation Based Classification for Structural Health Monitoring.

###  4.2 Classification Task

A synthetic train-test dataset is constructed using the two-level PR-RBC
method for the SBC task of the global system given in figure 3. To do so, the
parameters are sampled following their governing distributions introduced in
the previous subsection, and the corresponding two-level PR-RBC time-domain
solutions are computed. Moreover, for this classification task, components $8$
and $16$ can be damaged based on the existence or not of a crack, such that
$C^{d}=\\{8,16\\}$. For every component $c\in C^{d}$, $n_{\rm sensors}(c)=4$
sensors are chosen and located as shown in figure 5.

Figure 5: Note-to-scale representation of the global system for the bridge
with sensors location

Note that the potential cracks are located at $(8\times L,H)$ and $(16\times
L,H)$. The damage is defined as the existence of crack, thus ${\mu_{\rm
damage}}\equiv(\theta_{8},\theta_{16})$, where $\theta_{i}$ follows a uniform
distribution over the set $\\{1,2\\}$ and component $i$ has a crack if
$\theta_{i}=1$.

The T-T-Learning algorithm (Algorithm 1) is run for each component $c\in
C^{d}$ and for different features and classifiers. A particular focus will be
given to the effect of the size $n^{\textrm{t-t}}$ of the train-test
$\Xi^{\text{t-t}}$ on the accuracy of the classification. For the different
results given this section, $\phi$ is taken equal to $0.7$, which means that
for every train-test set $\Xi^{\text{t-t}}$ considered, $70\%$ of
$\Xi^{\text{t-t}}$ are used to train the classifier, while the remaining
$30\%$ are used to test it, and $n_{\rm part}=100$ random partitions are
considered.

Figure 6 shows the decrease of the expected misclassification for the
different components $c\in C^{d}$, the expected misclassification for the
binary structure state and the expected misclassification for the structure
state ($4$-class classification task in this example), with the increase of
the train-test dataset size $n^{\textrm{t-t}}$ and obtained using the ${\rm
IPV}$ and ${\rm IPV_{x}}$ features and ANN. The figure contains the results
obtained for noiseless test dataset and for test dataset perturbed with noise
with a factor equal to $\sigma=0.02$. This value was chosen based on error
estimates carried out for bridge displacement computed from measured
acceleration records [65]. Computer vision-based techniques for bridge
displacement measurements generally present higher noise levels [66]. As
$n^{\textrm{t-t}}$ gets close to $10^{4}$, the ${\rm IPV_{x}}$ feature gives
misclassification errors of orders around or below $10^{-4}$. Such observation
is still valid even for the structure state classification which is a
$4$-class classification task and even with a test dataset perturbed by noise.
Figure 6 also contains the misclassification error corresponding to only one
misclassified point out of the whole test dataset, which is plotted in black
dashed line. The misclassification errors obtained correspond to less than one
misclassified point as $n^{\textrm{t-t}}$ gets close to $10^{4}$ for all cases
considered except for the ${\rm IPV}$ feature with $\sigma=0.02$.

As mentioned above, the computation time for one simulation using the two-
level PR-RBC approach is equal to $3.73s$, while computing the features takes
$0.05s$ (for the $2$ different noise levels considered). All simulations
considered in this work were run on a 4-core laptop (with a 3.5 GHz Intel CPU
and 16 GB RAM). To obtain the results presented above, a train-test dataset
$\Xi^{\text{t-t}}$ of size $n^{\textrm{t-t}}=10^{4}$ needs to be constructed.
Conducting such a task using the two-level PR-RBC approach has a total
computation time of $10.52$ hours (taking into account the computational cost
of the offline stage, which is run only once), as opposed to an estimated
$24.88$ days using the full FE approximation.

Figure 6: Misclassification errors variation with $n^{\textrm{t-t}}$ using
${\rm IPV}$, ${\rm IPV_{x}}$ features and ANN with $\phi=0.7$

Sensors located further away from the potential cracks (centered at $(8\times
L,H)$ and $(16\times L,H)$) are considered and the new sensors’ positions are
shown in figure 7 where the sensors belonging to the lower boundary of the
decks have the same location as in the previous example, while the sensors
located at the upper boundary are now $0.5\ {\rm m}$ away from the potential
cracks instead of the $0.2\ {\rm m}$ chosen before.

Figure 7: Note-to-scale representation of the global system for the bridge
with different sensors location

Figure 8 shows the expected misclassifications obtained with the ${\rm
IPV_{x}}$ features and ANN with $\sigma=0$ and $\sigma=0.02$ for the two
different choices of sensors’ locations. As expected, the misclassification
error is higher for the new sensors’ positions since they are located further
away from the potential cracks. Nonetheless, for noiseless test dataset, the
change in expected misclassification is not as significant as for test dataset
with noise. Indeed, for $\sigma=0$ and as $n^{\textrm{t-t}}$ gets close to
$10^{4}$, the expected misclassification obtained for the new sensors’
positions is on the same order as the one obtained for the sensors located
closer to the cracks’ locations. However, for test dataset altered with noise,
the expected misclassification obtained for for sensors located further from
the cracks is higher; and even for $n^{\textrm{t-t}}$ close to $10^{4}$, the
expected misclassifications are one order higher than those obtained for
sensors located closer to the cracks. For instance, the expected
misclassification obtained with the closest sensors located at $0.5\ {\rm m}$
from the cracks is not less than one misclassified point out of the test
dataset as $n^{\textrm{t-t}}$ gets close to $10^{4}$, while it is the case
when the closest sensors are located at $0.2\ {\rm m}$ from the cracks.
Nevertheless, the expected misclassifications always presents an almost
monotonous decrease with the increase of the train-test dataset size
$n^{\textrm{t-t}}$. Hence, a longer distance between the sensors and the
cracks along with a reasonable noise level requires larger train-test dataset,
which emphasizes more the importance of pMOR and more specifically of the two-
level PR-RBC approach in the context of large structure with localized
excitations such that the sufficiently large dataset can be built in a
reasonable computation time. The same conclusions can be drawn regarding the
effect of the sensors’ locations when considering the ${\rm IPV}$ features as
shown in figure 9.

Figure 8: Misclassification errors variation with $n^{\textrm{t-t}}$ for
different sensors locations using ${\rm IPV_{x}}$ features, ANN and with
$\phi=0.7$: $d_{\rm min}^{c-s}$ refers to the minimum distance between the
potential crack location and the closest sensor Figure 9: Misclassification
errors variation with $n^{\textrm{t-t}}$ for different sensors locations using
${\rm IPV}$ features, ANN and with $\phi=0.7$: $d_{\rm min}^{c-s}$ refers to
the minimum distance between the potential crack location and the closest
sensor

Finally, in order to show its importance on the classification results, the
accuracy of the two-level PR-RBC-based solution is intentionally reduced by,
for instance, reducing the size $\mathcal{N}^{\rm RB}$ of the final reduced
space $X_{\rm RB}$ obtained at Level 2 reduction. The size needed to have a
relative error in time of order of magnitude equal to $10^{-3}$ is around $30$
as shown in figure 4. If $\mathcal{N}^{\rm RB}$ is limited to $15$, the
relative error in time has an order of magnitude equal to $10^{-2}$. Figure 10
gives the misclassification results obtained with the ${\rm IPV_{x}}$ features
and ANN, for $\sigma=0$ and $\sigma=0.02$, using the new sensors’ locations
(i.e. the closest sensor to any crack is $0.5\ m$ away from it) and with the
size of the final reduced space $\mathcal{N}^{\rm RB}$ equal to $30$ and $15$.
For any noise level considered, the misclassification accuracy suffers
significantly when the two-level PR-RBC-based solution was intentionally made
not accurate enough. The same conclusions can be drawn for the ${\rm IPV}$
features and ANN as shown in figure 11.

Figure 10: Effect of numerical solution accuracy on misclassification errors
using ${\rm IPV_{x}}$ features, ANN and with $\phi=0.7$ Figure 11: Effect of
numerical solution accuracy on misclassification errors using ${\rm IPV}$
features, ANN, with $\phi=0.7$

The accuracy of the numerical solution is crucial for any feature, any noise
level and any classifier considered. Indeed, figure 12 gives the
misclassification results obtained with the ${\rm IPV_{x}}$ features and SVM,
for $\sigma=0$ and $\sigma=0.02$, using the new sensors’ locations (i.e. the
closest sensor to the crack is $0.5\ {\rm m}$ away from it) and with
$\mathcal{N}^{\rm RB}$ equal to $30$ and $15$. Again, for any noise level
considered, the misclassification accuracy suffers significantly when the two-
level PR-RBC-based solution was made not accurate enough on purpose.
Therefore, sufficiently accurate numerical solutions that can be computed
within reasonable computation time to be able to build sufficiently large
train-test synthetic dataset are required to have satisfactory classification
results in the context of SHM for large structures with localized excitation,
and the two-level PR-RBC method fulfills all these requirements.

Figure 12: Effect of numerical solution accuracy on misclassification errors
using ${\rm IPV_{x}}$ and SVM, with $\phi=0.7$

##  5 Conclusions

###  5.1 Summary

This work presents the development of a SBC approach for SHM of large deployed
mechanical structures. The method takes advantage of the two-level PR-RBC
method as a pMOR approach to approximate the hyperbolic PDE of time-domain
elastodynamics for large domains with localized source terms and/or local
parameters variation. The usefulness of such approach was demonstrated by
considering a system with localized operational excitations and nuisance
parameters. For such systems, frequency analyses fail to capture the structure
response to such forces and hence the two-level PR-RBC approach developed for
PDE of time-domain elastodynamics is of great usefulness.

For the simulation task, features based on the two-point correlation function
were built and state-of-the-art machine learning algorithms were considered to
perform a damage detection on the structure. The features considered are based
on damage indices defined in previous works but which were never used within a
context of SBC. Thus, the damage indices were already shown to be sensitive to
the damage — thanks to the definition of the two-point correlation function —
but the generalization of their former definition showed that the features
used for the classification task in this work are relatively insensitive to
nuisance parameters and measurement noise. Within the generalization, the
reference sensor choice dilemma was solved since all possible combinations of
the two-point correlation functions are considered and imbedded in the feature
that will be sued for the classification. Moreover, instead of visual
inspection of the damage indices, classifiers are exploited to automatically
discriminate between the different damage cases. The two-level PR-RBC approach
was used to construct a sufficiently large dataset with rich description of
undamaged and damaged states. The system was modeled as a continuum governed
by the elastodynamics PDE, providing a more faithful characterization than
modeling the mechanical structure as an assembly of rigid beams. The two-level
PR-RBC approach lets us considerably reduce the degrees of freedom of the
numerical system and thus the computational cost related to the construction
of the dataset. In addition, the model order reduction technique and the
classification task inherently permit to account for the probabilistic nature
of the nuisance parameters. Finally, thanks to the great flexibility in
topology, the two-level PR-RBC approach can accommodate for more realistic
damage instances such as crack existence and is thus particularly well-suited
to SHM. Test classification errors below $0.1\%$ were reached for disjoint
training set of size $7\times 10^{3}$ and test set of size $3\times 10^{3}$,
thus showing the strong potential of the proposed approach in view of the
application to real-life systems.

For the SBC task, a bridge example was considered in which the goal is to
detect the existence of two potential cracks and such problem shows:

  * 1.

the characterization of localized operational excitation in terms of nuisance
parameters,

  * 2.

the merits of time-domain-based correlation function features in the context
of localized operational excitation, other nuisance variables and added noise,
and

  * 3.

the importance of the two-level PR-RBC approach which considerably reduces the
computational burden associated with the construction of synthetic training
dataset. Thus, it allows the generation of such sufficiently large dataset
with reasonable computation time, while guaranteeing the accuracy of the
numerical solution. These two conditions need to be satisfied in order to
obtain good classification results as shown by the significant decrease of the
misclassification errors with the increase of the size of the training
dataset.

###  5.2 Future Work

As a future work, it would be interesting to develop computational techniques
for automatic, or at least semi-automatic, feature identification which would
include sensor placement and two-point correlation functions selection [67].
These methods promise not only more effective classifiers but also more
effective deployment. For instance, they can reduce the number of sensors
needed. To this end, recent advances in Robust Optimization [68, 69] and Mixed
Integer Optimization [70] can be useful.

It would also be interesting to consider experimental data to use as a test
for the classifiers that are built using synthetic dataset. Impact hammers
represent an interesting option for SHM actuation in particular as regards
installation and operation. Experimental enrichment of the training procedure
could also be considered if some experimental data might be available at the
offline stage of the classification task.

Finally, in many applications it may not be feasible to provide an exhaustive
experimental description of all possible states of damage. In other
applications, a detailed characterization of damage may not in fact be
interesting — but rather only a determination of “undamaged” or “damaged” (or,
in the case of quality control, “different”) is required. From the learning or
statistical perspective, this two-way classification corresponds to the
problem of anomaly, or outlier, detection [71]. For the examples considered in
this work, the classification results obtained for such binary task were
presented. Within these two contexts, investigating unsupervised learning
techniques would be of great interest.

Acknowledgement:  This work was supported by the ONR Grant [N00014-17-1-2077]
and by the ARO Grant [W911NF1910098]. We would like to thank Professor Anthony
T. Patera, Dr. Tommaso Taddei and Professor Masayuki Yano for the helpful
comments and software they provided us with.

## References

## References

  * [1] S. Khatir, M. Wahab, Fast simulations for solving fracture mechanics inverse problems using pod-rbf xiga and jaya algorithm, Engineering Fracture Mechanics 205 (2019) 285–300.
  * [2] C. R. Farrar, K. Worden, Structural Health Monitoring: a Machine Learning Perspective, 1st Edition, John Wiley & Sons, Ltd., Chichester, West Sussex, UK, 2013.
  * [3] J. E. Hurtado, An examination of methods for approximating implicit limit state functions from the viewpoint of statistical learning theory, Structural Safety 26 (3) (2004) 271–293.
  * [4] J. E. Hurtado, D. A. Alvarez, Classification Approach for Reliability Analysis with Stochastic Finite-Element Modeling, Journal of Structural Engineering 129 (8) (2003) 1141–1149.
  * [5] M. Lecerf, D. Allaire, K. Willcox, Methodology for Dynamic Data-Driven Online Flight Capability Estimation, AIAA Journal 53 (10) (2015) 3073–3087.
  * [6] B. Peeters, J. Maeck, G. DeRoeck, Vibration-based damage detection in civil engineering: excitation sources and temperature effects, Smart Materials and Structures 10 (3) (2001) 518–527.
  * [7] J. J. Moughty, J. R. Casas, A State of the Art Review of Modal-Based Damage Detection in Bridges: Development, Challenges, and Solutions, Applied Sciences 7 (5) (2017) 510.
  * [8] A. Deraemaeker, E. Reynders, G. De Roeck, J. Kullaa, Vibration-based structural health monitoring using output-only measurements under changing environment, Mechanical Systems and Signal Processing 22 (1) (2008) 34–56.
  * [9] S. Au, F.-L. Zhang, Y.-C. Ni, Bayesian operational modal analysis: Theory, computation, practice, Computers & Structures 126 (2013) 3–14.
  * [10] L. Zhang, R. Brincker, P. Andersen, An overview of operational modal analysis: Major development and issues, in: Proceedings of the International Operational Modal Analysis Conference, Copenhagen, Denmark, 2005, pp. 26–27.
  * [11] G.-R. Gillich, H. Furdui, M. Wahab, Z.-I. Korka, A robust damage detection method based on multi-modal analysis in variable temperature conditions, Mechanical Systems and Signal Processing 115 (2019) 361–379.
  * [12] N. M. M. Maia, J. M. M. Silva, Theoretical and experimental modal analysis, Research Studies Press, Boston, MA, 2011.
  * [13] R. Brincker, L. Zhang, P. Andersen, Modal identification from ambient responses using frequency domain decomposition, in: IMAC 18 : Proceedings of the International Modal Analysis Conference (IMAC), San Antonio, TX, 2000, pp. 625–630.
  * [14] X. Jiang, H. Adeli, Pseudospectra, music, and dynamic wavelet neural network for damage detection of highrise buildings, International Journal for Numerical Methods in Engineering 71 (5) (2007) 606–629.
  * [15] J. Amezquita-Sanchez, H. H. Adeli, Synchrosqueezed wavelet transform-fractality model for locating, detecting, and quantifying damage in smart highrise building structures, Smart Materials and Structures 24 (6).
  * [16] V. M. Karbhari, F. Ansari, Structural Health Monitoring of Civil Infrastructure Systems, 1st Edition, Elsevier, Amsterdam, Netherlands, 2009.
  * [17] P. Cruz, R. Salgado, Performance of Vibration-Based Damage Detection Methods in Bridges, Computer-Aided Civil and Infrastructure Engineering 24 (2008) 62–79.
  * [18] I. Talebinejad, C. Fischer, F. Ansari, Numerical Evaluation of Vibration-Based Methods for Damage Assessment of Cable-Stayed Bridges, Computer-Aided Civil and Infrastructure Engineering 26 (2011) 239–251.
  * [19] W. Fan, P. Qiao, Vibration-based Damage Identification Methods: A Review and Comparative Study, Structural Health Monitoring 10 (2011) 83–111.
  * [20] T. Taddei, J. D. Penn, M. Yano, A. T. Patera, Simulation-Based Classification; a Model-Order-Reduction Approach for Structural Health Monitoring, Archives of Computational Methods in Engineering 25 (1) (2018) 23–45.
  * [21] Q. Mei, M. Gül, A Cost Effective Solution for Road Crack Inspection using Cameras and Deep Neural Networks, arXiv:1907.06014v2 [cs.CV] (2019).
  * [22] X. Zhang, D. Rajan, B. Story, Concrete crack detection using context-aware deep semantic segmentation network, Computer-Aided Civil and Infrastructure Engineering 34 (11) (2019) 951–971.
  * [23] F. Nayyeri, L. Hou, J. Zhou, H. Guan, Foreground-background separation technique for crack detection, Computer-Aided Civil and Infrastructure Engineering 34 (6) (2019) 457–470.
  * [24] F. Ni, J. Zhang, Z. Chen, Zernike-moment measurement of thin-crack width in images enabled by dual-scale deep learning, Computer-Aided Civil and Infrastructure Engineering 34 (5) (2019) 367–384.
  * [25] Y. Cha, W. Choi, O. Büyüköztürk, Deep Learning-Based Crack Damage Detection Using Convolutional Neural Networks, Computer-Aided Civil and Infrastructure Engineering 32 (5) (2018) 361–378.
  * [26] S. Li, X. Zhao, G. Zhou, Automatic pixel-level multiple damage detection of concrete structure using fully convolutional network, Computer-Aided Civil and Infrastructure Engineering 34 (7) (2019) 616–634.
  * [27] S. Bang, S. Park, H. Kim, H. Kim, Encoder-decoder network for pixel-level road crack detection in black-box images, Computer-Aided Civil and Infrastructure Engineering 34 (8) (2019) 713–727.
  * [28] M. Rafiei, H. Adeli, A novel machine learning-based algorithm to detect damage in high-rise building structures, Smart Materials and Structures 26 (18).
  * [29] M. Y. Zhang, R. Schmidt, B. Markert, Structural damage detection methods-based on the correlation functions, in: Proceedings of the 9th International Conference on Structural Dynamics, EURODYN, Porto, Portugal, 2014.
  * [30] Z. Yang, L. Wang, H. Wang, Y. Ding, X. Dang, Damage Detection in Composite Structures Using Vibration Response under Stochastic Excitation, Journal of Sound and Vibration 325 (4) (2009) 14–16.
  * [31] L.-S. Huo, X. Li, Y.-B. Yan, H.-N. Li, Damage Detection of Structures for Ambient Loading-Based on Cross-Correlation Function Amplitude and SVM, Shock and Vibration 2016 (2016) 1–12.
  * [32] Z. Yang, Z. Yu, H. Sun, On the Cross-Correlation Function Amplitude Vector and its Application to Structural Damage Detection, Mechanical Systems and Signal Processing 21 (2007) 2918–2932.
  * [33] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problem, Numerische Mathematik 90 (2001) 117–148.
  * [34] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics, SIAM Journal on Numerical Analysis 40 (2) (2003) 492–515.
  * [35] M. Rathinam, L. R. Petzold, A new look at proper orthogonal decomposition, SIAM Journal on Numerical Analysis 41 (5) (2003) 1893–1925.
  * [36] M. A. Grepl, Y. Maday, N. C. Nguyen, A. T. Patera, Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations, ESAIM Mathematical Modelling and Numerical Analysis 41 (3) (2007) 575–605.
  * [37] M. A. Grepl, A. T. Patera, A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations, ESAIM Mathematical Modelling and Numerical Analysis 39 (1) (2005) 157–181.
  * [38] M. A. Grepl, Certified Reduced Basis Methods for Nonaffine Linear Time-Varying and Nonlinear Parabolic Partial Differential Equations, Mathematical Models and Methods in Applied Sciences 22 (3), paper 1150015.
  * [39] Y. Masayuki, An optimization framework for adaptive higher-order discretizations of partial differential equations on anisotropic simplex meshes, Ph.D. thesis, Massachusetts Institute of Technology, phD Thesis, Published (2012).
  * [40] M. A. Bhouri, A. T. Patera, A two-level parameterized model-order reduction approach for time-domain elastodynamics, arXiv:2002.11084v2 [math.NA, CS.NA] (2020).
  * [41] D. Huynh, D. Knezevic, A. Patera, A static condensation reduced basis element method: Approximation and a posteriori error estimation, ESAIM Mathematical Modelling and Numerical Analysis 47 (1) (2013) 213–251.
  * [42] D. Huynh, D. Knezevic, A. Patera, A static condensation reduced basis element method: Complex problems, Computer Methods in Applied Mechanics and Engineering 259 (2013) 197–216.
  * [43] J. L. Eftang, A. Patera, Port Reduction in Component-Based Static Condensation for Parametrized Problems: Approximation and a Posteriori Eerror Estimation, International Journal for Numerical Methods in Engineering 96 (5) (2013) 269–302.
  * [44] K. Smetana, A new certification framework for the port reduced static condensation reduced basis element method, Computer Methods in Applied Mechanics and Engineering 283 (2015) 352–383.
  * [45] K. Smetana, A. T. Patera, Optimal local approximation spaces for component-based static condensation procedures, SIAM Journal on Scientific Computing 38 (5) (2016) A3318–A3356.
  * [46] A. C. Antoulas, C. A. Beattie, S. Gugercin, Interpolatory Model Reduction of Large-Scale Dynamical Systems, in: J. Mohammadpour, K. Grigoriadis (Eds.), Efficient Modeling and Control of Large-Scale Systems, Springer, Boston, MA, 2010, pp. 3–58.
  * [47] C. Beattie, S. Gugercin, Model Reduction by Rational Interpolation, arXiv:1409.2140v1 [math.NA] (2014).
  * [48] A. T. Patera, M. Yano, An lp empirical quadrature procedure for parametrized functions, Comptes Rendus Mathematique 355.
  * [49] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An empirical interpolation method: Application to efficient reduced-basis discretization of partial differential equations, C. R. Acad. Sci. Paris Series I 339 (9) (2004) 667–672.
  * [50] M. A. Bhouri, A two-step port-reduced reduced-basis component method for time domain elastodynamic pde with application to structural health monitoring, Ph.D. thesis, Massachusetts Institute of Technology, phD Thesis, Published (2020).
  * [51] F. Chiba, T. Kako, Stability and error analyses by energy estimate for Newmark’s method, Tech. Rep. NIFS-PROC–40, National Inst. for Fusion Science, Nagoya, Japan (1999).
  * [52] R. Nickell, On the Stability of Approximation Operators in Problems of Structural Dynamics, International Journal of Solids and Structures 7 (3) (1971) 301–319.
  * [53] J. L. Eftang, A. Patera, A port-reduced static condensation reduced basis element method for large component-synthesized structures: approximation and a posteriori error estimation, Advanced Modeling and Simulation in Engineering Sciences 1 (1) (2014) 3.
  * [54] C. M. Bishop, Neural Networks for Pattern Recognition, Oxford University Press, New York, 1995.
  * [55] N. Cristianini, J. Shawe-Taylor, An introduction to support vector machines and other kernel-based learning methods, Cambridge university press, Cambridge, England, 2000.
  * [56] C. Cortes, V. N. Vapnik, Support-Vector Networks, Machine Learning 20 (3) (1995) 273–297.
  * [57] M. Yu, G. Wu, L. Kong, Y. Tang, Tire-Pavement Friction Characteristics with Elastic Properties of Asphalt Pavements, Applied Sciences 7 (2017) 1123.
  * [58] P. Yap, Truck tire types and road contact pressures, in: Proceedings of the 2nd International Symposium on Heavy Vehicle Weights and Dimensions, The Roads and Transport Association of Canada, Canada, 1989.
  * [59] M. Q. Feng, S. C. Lee, Determining the effective system damping of highway bridges, Tech. Rep. CA-UCI-2009-001, California Department of Transportation, Sacramento, CA (2009).
  * [60] M. Musiał, J. Grosel, Determining the young’s modulus of concrete by measuring the eigenfrequencies of concrete and reinforced concrete beams, Construction and Building Materials, Elsevier 121 (2016) 44–52.
  * [61] M. P. Culmo, Connection Details for Prefabricated Bridge Elements and Systems, Tech. Rep. FHWA-IF-09-010, Federal Highway Administration, Washington, DC (2009).
  * [62] F. Barth, J. R. Frosch, M. Abou-Zeid, H. J. Allen, J. P. Barlow, M. E. Brander, K. Carlson, D. Darwin, F. H. Fouad, D. W. Fowler, G. T. Halvorsen, W. Hansen, M. N. Hassoun, H. Haynes, P. Hedli, T. C. Liu, E. G. Nawy, R. W. Poston, R. J. Rhoads, A. Scanlon, E. K. Schrader, W. Suaris, Z. A. Zielinski, Control of Cracking of Concrete Structures, Tech. Rep. ACI 224R-01, American Concrete Institute, Farmington Hills, MI (2001).
  * [63] S. S. G. Balakumaran, R. E. Weyers, M. C. Brown, Linear Cracking in Bridge Decks, Tech. Rep. FHWA/VTRC 18-R13, Virginia Transportation Research Council, Charlottesville, VA (2018).
  * [64] A. Logan, W. Choi, A. Mirmiran, S. Rizkalla, P. Zia, Short-Term Mechanical Properties of High-Strength Concrete, ACI Materials Journal 106 (2009) 413–418.
  * [65] M. Gindy, H. H. Nassif, J. Velde, Bridge Displacement Estimates from Measured Acceleration Records, Transportation Research Record 2028 (1) (2007) 136–145.
  * [66] B. W. Jo, Y. S. Lee, J. H. Jo, K. R. M. A., Computer Vision-based Bridge Displacement Measurements using Rotation-Invariant Image Processing Technique, Sustainability 10 (6) (2018) 1–16.
  * [67] H. Qarib, H. Adeli, A new adaptive algorithm for automated feature extraction in exponentially damped signals, Smart Materials and Structures 24 (12).
  * [68] A. Ben-Tal, L. E. I. Ghaoui, A. Nemirovski, Robust Optimization, Princeton University Press, Princeton, 2009.
  * [69] D. Bertsimas, D. B. Brown, C. Caramanis, Theory and applications of robust optimization, SIAM Review 53 (3) (2011) 464–501.
  * [70] D. Bertsimas, R. Weismantel, Optimization Over Integers, Vol. 13, Dynamic Ideas, Charlestown, 2011.
  * [71] V. Chandola, A. Banerjee, V. Kumar, Anomaly detection: A survey, Computer-Aided Civil and Infrastructure Engineering 41 (3) (2009) 1–72.