# An Adaptive Multivariate Functional EWMA Control Chart

Christian Capezza  Department of Industrial Engineering, University of Naples
Federico II, Piazzale Tecchio 80, 80125, Napoli, Italy  Giovanna Capizzi
Department of Statistical Sciences, University of Padua, Via Cesare Battisti
241, 35121, Padova, Italy  Fabio Centofanti  Department of Industrial
Engineering, University of Naples Federico II, Piazzale Tecchio 80, 80125,
Napoli, Italy  Antonio Lepore  Biagio Palumbo  Department of Industrial
Engineering, University of Naples Federico II, Piazzale Tecchio 80, 80125,
Napoli, Italy

###### Abstract

In many modern industrial scenarios, the measurements of the quality
characteristics of interest are often required to be represented as functional
data or profiles. This motivates the growing interest in extending traditional
univariate statistical process monitoring (SPM) schemes to the functional data
setting. This article proposes a new SPM scheme, which is referred to as
adaptive multivariate functional EWMA (AMFEWMA), to extend the well-known
exponentially weighted moving average (EWMA) control chart from the univariate
scalar to the multivariate functional setting. The favorable performance of
the AMFEWMA control chart over existing methods is assessed via an extensive
Monte Carlo simulation. Its practical applicability is demonstrated through a
case study in the monitoring of the quality of a resistance spot welding
process in the automotive industry through the online observations of dynamic
resistance curves, which are associated with multiple spot welds on the same
car body and recognized as the full technological signature of the process.

Keywords: Functional Data Analysis, Profile Monitoring, Statistical Process
Control

##  1 Introduction

The main aim of statistical process monitoring (SPM) is to detect special
causes of variation possibly acting on a process, which is then said to be out
of control (OC). Otherwise, it is said to be in control (IC). As is known, SPM
is currently implemented in two phases. The first one (Phase I) is concerned
with identifying a clean data set to be assumed as representative of the IC
state of the process, named Phase I sample, while the second one (Phase II)
aims at the prospective monitoring of new observations (Qiu, , 2014), often
referred to as Phase II observations or data set.

In modern industrial SPM, the capability of new data acquisition systems is
reshaping the size and the variety of signals and measurements available,
which are often functions of time or space and can be best characterized by
functional variables or profiles (Ramsay and Silverman, , 2005; Kokoszka and
Reimherr, , 2017). The simplest approach for the monitoring of one or multiple
quality characteristics in this form is the application of classical SPM
techniques for multivariate data (Montgomery, , 2019) on scalar features
extracted from them, although this simple approach risks being problem-
specific, arbitrary, and flattening useful information.

To avoid that, profile monitoring (Noorossana et al., , 2011) aims to directly
monitor quality characteristics in the form of one or multiple functional
variables. Some recent examples are those of Jin and Shi,  (1999); Colosimo
and Pacella,  (2010); Zou et al.,  (2012); Chou et al.,  (2014); Grasso et
al.,  (2014); Paynabar et al.,  (2016); Grasso et al.,  (2016, 2017);
Menafoglio et al.,  (2018); Maleki et al.,  (2018); Wang et al.,  (2018); Ren
et al.,  (2019); Jones et al.,  (2021); Centofanti et al.,  (2021); Capezza et
al., 2022a ; Capezza et al., 2022b . These contributions are based on
functional extensions of the Shewart control chart and allow the process state
to be assessed only on the information acquired at the current point in time,
without considering previous data points, which may contain decisive
information about process condition trends. For this reason, the Shewart
control chart may be particularly ineffective in detecting small persistent
shifts.

The exponentially weighted moving average (EWMA) control chart is instead
known to overcome this limitation by evaluating the current state of the
process through a statistic calculated as a weighted average of the
observation at the current point in time and the statistic itself at the
previous point. The weighted average is determined by a single weighting
parameter, usually denoted with $\lambda$ (Qiu, , 2014), such that
$0\leq\lambda\leq 1$. Although the design, implementation, and properties of
the EWMA control chart have been widely discussed for both univariate and
multivariate quality characteristics (Crowder, , 1989; Gan, , 1993; Serel, ,
2009; Chan and Zhang, , 2000; Lucas and Saccucci, , 1990; Capizzi and
Masarotto, , 2003; Knoth, , 2007; Lowry et al., , 1992; Jones et al., , 2001),
few extensions have been proposed in the profile monitoring literature in this
regard. It is worth mentioning Fassò et al.,  (2016), who developed a
functional EWMA control chart for both univariate and multivariate functional
data. In particular, they applied a multivariate EWMA control chart to the
basis coefficients resulting from a chosen finite-dimensional representation
of the functional data in the reference sample. Ren et al.,  (2019) move
forward with the integration of the EWMA scheme with the multichannel
functional principal component analysis, developed by Paynabar et al.,  (2016)
for the Phase II monitoring of multichannel profiles.

The major challenge for any EWMA control chart is related to the choice of the
weighting parameter $\lambda$, which should minimize the average run length
(ARL) defined as the average number of observations required to signal an OC
state with respect to a specified shift (Lucas and Saccucci, , 1990). The ARL
corresponding to no shift is usually referred to as IC ARL. It is known that
the smaller the value of $\lambda$, the smaller the shift with respect to
which the ARL is the minimum for a given IC ARL. However, in most real cases,
there is no prior information available on the magnitude of the shift that may
actually occur in the OC scenario.

Therefore, methods are needed for adaptively selecting $\lambda$ such that the
resulting EWMA control chart performs adequately in a wide range of OC
conditions. In the univariate scalar setting, this issue is addressed by the
adaptive EWMA (AEWMA) control chart proposed by Capizzi and Masarotto,  (2003)
that smoothly combines EWMA and Shewhart control charting schemes, by
adaptively selecting the weighting parameter through a suitable function of
the magnitude of the difference between the current observation and the
monitoring statistic observed up to the previous point in time. This idea has
been extended to the multivariate setting (Mahmoud and Zahran, , 2010; Haq and
Khoo, , 2019), but has never been addressed in the profile monitoring
literature.

In this work, we propose a new method for the monitoring of a process through
observations of a multivariate functional quality characteristic, which is
named adaptive multivariate functional EWMA (AMFEWMA) control chart and
extends the underlying idea of Capizzi and Masarotto,  (2003) to the
multivariate profile monitoring setting. As a case study, we consider the
monitoring of a resistance spot welding (RSW) process in the automotive body-
in-white manufacturing (Zhang and Senkara, , 2011). RSW is commonly used to
join overlapping steel galvanized sheets by means of two copper electrodes to
guarantee the structural integrity and solidity of welded assemblies in each
vehicle (Martín et al., , 2014). In this context, the dynamic resistance curve
(DRC) is considered as highly representative of the physical and metallurgical
development of a spot weld (Capezza et al., , 2021), and therefore, of the
quality of the joint produced. In particular, there is a consensus among
experts regarding the correlation between the latter and the degradation
status or wear level of the electrodes (Manladan et al., , 2017). A control
chart capable of accommodating varying degrees of mean shift, contingent upon
different wear levels, without prior knowledge about the current wear level,
can prove immensely advantageous.

The paper is structured as follows. Section 2 introduces the AMFEWMA control
chart. Section 3 reports a Monte Carlo simulation study of the performance of
the AMFEWMA control chart in identifying process mean shifts of a multivariate
functional quality characteristic with respect to other non-adaptive SPM
approaches that have already been presented in the literature. The practical
applicability of the proposed method is demonstrated in Section 4 through a
case study in the SPM of multivariate DRCs from a RSW process in automotive
body-in-white manufacturing. Section 5 concludes the paper. Supplementary
Materials that provide details on data generation in the simulation study are
available online. All computations and plots have been obtained using the R
programming language (R Core Team, , 2023).

##  2 The AMFEWMA control chart

The elements of the AMFEWMA control chart are developed in the following
sections. Section 2.1 shows how smooth functional data are obtained from the
original discrete observations for each curve; Section 2.2 illustrates the
multivariate functional principal component analysis (MFPCA); Section 2.3
introduces the AMFEWMA scheme, while its design and implementation details are
referred to Section 2.4.

###  2.1 Data Smoothing

For each profile, measurements $y_{i}$, $i=1,\dots,m$, of the quality
characteristic are supposed to be collected in a discrete fashion with the
observation points $t_{1},\dots,t_{m}\in\mathcal{T}\subset\mathbb{R}$, into a
realization $X$ of a random function with values in $L^{2}(\mathcal{T})$,
i.e., the Hilbert space of square integrable functions defined on the compact
set $\mathcal{T}\subset\mathbb{R}$, with inner product between $a,b\in
L^{2}(\mathcal{T})$ given by $\langle
a,b\rangle=\int_{\mathcal{T}}a(t)b(t)dt$. Hence, methods are required to
convert discrete raw data $\\{(t_{i},y_{i})\\}_{i=1,\dots,m}$ into functional
data $X$. If the discrete data are assumed without any measurement error,
functional data can be theoretically drawn up by merely connecting the whole
set of points $\\{(t_{i},y_{i})\\}$. However, this does not represent the
ordinary situation. When measurement error is present, each discrete
observation is expressed as

$y_{i}=X\left(t_{i}\right)+\varepsilon_{i},\quad i=1,\dots,m,$ (1)

for $i=1,\dots,m$, where $\varepsilon_{i}$ are independent identically
distributed random errors with zero mean. The functional variable $X$ is
intrinsically infinite dimensional, that is infinite features would be needed
to completely specify it at each possible argument $t\in\mathcal{T}$.
Therefore, from Equation (1), data smoothing techniques are commonly used to
recover $X$ by discarding exogenous perturbation due to error terms
$\varepsilon_{i}$. A common approach consists of representing $X$ through a
linear combination of $K$ known basis functions
$\bm{\phi}=\left(\phi_{1},\dots,\phi_{K}\right)^{T}$. That is,

$X\left(t\right)=\sum_{k=1}^{K}c_{k}\phi_{k}\left(t\right)=\bm{c}^{T}\bm{\phi}\left(t\right)\quad
t\in\mathcal{T},$ (2)

where $\bm{c}=\left(c_{1},\dots,c_{K}\right)^{T}$ is the vector of basis
coefficients. The value of $K$ is not crucial unless it is sufficiently large
to capture the local behavior of the functional data (Cardot et al., , 2003).
Then, the problem of recovering $X$ from the discrete raw data reduces to
finding the estimate $\hat{\bm{c}}$ of the unknown coefficient vector
$\bm{c}$, which is obtained by minimizing the following penalized sum of
squares error

$\hat{\bm{c}}=\operatorname*{argmin}_{\bm{c}\in\mathbb{R}^{K}}\sum_{i=1}^{m}\left(y_{i}-\bm{c}^{T}\bm{\phi}\left(t_{i}\right)\right)^{2}+\lambda_{s}\bm{c}^{T}\bm{R}\bm{c},$
(3)

where $\lambda_{s}>0$ is a smoothing parameter and $\bm{R}$ is a matrix whose
entries are
$\langle\phi^{\left(d\right)}_{i}\left(t\right),\phi^{\left(d\right)}_{j}\rangle$,
with $\phi^{\left(d\right)}$ denoting the $d$-th derivative of $\phi$. The
smoothing parameter $\lambda_{s}$ is chosen by minimizing the generalized
cross-validation (GCV) criterion, which is a well-known method to achieve a
trade-off between variance and bias. This criterion takes into account the
degrees of freedom of the estimated curve that vary according to
$\lambda_{s}$. The reader may want to refer to Ramsay and Silverman,  (2005)
for further details. The matrix $\bm{R}$, which controls the penalty on the
right-hand side of Equation (3), is usually computed by setting $d=2$, i.e.,
by penalizing the function roughness.

The functional datum $\hat{X}(t)$ so estimated is

$\hat{X}\left(t\right)=\hat{\bm{c}}^{T}\bm{\phi}\left(t\right)\quad
t\in\mathcal{T}.$ (4)

The most common choice for $\bm{\phi}$ in case of non-periodic functional data
is the B-spline basis system, owing good computational properties and great
flexibility (Ramsay and Silverman, , 2005). More in general, splines are known
to be optimal in the sense of being the smoothest functions interpolating the
data (Green and Silverman, , 1993). Spline functions divide the functional
domain into subintervals, by means of breakpoints. Over any subinterval, the
spline is a polynomial of order $q$, with $q-1$ non-zero derivatives and
matching proper derivative constraints between adjacent polynomials (De Boor,
, 1978).

###  2.2 Multivariate Functional Principal Component Analysis

Let us consider a random vector $\bm{X}=(X_{1},\dots,X_{p})^{T}$ with
realizations in $\mathbb{H}$, the Hilbert space of $p$-dimensional vectors of
functions in $L^{2}(\mathcal{T})$. Given two elements
${\bm{a}},{\bm{b}}\in\mathbb{H}$, where $\bm{a}=(a_{1},\dots,a_{p})^{T}$ and
$\bm{b}=(b_{1},\dots,b_{p})^{T}$, the inner product of $\mathbb{H}$ is defined
as $\langle{\bm{a}},{\bm{b}}\rangle_{\mathbb{H}}=\sum_{j=1}^{p}\langle
a_{j},b_{j}\rangle$. We assume that ${\bm{X}}$ has mean
$\bm{\mu}=\left(\mu_{1},\dots,\mu_{p}\right)^{T}$, where
$\mu_{j}(t)=\operatorname{E}(X_{j}(t))$, $j=1,\dots,p$, $t\in\mathcal{T}$ and
covariance $\bm{C}=\\{C_{jk}\\}_{1\leq j,k\leq p}$,
$C_{jk}(s,t)=\operatorname{Cov}(X_{j}(s),X_{k}(t))$, $s,t\in\mathcal{T}$. From
the multivariate Karhunen-Loève’s Theorem (Happ and Greven, , 2018) it follows
that

$\bm{X}(t)-\bm{\mu}(t)=\sum_{l=1}^{\infty}\xi_{l}\bm{\psi}_{l}(t),\quad
t\in\mathcal{T},$

where $\xi_{l}=\langle\bm{\psi}_{l},\bm{Z}\rangle_{\mathbb{H}}$ are random
variables, say principal components scores or simply scores, such that
$\operatorname{E}\left(\xi_{l}\right)=0$ and
$\operatorname{E}\left(\xi_{l}\xi_{m}\right)=\rho_{l}\delta_{lm}$, with
$\delta_{lm}$ denoting the Kronecker delta. The elements of the orthonormal
set $\\{\bm{\psi}_{l}\\}_{l\geq 1}$,
$\bm{\psi}_{l}=\left(\psi_{l1},\dots,\psi_{lp}\right)^{T}$, with
$\langle\bm{\psi}_{l},\bm{\psi}_{m}\rangle_{\mathbb{H}}=\delta_{lm}$, are
referred to as principal components, and are the eigenfunctions of the
covariance $\bm{C}$ of $\bm{Z}$ corresponding to the eigenvalues
$\rho_{1}\geq\rho_{2}\geq\dots\geq 0$.

MFPCA aims to estimate $\bm{\psi}_{l}$ and $\rho_{l}$, $l\geq 1$, starting
from $N$ independent realizations, $\bm{X}_{1},\dots,\bm{X}_{N}$, of $\bm{X}$,
where $\bm{X}_{i}=(X_{i1},\dots,X_{ip})^{T}$, which can be obtained through
the data smoothing technique described in Section 2.1 and denoted by
$\hat{\bm{X}}_{1},\dots,\hat{\bm{X}}_{N}$. That is, each $X_{ij}$,
$j=1,\dots,p$, is estimated as
$\hat{X}_{ij}(t)=\hat{\bm{c}}_{ij}^{T}\bm{\phi}(t)$, where
$\hat{\bm{c}}_{ij}=(\hat{c}_{ij1},\dots,\hat{c}_{ijK})^{T}$. Following the
basis function expansion approach of Ramsay and Silverman,  (2005), each
component $\psi_{lj}$ of $\bm{\psi}_{l}$ is represented as the following
linear combination of the $K$ basis functions $\phi_{1},\dots,\phi_{K}$ used
to obtain $\hat{\bm{X}}_{1},\dots,\hat{\bm{X}}_{N}$

$\psi_{lj}(t)=\sum_{k=1}^{K}b_{ljk}\phi_{k}(t),\quad t\in\mathcal{T},\quad
j=1,\dots,p,\quad l=1,2,\dots,$ (5)

where $\bm{b}_{lj}=\left(b_{lj1},\dots,b_{ljK}\right)^{T}$ are the
eigenfunction coefficient vectors. Under these assumptions, MFPCA (Ramsay and
Silverman, , 2005; Chiou et al., , 2014) consists of performing standard
multivariate principal component analysis of the random vectors
$\bm{W}^{1/2}(\hat{\bm{c}}_{i}-\overline{\bm{c}})$, where
$\hat{\bm{c}}_{i}=\left(\hat{\bm{c}}_{i1}^{T},\dots,\hat{\bm{c}}_{ip}^{T}\right)^{T}$,
$\overline{\bm{c}}=\sum_{i=1}^{N}\hat{\bm{c}}_{i}/N$ and $\bm{W}$ is a block-
diagonal matrix with diagonal blocks $\bm{W}_{j}$, $j=1,\dots,p$, whose
entries are $w_{k_{1}k_{2}}=\langle\phi_{k_{1}},\phi_{k_{2}}\rangle$,
$k_{1},k_{2}=1,\dots,K$. Specifically, the estimated eigenvalues
$\hat{\rho}_{l}$ of $\bm{C}$ are the eigenvalues of the sample covariance
matrix $\bm{S}$ obtained from
$\bm{W}^{1/2}(\hat{\bm{c}}_{i}-\overline{\bm{c}})$, whereas
$\hat{\bm{\psi}}_{l}=\left(\hat{\psi}_{l1},\dots,\hat{\psi}_{lp}\right)^{T}$
is obtained through Equation (5), where $\bm{b}_{lj}=\bm{W}^{-1/2}\bm{u}_{lj}$
and $\bm{u}_{l}=\left(\bm{u}_{l1}^{T},\dots,\bm{u}_{lp}^{T}\right)^{T}$ is the
$l$-th eigenvector of $\bm{S}$.

In practice, it is assumed that, as the eigenvalues decrease toward zero, the
leading eigenfunctions tend to reflect the most important features of
$\bm{X}$. That is, $\hat{\bm{X}}_{i}$ are approximated by
$\hat{\bm{X}}_{i}^{L}$ through the following truncated principal component
decomposition

$\hat{\bm{X}}_{i}^{L}(t)=\hat{\bm{\mu}}(t)+\hat{\bm{D}}(t)\sum_{l=1}^{L}\hat{\xi}_{il}\hat{\bm{\psi}}_{l}(t)\quad
t\in\mathcal{T},$ (6)

where $\hat{\bm{D}}$ is a $p\times p$ diagonal matrix whose diagonal entries
are $\hat{v}_{j}^{1/2}$ and
$\hat{\xi}_{il}=\langle\hat{\bm{\psi}}_{l},\hat{\bm{X}}_{i}\rangle_{\mathbb{H}}$.
In the profile monitoring literature, the parameter $L$ is generally chosen
such that the retained principal components explain at least a given
percentage of the total variability, which is usually in the range 70-90$\%$
(Paynabar et al., , 2016; Ren et al., , 2019; Centofanti et al., , 2022;
Capezza et al., 2022b, ), even though more sophisticated methods can be used
as well (Jolliffe and Cadima, , 2016; Capezza et al., , 2020).

###  2.3 The monitoring scheme

Let us consider a sequence of independent realizations in $\mathbb{H}$ of the
random vector $\bm{X}_{n}=(X_{n1},\dots,X_{np})^{T}$, $n=1,2,\dots$. Suppose
that (a) before an unknown observation, say $n_{0}$, the process is in control
and the mean is equal to a target value $\bm{\mu}_{0}$, (b) after $n_{0}$, the
mean becomes $\bm{\mu}_{1}$ = $\bm{\mu}_{0}$ \+ $\bm{\delta}$, where
$\bm{\delta}$ is an unknown shift, and (c) we want to detect the occurrence of
this shift as soon as possible. Without loss of generality, it is assumed
that, when the process is in control, the data are mean centered, i.e.,
$\bm{\mu}_{0}$ is a vector of zero functions. As in Fassò et al.,  (2016), an
EWMA control chart for multivariate functional data can be built based on the
statistic $\tilde{\bm{Y}}_{n}=(\tilde{Y}_{n1},\dots,\tilde{Y}_{np})^{T}$,
defined as

$\tilde{\bm{Y}}_{n}(t)=({\bm{I}}-{\bm{\Lambda}})\tilde{\bm{Y}}_{n-1}(t)+{\bm{\Lambda}}{\bm{X}}_{n}(t)=\tilde{\bm{Y}}_{n-1}(t)+{\bm{\Lambda}}\tilde{\bm{E}}_{n}(t),\quad
t\in\mathcal{T},n=1,2,\dots,$ (7)

where
$\tilde{\bm{E}}_{n}(t)=(\tilde{E}_{n1}(t),\dots,\tilde{E}_{np}(t))^{T}={\bm{X}}_{n}(t)-\tilde{\bm{Y}}_{n-1}(t)$,
${\bm{I}}$ is the $p\times p$ identity matrix and
${\bm{\Lambda}}=\text{diag}(\lambda_{1},\lambda_{2},...,\lambda_{p})$, with
$\lambda_{j}\in(0,1]$, $j=1,\dots,p$, the weighting parameters. For $n=1$,
$\tilde{\bm{Y}}_{0}$ represents the starting value, often set equal to the
target value $\bm{\mu}_{0}$.

By elaborating on Equation (7), we propose the AMFEWMA statistic
${\bm{Y}}_{n}=(Y_{n1},\dots,Y_{np})^{T}$, defined as follows

${\bm{Y}}_{n}(t)=({\bm{I}}-{\bm{\Lambda}}_{n}(t)){\bm{Y}}_{n-1}(t)+{\bm{\Lambda}}_{n}(t){\bm{X}}_{n}(t)={\bm{Y}}_{n-1}(t)+{\bm{\Lambda}}_{n}(t)\bm{E}_{n}(t),\quad
t\in\mathcal{T},n=1,2,\dots.$ (8)

Similarly to $\tilde{\bm{E}}_{n}(t)$,
${\bm{E}}_{n}(t)=(E_{n1}(t),\dots,E_{np}(t))^{T}={\bm{X}}_{n}(t)-{\bm{Y}}_{n-1}(t)$;
while, differently from $\bm{\Lambda}$,
${\bm{\Lambda}}_{n}(t)=\text{diag}(w(E_{n1}(t)),\dots,w(E_{np}(t)))$, where
$w(E_{nj}(t))=\eta(E_{nj}(t))/E_{nj}(t)$ and $\eta(E_{nj}(t))$ is a score
function evaluated at the error $E_{nj}(t)$, $j=1,\dots,p$. In this way, the
AMFEWMA statistic ${\bm{Y}}_{n}$ results in a weighted average of the current
observation ${\bm{X}}_{n}$ and the charting statistic value at the previous
time point ${\bm{Y}}_{n-1}$, as in a conventional EWMA chart, but with weights
$w(E_{nj}(t))$ changing over time. That is, the score functions $\eta$ to be
used in Equation (8) are proposed by extending the score functions for a
univariate scalar quality characteristic, proposed by Capizzi and Masarotto,
(2003), to the multivariate functional data setting. A score function
$\eta(\cdot)$ must be strictly increasing, odd ($\eta(x)=\eta(-x)$), and such
that $\eta(E_{nj}(t))$ is close to $\lambda E_{nj}(t)$ (resp. to $E_{nj}(t)$)
when $|E_{nj}(t)|$ is small (resp. is large). Based on these properties, we
propose to use as the score function either

$\eta_{1}(E_{nj}(t))=\begin{cases}E_{nj}(t)+(1-\lambda)C_{j}(t)&\text{if
}E_{nj}(t)<-C_{j}(t)\\\ \lambda E_{nj}(t)&\text{if }|E_{nj}(t)|\leq
C_{j}(t)\\\ E_{nj}(t)-(1-\lambda)C_{j}(t)&\text{if
}E_{nj}(t)>C_{j}(t)\end{cases},\quad t\in\mathcal{T},$ (9)

or

$\eta_{2}(E_{nj}(t))=\begin{cases}E_{nj}(t)[1-(1-\lambda)(1-(E_{nj}(t)/C_{j}(t))^{2})^{2}]&\text{if
}|E_{nj}(t)|\leq C_{j}(t)\\\ E_{nj}(t)&\text{otherwise}\end{cases},\quad
t\in\mathcal{T},$ (10)

where $\lambda\in(0,1]$ and $C_{j}(t)=k\sigma_{j}(t)$, with $k$ denoting a
real positive constant and $\bm{\sigma}=(\sigma_{1},\dots,\sigma_{p})^{T}$
representing the standard deviation function of the current observation
$\bm{X}_{n}$. From Equation (8), it is worth noting that, when ${\bm{X}}_{n}$
approaches ${\bm{Y}}_{n-1}$, i.e., $E_{nj}(t)$ tends to zero through a given
score function, the AMFEWMA chart is expected to perform like the MFEWMA
control chart based on the statistic $\tilde{\bm{Y}}_{n}$ defined in Equation
(7),. Instead, it is expected to behave as a Shewhart control chart when
$|E_{nj}(t)|$ is large. Hence, the calculation of the AMFEWMA statistic
requires the careful selection of the parameters $\lambda$ and $k$, which will
be addressed in Section 2.4. The functional extension of the classical
Hotelling’s statistic, based on the AMFEWMA statistic can be defined as

$V_{n}^{2}=\sum_{i=1}^{p}\sum_{j=1}^{p}\int_{\mathcal{T}}\int_{\mathcal{T}}Y_{ni}(s)K_{ij}^{*}(s,t)Y_{nj}(t)dsdt,$
(11)

where the function $K^{*}_{ij}$ is defined as

$K_{ij}^{*}(s,t)=\sum_{l=1}^{L}\frac{1}{\rho_{l}}\psi_{li}(s)\psi_{lj}(t),\quad
s,t\in\mathcal{T},i,j=1,\dots,p,$ (12)

and ${\bm{\psi}}_{l}=(\psi_{l1},\dots,\psi_{lp})^{T}$ denotes the
eigenfunctions corresponding to the eigenvalues $\rho_{l}$ of the asymptotic
covariance $\bm{C}(s,t)$ of ${\bm{Y}}_{n}$, calculated by applying the MFPCA
presented in Section 2.2. The number of principal components $L$ is usually
chosen such that the retained principal components explain a given fraction of
the total variability.

###  2.4 Design and implementation

The AMFEWMA control chart will be designed to signal when the statistic
$V_{n}^{2}$ given in Equation (11) exceeds the upper control limit $h$ chosen
as the smallest value achieving an IC ARL larger than or equal to a pre-
specified value, usually denoted by $ARL_{0}$. However, to calculate the
control limit $h$, we need to design the AMFEWMA statistics through the
selection of the parameters $\bm{\theta}=(\lambda,k)^{T}$ involved in the
score functions suggested in Equation (9) and Equation (10). In this paper,
the optimal parameters are obtained by elaborating on the idea of Capizzi and
Masarotto,  (2003) to provide, to the best possible extent, the AMFEWMA
control chart with good performance at different shift magnitudes.

First, a desired $ARL_{0}$, together with two shift values, say
$\bm{\delta}_{1}$ and $\bm{\delta}_{2}$, with
$\lVert\bm{\delta}_{1}\rVert\ll\lVert\bm{\delta}_{2}\rVert$, and a small
positive constant $\varepsilon$ (e.g., $\varepsilon=0.05$) must be chosen. In
general $\bm{\delta}_{1}$ and $\bm{\delta}_{2}$ can represent any “small” and
“large” split, respectively, as in Capizzi and Masarotto,  (2003). Then, the
optimal parameters $\bm{\theta}$, denoted by $\bm{\theta}^{*}$, are found as
the solutions of the following problem

$\displaystyle\bm{\theta}_{2}$
$\displaystyle=\operatorname*{argmin}_{\bm{\theta}}ARL(\bm{\delta}_{2},\bm{\theta})\quad\text{s.
t.}\quad ARL(0,\bm{\theta})=ARL_{0},$ (13) $\displaystyle\bm{\theta}^{*}$
$\displaystyle=\operatorname*{argmin}_{\bm{\theta}}ARL(\bm{\delta}_{1},\bm{\theta})\quad\text{s.
t.}\quad ARL(0,\bm{\theta})=ARL_{0}\quad\text{and}\quad
ARL(\bm{\delta}_{2},\bm{\theta})\leq(1+\varepsilon)ARL(\bm{\delta}_{2},\bm{\theta}_{2}),$
(14)

where $ARL(\bm{\delta},\bm{\theta})$ denotes the ARL achieved by the AMFEWMA
control chart with parameters $\bm{\theta}$ when the mean shift is equal to
$\bm{\delta}$. The constant $\varepsilon$ allows the AMFEWMA scheme to achieve
an ARL at $\bm{\delta}_{2}$ that is close to the minimum
$ARL(\bm{\delta}_{2},\bm{\theta}_{2})$ while searching for the minimum
$ARL(\bm{\delta}_{1},\bm{\theta})$ with respect to $\bm{\theta}$. Hereinafter,
we will denote the AMFEWMA control chart with parameters $\bm{\theta}^{*}$
chosen with the above procedure as $\text{AMFEWMA}^{*}$. It is worth
highlighting that the parameters included in the vector $\bm{\theta}$, and
thus $\bm{\theta}^{*}$, may vary with respect to the score function used.

To calculate the control limit $h$, $n_{seq}$ sequences of $n_{obs}$
observations are generated through a bootstrap procedure as proposed by Gandy
and Kvaløy,  (2013) from the Phase I sample. On each sequence, the AMFEWMA
statistic and the function $K^{*}_{ij}$ are calculated as defined in Equation
(8) and Equation (12). Then, the monitoring statistic $V_{n}^{2}$ is
calculated as defined in Equation (11), respectively. For each sequence, the
number of observations acquired up to the first signal of an OC state is
stored as a run length (RL), then the ARL is calculated as the average of the
RL values over the $n_{seq}$ sequences. The control limit $h$ is then chosen
to reach a pre-specified $ARL_{0}$.

Overfitting issues may arise when calculating $h$ (Kruger and Xie, , 2012)
because the asymptotic covariance of the AMFEWMA statistic $\bm{Y}_{n}$, used
in Equation (11) to calculate $V_{n}^{2}$, depends on $\bm{Y}_{n}$ itself.
Thus, to mitigate this problem, the control limit $h$ will be customarily
based on a random subset of the Phase I sample, referred to as the tuning set,
which is separate from the remaining one used to estimate the covariance of
${\bm{Y}}_{n}$, referred to as the training set.

In the Phase II monitoring, given the current sequence ${\bm{X}}_{n}$, the
AMFEWMA statistic ${\bm{Y}}_{n}$ and the corresponding monitoring statistic
${V}_{n}^{2}$ are obtained from Equation (8) and Equation (11), respectively,
according to the asymptotic covariance function $\bm{C}(s,t)$ estimated in
Phase I. An alarm signal is issued if ${V}_{n}^{2}$ exceeds the control limit
$h$ calculated in Phase I. A similar procedure is used in Phase II, where
$n^{II}_{seq}$ bootstrap sequences are generated from the Phase II data set,
and the ARL is obtained as the average of the RL values calculated for each
sequence.

##  3 Simulation study

The performance of the proposed AMFEWMA control chart in identifying mean
shifts in multivariate functional data is evaluated by means of an extensive
Monte Carlo simulation study. The data generation process, which is detailed
in Supplementary Material A, is performed as in Capezza et al., 2022b and is
inspired by the typical functional form of DRCs in a RSW process, as that
introduced in Section 1 and studied in Section 4. Without loss of generality,
multivariate functional quality data are generated with $p=5$ components. Two
OC scenarios are considered where the OC observations simulate typical
outlying DRCs, based on the work of Xia et al.,  (2019). In Scenario 1, OC
observations mimic a splash weld, which is also known as the expulsion
phenomenon and is caused by an excessive welding current flowing through the
electrodes. In Scenario 2, they mimic a phase shift of the peak time caused by
an overlay large electrode force. In each scenario, six increasing severity
levels $SL\in\\{0,1,2,3,4,5,6\\}$ are explored, where $SL=0$ means that the
process is IC. The severity level quantifies the extent to which Phase II data
set profiles deviate from their Phase I counterparts. A higher severity level
indicates a greater departure from the baseline and, consequently, a more
pronounced OC state for the process. Supplementary Material A defines the
severity levels, which can be visually interpreted in Figure A.1 and A.3 that
superimpose, on some IC observations, the shifted mean under each shift type,
for severity levels $SL\in\\{1,6\\}$.

The proposed chart is compared with two competing approaches: the EWMA control
chart for multivariate functional data, which is referred to as MFEWMA and
based on the Phase II monitoring scheme of Ren et al.,  (2019) with weighting
parameter $\lambda\in\\{0.1,0.2,0.3,0.5\\}$; and the classical multivariate
Shewhart control chart based on the Hotelling $T^{2}$ statistic applied to the
coefficients obtained from the MFPCA (Capezza et al., , 2023), referred to as
SHEWHART.

For each scenario and $SL$, 30 simulation runs are performed. Each run
considers a Phase I sample composed of 2500 observations, with training and
tuning samples of 1000 and 1500 observations, respectively. The Phase II data
set is composed of 200 sequences of i.i.d. observations with shift locations
equal to 100. The AMFEWMA control chart is implemented as detailed in Section
2 with $n_{seq}$ and $n_{obs}$ set equal to 500 and 300, respectively. The
parameter $L$ of Equation (12) is chosen such that the retained principal
components explain 90% of the total variability.

The AMFEWMA and the competing methods’ performance is reported in Table 1 and
Table 2 in terms of the estimated $ARL$, denoted by $\widehat{ARL}$, as a
function of $SL$ for Scenario 1 and Scenario 2, respectively. When the process
is IC, the $\widehat{ARL}$ should be as close as possible to $ARL_{0}$, which
is set equal to 20, while it should be as small as possible when the process
is OC. The AMFEWMA control chart is implemented for
$\lambda\in\\{0.1,0.2,0.3,0.5\\}$ and score function $\eta_{1}(\cdot)$ as
reported in Equation (9) with $k\in\\{2,3,4\\}$, whereas the
$\text{AMFEWMA}^{*}$ is obtained with optimal parameters $\bm{\theta}^{*}$
chosen through the algorithm presented in Section 2.4. In the latter, the
values of $\bm{\delta}_{1}$ and $\bm{\delta}_{2}$ are set as
$\bm{\delta}_{1}=\hat{\bm{\mu}}_{0}+0.5\hat{\bm{\sigma}}$ and
$\bm{\delta}_{2}=\hat{\bm{\mu}}_{0}+2\hat{\bm{\sigma}}$, where the sample mean
$\bm{\hat{\mu}_{0}}$ and standard deviation $\hat{\bm{\sigma}}$ functions are
estimated based on 100 bootstrap samples drawn from the Phase I sample.

In agreement with the literature, Table 1 and Table 2 reveal that, among the
competing approaches, the SHEWHART control chart achieves the best performance
for large $SL$s, whereas the MFEWMA with small $\lambda$ values turn out to be
the best approach for small $SL$s. From these tables, it is clear that the
parameters $\lambda$ and $k$ have a marked impact on the performance of the
AMFEWMA control chart in both scenarios. As an example, for $\lambda=0.2$ and
$k=2$, the proposed control chart, while being optimal for large shifts,
achieves unsatisfactory $\widehat{ARL}$ for small shifts, performing very
similarly to SHEWART. On the other hand, keeping fixed $\lambda=0.2$ but
setting $k=4$, the AMFEWMA performs similarly to the MFEWMA. In all cases, for
each $SL$, there is a combination of parameters, e.g., $\lambda=0.3$ and
$k=3$, for which the AMFEWMA performance is the best approach, or is very
close to it. Therefore, the problem of finding optimal parameters
$\bm{\theta}^{*}$ remains indispensable. The $\widehat{ARL}$s values achieved
by the proposed control chart with optimal parameters $\bm{\theta}^{*}$,
obtained through the score function $\eta_{1}$ in Equation (9) and referred to
as $\text{AMFEWMA}^{*}$, are reported in the last columns of Table 1 and Table
2 and indicate adequate performance for every OC condition considered. Indeed,
the $\text{AMFEWMA}^{*}$ performs similarly to the MFEWMA control chart with
small $\lambda$ for small $SL$s, and comparably to the SHEWHART control chart
for large $SL$s. This behavior is graphically shown in Figure 1, which reports
the estimated $\widehat{ARL}$ values for the SHEWHART, MFEWMA (with
$\lambda\in\\{0.1,0.5\\}$) and the $\text{AMFEWMA}^{*}$ control charts in
Scenario 1 and Scenario 2. From these figures, the $\widehat{ARL}$s for the
$\text{AMFEWMA}^{*}$ control chart appears always close to the smallest
estimated $ARL$ value, for each $SL$ and scenario considered.

Figure 1: The estimated $ARL$ values for the SHEWHART, MFEWMA (with
$\lambda\in\\{0.1,0.5\\}$) and the $\text{AMFEWMA}^{*}$ control charts in
Scenario 1 and Scenario 2.

These results are further demonstrated by means of the relative mean index
(RMI) suggested by Zou et al.,  (2012) and defined as

$RMI=\frac{1}{6}\sum_{SL=1}^{6}\frac{ARL_{SL}-sARL_{SL}}{sARL_{SL}},$

where, at given SL, $ARL_{SL}$ denotes the estimated $ARL$ value of the
considered control chart and $sARL_{SL}$ the smallest estimated $ARL$ value
among the competitors. The ratio $\frac{ARL_{SL}-sARL_{SL}}{sARL_{SL}}$ can be
regarded as a relative efficiency measure of the considered control chart with
respect to the best competitor, at given a $SL$, while the $RMI$ represents
its average over different SLs. The smaller the $RMI$, the better the overall
performance of the control chart. Table 3 shows the $RMI$ achieved by the
SHEWHART, the MFEWMA (with $\lambda\in\\{0.1,0.2,0.3,0.5\\}$), and the
$\text{AMFEWMA}^{*}$ control charts in Scenario 1 and Scenario 2. In both
scenarios, the $\text{AMFEWMA}^{*}$ confirms to outperform the competing
methods also in terms of $RMI$. It is worth noting that the SHEWHART control
chart and MFEWMA scheme with $\lambda=0.1$, which represent practitioner’s
standard choices, show very poor performance in both Scenario 1 and Scenario 2
and thus they should be used with care unless the size of the potential OC
condition is known a priori.

Table 1: The estimated $ARL$ values as a function of the severity level $SL$
for Scenario 1. The AMFEWMA control chart is implemented with
$\lambda\in\\{0.1,0.2,0.3,0.5\\}$ and $k\in\\{2,3,4\\}$. The smallest
$\widehat{ARL}$ in each row is highlighted in bold when $SL\neq 0$.

SL | SHEWHART | MFEWMA | AMFEWMA ($k=2$) | AMFEWMA ($k=3$) | AMFEWMA ($k=4$) | $\text{AMFEWMA}^{*}$
---|---|---|---|---|---|---
|  | $\lambda$ | $\lambda$ | $\lambda$ | $\lambda$ |
|  | 0.1 | 0.2 | 0.3 | 0.5 | 0.1 | 0.2 | 0.3 | 0.5 | 0.1 | 0.2 | 0.3 | 0.5 | 0.1 | 0.2 | 0.3 | 0.5 |
0 | 20.05 | 23.30 | 21.95 | 21.92 | 21.11 | 21.09 | 21.10 | 20.84 | 20.98 | 22.50 | 21.54 | 21.63 | 21.08 | 23.38 | 21.90 | 21.67 | 21.02 | 21.92
1 | 14.38 | 7.97 | 8.04 | 8.67 | 10.28 | 10.87 | 10.60 | 10.75 | 11.81 | 8.36 | 8.29 | 8.92 | 10.68 | 8.00 | 8.09 | 8.69 | 10.35 | 8.73
2 | 6.44 | 4.09 | 3.71 | 3.62 | 3.82 | 4.08 | 3.86 | 3.84 | 4.23 | 3.81 | 3.60 | 3.57 | 3.84 | 4.03 | 3.69 | 3.60 | 3.83 | 3.63
3 | 2.75 | 2.82 | 2.44 | 2.27 | 2.13 | 1.97 | 1.93 | 1.93 | 2.01 | 2.21 | 2.14 | 2.07 | 2.05 | 2.71 | 2.39 | 2.23 | 2.11 | 2.12
4 | 1.49 | 2.20 | 1.88 | 1.72 | 1.52 | 1.29 | 1.27 | 1.27 | 1.30 | 1.45 | 1.47 | 1.45 | 1.39 | 1.98 | 1.78 | 1.65 | 1.49 | 1.49
5 | 1.10 | 1.84 | 1.58 | 1.42 | 1.22 | 1.06 | 1.05 | 1.06 | 1.06 | 1.13 | 1.15 | 1.14 | 1.11 | 1.50 | 1.39 | 1.31 | 1.18 | 1.16
6 | 1.01 | 1.61 | 1.36 | 1.20 | 1.06 | 1.01 | 1.01 | 1.01 | 1.01 | 1.02 | 1.03 | 1.03 | 1.02 | 1.19 | 1.15 | 1.10 | 1.04 | 1.03

Table 2: The estimated $ARL$ values as a function of the severity level $SL$
for Scenario 2. The AMFEWMA control chart is implemented with
$\lambda\in\\{0.1,0.2,0.3,0.5\\}$ and $k\in\\{2,3,4\\}$. The smallest
$\widehat{ARL}$ in each row is highlighted in bold when $SL\neq 0$.

SL | SHEWHART | MFEWMA | AMFEWMA ($k=2$) | AMFEWMA ($k=3$) | AMFEWMA ($k=4$) | $\text{AMFEWMA}^{*}$
---|---|---|---|---|---|---
|  | $\lambda$ | $\lambda$ | $\lambda$ | $\lambda$ |
|  | 0.1 | 0.2 | 0.3 | 0.5 | 0.1 | 0.2 | 0.3 | 0.5 | 0.1 | 0.2 | 0.3 | 0.5 | 0.1 | 0.2 | 0.3 | 0.5 |
0 | 20.05 | 23.30 | 21.95 | 21.92 | 21.11 | 21.09 | 21.10 | 20.84 | 20.98 | 22.50 | 21.54 | 21.63 | 21.08 | 23.38 | 21.90 | 21.67 | 21.02 | 21.92
1 | 16.87 | 10.36 | 10.75 | 11.73 | 13.57 | 13.69 | 13.49 | 13.97 | 14.77 | 10.94 | 11.17 | 12.17 | 14.04 | 10.38 | 10.88 | 11.79 | 13.64 | 11.82
2 | 10.09 | 5.26 | 5.04 | 5.09 | 5.82 | 6.20 | 5.95 | 6.03 | 6.87 | 5.24 | 5.06 | 5.16 | 6.03 | 5.26 | 5.04 | 5.10 | 5.84 | 5.19
3 | 5.02 | 3.57 | 3.17 | 3.02 | 3.06 | 3.21 | 3.08 | 3.07 | 3.26 | 3.29 | 3.07 | 2.97 | 3.08 | 3.53 | 3.16 | 3.02 | 3.07 | 3.04
4 | 2.57 | 2.72 | 2.36 | 2.19 | 2.04 | 1.96 | 1.91 | 1.89 | 1.93 | 2.29 | 2.16 | 2.06 | 1.98 | 2.66 | 2.34 | 2.17 | 2.02 | 2.11
5 | 1.54 | 2.22 | 1.91 | 1.75 | 1.56 | 1.40 | 1.39 | 1.38 | 1.38 | 1.68 | 1.63 | 1.57 | 1.48 | 2.13 | 1.86 | 1.72 | 1.55 | 1.62
6 | 1.16 | 1.90 | 1.63 | 1.49 | 1.29 | 1.13 | 1.13 | 1.13 | 1.13 | 1.33 | 1.31 | 1.28 | 1.20 | 1.75 | 1.56 | 1.44 | 1.27 | 1.31

Table 3: The $RMI$ for the SHEWHART, MFEWMA (with
$\lambda\in\\{0.1,0.2,0.3,0.5\\}$), and the $\text{AMFEWMA}^{*}$ control
charts in Scenario 1 and Scenario 2.

| SHEWHART | MFEWMA | $\text{AMFEWMA}^{*}$
---|---|---|---
|  | $\lambda$ |
|  | 0.1 | 0.2 | 0.3 | 0.5 |
Scenario 1 | 1.88 | 2.36 | 1.32 | 0.88 | 0.59 | 0.27
Scenario 2 | 2.55 | 1.81 | 0.99 | 0.73 | 0.66 | 0.49

##  4 Case study

The case study on the SPM of an RSW process in automotive body-in-white
manufacturing, mentioned in the introduction, is presented to demonstrate the
practical applicability of the proposed AMFEWMA control chart as well as to
confirm its superior performance over competitors. Data analyzed are courtesy
of Centro Ricerche Fiat and are recorded at the Mirafiori Factory during lab
tests on different bodies of the same car model. Each body is characterized by
the same large number of spot welds and different characteristics, e.g., metal
sheet thickness and material and the welding time.

The case study focuses on a set of $p=10$ spot welds made by the same welding
gun on $n=1340$ car bodies. That is, for each car body, we consider the
multivariate functional quality characteristic represented by the vector of
$p=10$ DRCs relative to the same ten spot welding points. For each DRC, raw
measurements were collected at a regular grid of points equally spaced by 1
ms, normalized on the time domain $[0,1]$. To provide a sketch of the data,
the 10-dimensional vectors of DRCs corresponding to 20 random car bodies out
of $n=1340$ are displayed in Figure 2.

Figure 2: A sketch of DRCs pertaining to the ten spot welds measured on 20 car
bodies randomly selected from the dataset analyzed in the case study.

Even though this is not meant to be the routine application, to allow for
further comparison with competing control charting schemes based on this case
study, the Phase I sample has been formed by selecting 568 multivariate
profiles observed immediately after the renewal, which the electrodes are
routinely subject to for neutralizing the wear effect that is known, as
mentioned in the introduction, to impact on the quality of the final welded
point. The remaining 772 observations are grouped into three data sets of 220,
368, and 184 observations, corresponding to three incremental electrode wear
levels, namely wear levels 1, 2, and 3, respectively. These three data sets
are then used as Phase II observations to evaluate the proposed chart
performance with respect to different levels of shift connected with wear
level. It is important to emphasize that this is a special case where
information regarding electrode wear is available because the data are
acquired during lab tests, while such information is commonly out of reach.
Hence, a control chart capable of accommodating varying degrees of mean shift,
contingent upon different levels of wear, without prior knowledge of the wear
levels, can prove immensely advantageous.

The AMFEWMA control chart is implemented as in Section 3 and the 568
observations of the Phase I sample are randomly split into a training and a
tuning set, with equal sample sizes. The optimal parameters $\lambda^{*}=0.5$
and $k^{*}=4$ used in the $\text{AMFEWMA}^{*}$ are obtained from Equation (14)
through the score function $\eta_{1}$ reported in Equation (9). The bootstrap
method presented in Section 2.4 is applied to each of the three Phase II data
sets to generate 200 sequences of 200 observations.

Figure 3 shows the $\text{AMFEWMA}^{*}$ application to the Phase II monitoring
of the three Phase II data sets corresponding to three increasing wear levels.
In each plot, the first 100 observations are randomly sampled from the tuning
set, whereas the Phase II observations for each wear level considered are
displayed on the right of the dashed vertical line. From this figure, the
$\text{AMFEWMA}^{*}$ scheme clearly appears to adequately detect the mean
shift under all shift severity levels.

Figure 3: $\text{AMFEWMA}^{*}$ control chart under wear levels 1, 2, and 3. In
each plot, the first 100 observations are randomly sampled from the tuning
set, while the Phase II observations are reported after the dashed vertical
line.

The AMFEWMA performance is further compared in this case study with the
competing methods, in terms of ARL, and the results are shown in Table 4.

Table 4:  ARL calculated on the three Phase II data sets (corresponding to wear level 1, 2, and 3, respectively) in the case study, for the proposed AMFEWMA and the competing control charts. In bold, the lowest ARL value for each row is reported. | SHEWHART | MFEWMA |  $\text{AMFEWMA}^{*}$ $\lambda=0.5,k=4$
---|---|---|---
Wear level |  |  $\lambda$=0.1 |  $\lambda$=0.2 |  $\lambda$=0.3 |  $\lambda$=0.5
1 | 5.68 | 3.69 | 3.095 | 2.91 | 2.905 | 2.92
2 | 1.25 | 2.185 | 1.715 | 1.52 | 1.315 | 1.3
3 | 1.025 | 1.695 | 1.415 | 1.295 | 1.12 | 1.085

For wear level 1, the MFEWMA control chart with $\lambda=0.5$ achieves the
lowest ARL value, whereas in the other cases, the SHEWHART control chart
achieves the best performance. It is easy to note that, as expected, the
performance of the AMFEWMA chart is always very close to the best-performing
competitor. Therefore, it can be regarded as the best overall choice in all
cases where there is no prior information about the severity of the OC
scenario.

##  5 Conclusions

A new control chart for the statistical process control of multivariate
functional data is proposed and referred to as adaptive multivariate
functional exponentially weighted moving average (EWMA) control chart
(AMFEWMA). The AMFEWMA chart is the first EWMA scheme that is able to
adaptively monitor a multivariate functional quality characteristic. Whereas,
methods already presented in the literature either apply the conventional EWMA
approach to multivariate or univariate functional data. Through an extensive
Monte Carlo simulation, the performance of the AMFEWMA chart is compared with
two competing methods, which proved the AMFEWMA to outperform competing
methods in the correct identification of out-of-control (OC) scenarios when
the magnitude of the mean shift is unknown. For small shifts, the AMFEWMA
performs as the conventional MFEWMA control chart, whereas, for large shifts,
it performs as the Shewart control chart. The ability of the proposed method
to deal with multivariate functional data and effectively signal process mean
shifts in a wide range of severity levels is of paramount importance to
protect modern industrial processes where many functional quality
characteristics are routinely available in line against OC scenarios with
different severity levels.

## Supplementary Materials

The Supplementary Materials contain additional details about the data
generation process in the simulation study as well as the R code to reproduce
graphs and results over competing methods.

## Data Availability Statement

The authors confirm that the data supporting the findings of this study are
available within the article and its supplementary materials.

## Acknowledgements

The authors are extremely grateful to CRFWCM R&I (World Class Manufacturing
Research and Innovation) sites of Orbassano and CampusMelfi, for the access to
experimental data and the technological insights in the interpretation of the
results.

## Funding

The contribution of A. Lepore and F. Centofanti has been economically
supported by Piano Nazionale di Ripresa e Resilienza (PNRR) - Missione 5
Componente 2, Investimento 1.3-D.D. 1551.11-10-2022, PE00000004 within the
Extended Partnership MICS (Made in Italy - Circular and Sustainable). The
activity of B. Palumbo is supported by the MOST - Sustain- able Mobility
National Research Center and received funding from the European Union Next-
GenerationEU (PIANO NAZIONALE DI RIPRESA E RESILIENZA (PNRR) - MIS- SIONE 4
COMPONENTE 2, INVESTIMENTO 1.4 - D.D. 1033 17/06/2022, CN00000023). This
manuscript reflects only the authors’ views and opinions, neither the European
Union nor the European Commission can be considered responsible for them.

## References

  * (1) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2022a).  Functional regression control chart for monitoring ship CO2 emissions.  Quality and Reliability Engineering International, 38(3):1519–1537.
  * Capezza et al., (2023) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2023).  funcharts: Control charts for multivariate functional data in R.  Journal of Quality Technology.
  * Capezza et al., (2021) Capezza, C., Centofanti, F., Lepore, A., and Palumbo, B. (2021).  Functional clustering methods for resistance spot welding process data in the automotive industry.  Applied Stochastic Models in Business and Industry, 37(5):908–925.
  * (4) Capezza, C., Centofanti, F., Lepore, A., and Palumbo, B. (2022b).  Robust multivariate functional control charts.  arXiv preprint arXiv:2207.07978.
  * Capezza et al., (2020) Capezza, C., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2020).  Control charts for monitoring ship operating conditions and CO2 emissions based on scalar-on-function regression.  Applied Stochastic Models in Business and Industry, 36(3):477–500.
  * Capizzi and Masarotto, (2003) Capizzi, G. and Masarotto, G. (2003).  An adaptive exponentially weighted moving average control chart.  Technometrics, 45(3):199–207.
  * Cardot et al., (2003) Cardot, H., Ferraty, F., and Sarda, P. (2003).  Spline estimators for the functional linear model.  Statistica Sinica, pages 571–591.
  * Centofanti et al., (2022) Centofanti, F., Lepore, A., Kulahci, M., and Spooner, M. P. (2022).  Real-time monitoring of functional data.  arXiv preprint arXiv:2205.06256.
  * Centofanti et al., (2021) Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2021).  Functional regression control chart.  Technometrics, 63(3):281–294.
  * Chan and Zhang, (2000) Chan, L. K. and Zhang, J. (2000).  Some issues in the design of EWMA charts.  Communications in Statistics-Simulation and Computation, 29(1):201–217.
  * Chiou et al., (2014) Chiou, J.-M., Chen, Y.-T., and Yang, Y.-F. (2014).  Multivariate functional principal component analysis: A normalization approach.  Statistica Sinica, pages 1571–1596.
  * Chou et al., (2014) Chou, S.-H., Chang, S. I., and Tsai, T.-R. (2014).  On monitoring of multiple non-linear profiles.  International Journal of Production Research, 52(11):3209–3224.
  * Colosimo and Pacella, (2010) Colosimo, B. M. and Pacella, M. (2010).  A comparison study of control charts for statistical monitoring of functional data.  International Journal of Production Research, 48(6):1575–1601.
  * Crowder, (1989) Crowder, S. V. (1989).  Design of exponentially weighted moving average schemes.  Journal of Quality technology, 21(3):155–162.
  * De Boor, (1978) De Boor, C. (1978).  A practical guide to splines, volume 27.  springer-verlag New York.
  * Fassò et al., (2016) Fassò, A., Toccu, M., and Magno, M. (2016).  Functional control charts and health monitoring of steam sterilizers.  Quality and Reliability Engineering International, 32(6):2081–2091.
  * Gan, (1993) Gan, F. (1993).  An optimal design of EWMA control charts based on median run length.  Journal of Statistical Computation and Simulation, 45(3-4):169–184.
  * Gandy and Kvaløy, (2013) Gandy, A. and Kvaløy, J. T. (2013).  Guaranteed conditional performance of control charts via bootstrap methods.  Scandinavian Journal of Statistics, 40(4):647–668.
  * Grasso et al., (2014) Grasso, M., Colosimo, B. M., and Pacella, M. (2014).  Profile monitoring via sensor fusion: the use of PCA methods for multi-channel data.  International Journal of Production Research, 52(20):6110–6135.
  * Grasso et al., (2017) Grasso, M., Colosimo, B. M., and Tsung, F. (2017).  A phase I multi-modelling approach for profile monitoring of signal data.  International Journal of Production Research, 55(15):4354–4377.
  * Grasso et al., (2016) Grasso, M., Menafoglio, A., Colosimo, B. M., and Secchi, P. (2016).  Using curve-registration information for profile monitoring.  Journal of Quality Technology, 48(2):99–127.
  * Green and Silverman, (1993) Green, P. J. and Silverman, B. W. (1993).  Nonparametric regression and generalized linear models: a roughness penalty approach.  Crc Press.
  * Happ and Greven, (2018) Happ, C. and Greven, S. (2018).  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association, 113(522):649–659.
  * Haq and Khoo, (2019) Haq, A. and Khoo, M. B. (2019).  An adaptive multivariate ewma chart.  Computers & Industrial Engineering, 127:549–557.
  * Jin and Shi, (1999) Jin, J. and Shi, J. (1999).  Feature-preserving data compression of stamping tonnage information using wavelets.  Technometrics, 41(4):327–339.
  * Jolliffe and Cadima, (2016) Jolliffe, I. T. and Cadima, J. (2016).  Principal component analysis: a review and recent developments.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2065):20150202.
  * Jones et al., (2021) Jones, C. L., Abdel-Salam, A.-S. G., and Mays, D. (2021).  Practitioners guide on parametric, nonparametric, and semiparametric profile monitoring.  Quality and Reliability Engineering International, 37(3):857–881.
  * Jones et al., (2001) Jones, L. A., Champ, C. W., and Rigdon, S. E. (2001).  The performance of exponentially weighted moving average charts with estimated parameters.  Technometrics, 43(2):156–167.
  * Knoth, (2007) Knoth, S. (2007).  Accurate ARL calculation for EWMA control charts monitoring normal mean and variance simultaneously.  Sequential Analysis, 26(3):251–263.
  * Kokoszka and Reimherr, (2017) Kokoszka, P. and Reimherr, M. (2017).  Introduction to functional data analysis.  Chapman and Hall/CRC.
  * Kruger and Xie, (2012) Kruger, U. and Xie, L. (2012).  Statistical monitoring of complex multivatiate processes: with applications in industrial process control.  John Wiley & Sons.
  * Lowry et al., (1992) Lowry, C. A., Woodall, W. H., Champ, C. W., and Rigdon, S. E. (1992).  A multivariate exponentially weighted moving average control chart.  Technometrics, 34(1):46–53.
  * Lucas and Saccucci, (1990) Lucas, J. M. and Saccucci, M. S. (1990).  Exponentially weighted moving average control schemes: properties and enhancements.  Technometrics, 32(1):1–12.
  * Mahmoud and Zahran, (2010) Mahmoud, M. A. and Zahran, A. R. (2010).  A multivariate adaptive exponentially weighted moving average control chart.  Communications in Statistics—Theory and Methods, 39(4):606–625.
  * Maleki et al., (2018) Maleki, M. R., Amiri, A., and Castagliola, P. (2018).  An overview on recent profile monitoring papers (2008–2018) based on conceptual classification scheme.  Computers & Industrial Engineering, 126:705–728.
  * Manladan et al., (2017) Manladan, S., Yusof, F., Ramesh, S., Fadzil, M., Luo, Z., and Ao, S. (2017).  A review on resistance spot welding of aluminum alloys.  The International Journal of Advanced Manufacturing Technology, 90(1):605–634.
  * Martín et al., (2014) Martín, Ó., Pereda, M., Santos, J. I., and Galán, J. M. (2014).  Assessment of resistance spot welding quality based on ultrasonic testing and tree-based techniques.  Journal of Materials Processing Technology, 214(11):2478–2487.
  * Menafoglio et al., (2018) Menafoglio, A., Grasso, M., Secchi, P., and Colosimo, B. M. (2018).  Profile monitoring of probability density functions via simplicial functional pca with application to image data.  Technometrics, 60(4):497–510.
  * Montgomery, (2019) Montgomery, D. C. (2019).  Introduction to statistical quality control.  John Wiley & Sons.
  * Noorossana et al., (2011) Noorossana, R., Saghaei, A., and Amiri, A. (2011).  Statistical analysis of profile monitoring.  John Wiley & Sons.
  * Paynabar et al., (2016) Paynabar, K., Zou, C., and Qiu, P. (2016).  A change-point approach for phase-I analysis in multivariate profile monitoring and diagnosis.  Technometrics, 58(2):191–204.
  * Qiu, (2014) Qiu, P. (2014).  Introduction to statistical process control.  CRC press.
  * R Core Team, (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Ramsay and Silverman, (2005) Ramsay, J. and Silverman, B. (2005).  Functional data analysis.  Springer.
  * Ren et al., (2019) Ren, H., Chen, N., and Wang, Z. (2019).  Phase-II monitoring in multichannel profile observations.  Journal of Quality Technology, 51(4):338–352.
  * Serel, (2009) Serel, D. A. (2009).  Economic design of EWMA control charts based on loss function.  Mathematical and Computer Modelling, 49(3-4):745–759.
  * Wang et al., (2018) Wang, Y., Mei, Y., and Paynabar, K. (2018).  Thresholded multivariate principal component analysis for phase I multichannel profile monitoring.  Technometrics, 60(3):360–372.
  * Xia et al., (2019) Xia, Y.-J., Su, Z.-W., Li, Y.-B., Zhou, L., and Shen, Y. (2019).  Online quantitative evaluation of expulsion in resistance spot welding.  Journal of Manufacturing Processes, 46:34–43.
  * Zhang and Senkara, (2011) Zhang, H. and Senkara, J. (2011).  Resistance welding: fundamentals and applications.  CRC press.
  * Zou et al., (2012) Zou, C., Ning, X., and Tsung, F. (2012).  Lasso-based multivariate linear profile monitoring.  Annals of operations research, 192(1):3–19.