[orcid=0000-0003-4222-9631]

1]organization=Technical University of Denmark, Department of Applied
Mathematics and Computer Science, city=Kongens Lyngby, country=Denmark

2]organization=Luleå University of Technology, Department of Business
Administration, Technology and Social Sciences, city=Luleå, country=Sweden

[orcid=0000-0002-8739-5310]

[1]

[orcid=0000-0003-1036-8127]

[orcid=0000-0002-4885-1635]

3]organization=University of Naples Federico II, Department of Industrial
Engineering, city=Naples, country=Italy

[cor1]Corresponding author

# Functional neural network control chart

Murat Kulahci<EMAIL_ADDRESS>[  [  Antonio Lepore<EMAIL_ADDRESS>Biagio
Palumbo<EMAIL_ADDRESS>Gianluca Sposito<EMAIL_ADDRESS>[

###### Abstract

In many Industry 4.0 data analytics applications, quality characteristic data
acquired from manufacturing processes are better modeled as functions, often
referred to as profiles. In practice, there are situations where a scalar
quality characteristic, referred to also as the response, is influenced by one
or more variables in the form of functional data, referred to as functional
covariates. To adjust the monitoring of the scalar response by the effect of
this additional information, a new profile monitoring strategy is proposed on
the residuals obtained from the functional neural network, which is able to
learn a possibly nonlinear relationship between the scalar response and the
functional covariates. An extensive Monte Carlo simulation study is performed
to assess the performance of the proposed method with respect to other control
charts that appeared in the literature before. Finally, a case study in the
railway industry is presented with the aim of monitoring the heating,
ventilation and air conditioning systems installed onboard passenger trains.

###### keywords:

Quality control Machine learning Data science

##  1 Introduction

In many modern production processes, advanced data acquisition systems
generate massive amounts of data in the form of curves or surfaces varying
over a continuum, such as time or space. Such data can be best modeled as
functions that are more generally defined multidimensionally, and are usually
referred to as profiles or functional data (Ramsay and Silverman,, 2005;
Menafoglio and Secchi,, 2017; Hsing and Eubank,, 2015; Ferraty,, 2006). The
statistical monitoring of a process that is best characterized by functional
data is known as profile monitoring. Some examples of profile monitoring
applications can be found in the works of Jin and Shi, (1999); Woodall et al.,
(2004); Zou et al., (2007); Williams et al., (2007); Colosimo and Pacella,
(2010); Saghaei et al., (2013). As an application area within statistical
process monitoring (SPM) (Stuart et al.,, 1996), profile monitoring aims at
the detection of special causes of variation acting on the process, which, in
such cases, is said to be out of control (OC). Otherwise, the process is said
to be affected by only random causes of variation and hence in control (IC).
More specifically, profile monitoring focuses on testing the stability of the
functional relationship between the quality characteristic of interest, also
referred to as the response, and one or more exploratory variables, referred
to as covariates. That is, the quality characteristic of interest is monitored
conditionally on the covariate levels. The regression control chart (Mandel,,
1969) is the first scheme proposed in statistical literature to address this
issue by monitoring the residuals of the regression of the quality
characteristic on the covariates. In this way, the monitoring strategy is able
to consider the variance explained by the covariates and to leverage their
additional information to improve the power of the monitoring scheme expressed
as the probability of correctly detecting a shift in the process being
monitored (Montgomery,, 2020). Recent and more advanced extensions of this
idea allow the covariates or the response to be described by functional data
themselves, such as the functional regression control chart (FRCC) framework
proposed by Centofanti et al., (2021) where functional models (Morris,, 2015)
can be used to map the influence of one or more functional covariates on a
scalar or functional response. However, all the instances that can be
reconducted to the FRCC framework, such as Capezza et al., (2020); Capezza et
al., 2022a ; Capezza et al., 2022b ; Centofanti et al., (2022), are
implemented by means of functional linear models, only.

In the meantime, neural networks (NNs) and deep learning (DL) techniques have
been receiving increasing attention in time series, computer vision, speech
recognition, and genetics. However, despite their success, their use for
functional data is not completely investigated. Rossi et al., (2005) discussed
how to incorporate functional pre-processing, e.g., functional principal
component analysis, into a multilayer perceptron (MLP), which is a particular
type of NN. Rossi et al., (2002); Conan-Guez and Rossi, (2002) introduced the
functional multilayer perceptron (FMLP), extending the MLP to functional data
through a proper reparameterization of the weight matrices. In particular, the
former proposes an FMLP that simply requires a discretization of the
functional covariates and, thus it cannot directly handle profiles with a
different number of observations or measured at different time points. To
overcome this issue, the latter uses the basis expansions (Wang et al.,, 2016)
of the functional covariates. The theoretical properties of the FMLP are
firstly studied by Rossi and Conan-Guez, (2006), where the authors revealed
that an FMLP is a universal approximator following the definition of Hornik et
al., (1989). FMLP has also been further investigated by Wang et al., 2019a ;
Wang et al., 2019b .

Another common approach to modeling functional data is through convolutional
(CNN) (LeCun et al.,, 1998) and recurrent NN (RNN) (Rumelhart et al.,, 1986).
In particular, long short-term memory variant (Hochreiter and Schmidhuber,,
1997) is the most popular, due to the ability to recognize patterns for a long
duration of time. However, these may fail to learn the underlying smoothness
of the functional data from raw noisy measurements, which are very common in
real-data applications. Yao et al., (2021) proposed an adaptive functional NN
(AdaFNN) that is a NN specifically designed for functional data through the
definition of a new basis layer which implements a micro NN (Lin et al.,,
2013) to directly learn the most relevant base functions to represent the
response value, thus avoiding a pre-specified choice. The authors demonstrate
the superior performance of their proposal over the MLPs through an extensive
simulation study, but they do not compare their proposal with any of the
functional regression models already present in the functional data (FDA)
literature. Additionally, even though the authors claimed that the functional
coefficients learned by AdaFNN are interpretable, it is not true in general as
each micro NN learns a different functional coefficient and it is impossible
to priorly know how to combine them and resemble the true functional
coefficients. Recently, Thind et al., (2023) introduced a new DL architecture,
known as functional neural network (FNN), which allows for deep architectures
for scalar responses with multiple functional and scalar covariates through
the definition of smooth weight functions. The functional weight can be
visualized during the training, increasing the interpretability of the FNN
while maintaining the nonlinear predictive power of traditional NNs.
Furthermore, the number of parameters of the FNN can be lower than the amount
needed in traditional MLPs, CNNs, and RNNs.(see Section 2.2 for more details).
Through simulated and benchmark data sets, Thind et al., (2023) demonstrated
that FNN outperforms state-of-the-art DL architectures, functional linear
models, and a number of other multivariate methods, e.g., least square
regression and random forest, in terms of prediction accuracy (Hastie et al.,,
2009). FNN is also proven to be a universal approximator (Cybenko,, 1989),
that is it can be used to learn any continuous function to a desired degree of
accuracy. In recent years, different DL-based have been proposed for the SPM
of multivariate processes and functional data, and some relevant examples can
be found in Stuart et al., (1996); Sergin and Yan, (2021); Pacella and
Semeraro, (2011); Chen et al., (2020); Howard et al., (2018); Yeganeh et al.,
(2022); Yeganeh and Shadman, (2021).

In this paper, we propose a novel control chart, named FNN control chart
(FNNCC), that exploits the FNN architecture advantages for the monitoring of a
scalar quality characteristic when functional and scalar covariates are
available. The proposed FNNCC can be regarded as an implementation of the FRCC
framework where the influence of functional or scalar covariates on a scalar
response does not need to be necessarily linear as in Capezza et al., (2020);
Capezza et al., 2023b , which is based on a linear scalar-on-function (SOF)
regression model (Reiss et al.,, 2017).

The article is structured as follows. Section 2 introduces some background on
the SOF regression model and FNNs. Then, we describe the proposed FNNCC in
detail. An extensive Monte Carlo simulation study is performed in Section 3 to
quantify the FNNCC OC average run length ($ARL_{1}$) at a given IC average run
length ($ARL_{0}$) (Qiu,, 2013) in identifying a mean shift in the scalar
response in the presence or absence of drifts in the covariate mean and to
compare it with other competing control charting schemes that have already
appeared in the literature before. The practical applicability of the proposed
method is illustrated in Section 4 by means of a case study in the monitoring
of heating, ventilation and air conditioning (HVAC) systems installed on
passenger trains (Lepore et al., 2022b, ). The HVAC data set, courtesy of the
rail transportation company Hitachi Rail Italy, and the analysis code are
available online at https://github.com/unina-sfere/FNNCC. The final section
contains concluding remarks and outlines new directions for future research.
Additionally, in Appendix A, we implement two other NN-based monitoring
strategies and compare them to the proposed control chart. All computations
and plots have been obtained using the programming language R (R Core Team,,
2021).

##  2 Methodology

In this section, we briefly review the SOF linear regression model in Section
2.1. Then, we introduce the FNN in Section 2.2 and provide details of the
monitoring strategy based on the FNN in Section 2.3.

###  2.1 Functional linear regression

We begin by providing the necessary notation used in the article, followed by
a summary of the linear model for SOF regression before describing the
proposed methodology.

Let $y_{i}$ and
$\boldsymbol{X}_{i}=\left\\{\left(X_{i1},\ldots,X_{iP}\right)^{\top}\right\\}$,
$i=1,\ldots,n$, denote $n$ observations of a scalar response covariates and a
vector of $P$ functional variables, respectively. $\boldsymbol{X}_{i}$ is a
random element that takes values in the Hilbert space
$L^{2}(\mathcal{T})^{P}$, i.e., $X_{i,1},\ldots X_{i,P}$ belong to
$L^{2}(\mathcal{T})$, the space of square-integrable functions defined on the
compact interval $\mathcal{T}$. We also assume that $\mathbf{X}_{i}$ is fully
observed (Kokoszka and Reimherr,, 2017) i.e., is densely observed on a set of
discrete grid points. The aim is to learn the mapping
$G:\mathcal{L}^{2}(\mathcal{T})\times$
$\ldots\times\mathcal{L}^{2}(\mathcal{T})\rightarrow\mathcal{R}$ from the
functional covariates $\mathbf{X}_{i}$ to the scalar response $y_{i}$

$y_{i}=G\left(\mathbf{X}_{i}\right)+\epsilon_{i}\quad i=1,\ldots,n,$ (1)

where $\epsilon_{i}$ is a scalar error term. The SOF regression has been
extensively studied in the FDA literature and, in this section, we give a
succinct description of the SOF linear regression model used by Capezza et
al., (2020), defined as

$y_{i}=\alpha+\sum_{p=1}^{P}\int_{\mathcal{T}}\beta_{p}(t)X_{ip}(t)dt+\epsilon_{i}\quad
i=1,\ldots,n,$ (2)

where $\alpha\in\mathbb{R}$ is the scalar intercept,
$\boldsymbol{\beta}=\left(\beta_{1},\ldots,\beta_{p},\ldots,\beta_{P}\right)^{\top}\in
L^{2}(\mathcal{T})^{P}$ are the functional coefficients to be estimated and
$\epsilon_{i}$ are the error terms, which are assumed to be independent and
identically distributed normal random variables with mean zero and variance
$\sigma^{2}$. Comparing Equation (1) and Equation (2), it is trivial to note
that in the linear SOF regression model
$G\left(\mathbf{X}_{i}\right)=\alpha+\sum_{p=1}^{P}\int_{\mathcal{T}}\beta_{p}(t)X_{ip}(t)dt$.
Without loss of generality, the functional covariates are assumed to be
empirically standardized (Chou et al.,, 2014), that is each covariate
observation is standardized by subtracting pointwise the corresponding sample
mean and dividing the result by the relative standard deviation function. The
coefficients $\alpha$ and $\boldsymbol{\beta}$ in Equation (2) can be
estimated by solving the following minimization problem

$\min_{\alpha\in\mathbb{R},\boldsymbol{\beta}\in
L^{2}(\mathcal{T})^{P}}\sum_{i=1}^{n}\left(y_{i}-\alpha-\sum_{r=1}^{P}\int_{\mathcal{T}}\beta_{p}(t)X_{ip}(t)dt\right)^{2}.$
(3)

Because of the infinite dimensionality of the functional data, the above
minimization problem is not well-posed and the model cannot be estimated using
the least squares approach James et al., (2013) directly. However, being
square integrable, the functional covariates can be represented through the
Karhunen-Loéve expansion as follows

$\boldsymbol{X}_{i}(t)=\sum_{m=1}^{\infty}\xi_{im}\boldsymbol{\psi}_{m}(t),\quad
t\in\mathcal{T},\quad i=1\ldots,n,$ (4)

where the
$\boldsymbol{\psi}_{m}=\left(\psi_{m1},\ldots,\psi_{mP}\right)_{m\in\mathbb{N}}\in
L^{2}(\mathcal{T})^{P}$ are the multivariate functional principal components
(MFPCs) defined as the eigenfunctions of the covariance function
$\mathbf{c}(s,t)$ of the multivariate functional data i.e. they are the
solutions to the equation

$\int_{\mathcal{T}}\mathbf{c}(s,t)\boldsymbol{\psi}_{m}(s)ds=\lambda_{m}\boldsymbol{\psi}_{m}(t),\quad\forall
t\in\mathcal{T},$ (5)

where $\lambda_{m}$ are called the eigenvalues of $\mathbf{c}(s,t)$. The
eigenvalues $\lambda_{m}$, and so the corresponding eigenfunctions
$\boldsymbol{\psi}_{m}(t)$, are arranged in non-increasing order
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$. The eigenfunctions
$\boldsymbol{\psi}_{m}(t)$ are by construction such that

$\sum_{p=1}^{P}\int_{\mathcal{T}}\psi_{m_{1},p}(t)\psi_{m_{2},p}(t)dt=\begin{cases}1,&\mbox{if
}m_{1}=m_{2}\\\ 0,&\mbox{if }m_{1}\neq m_{2},\end{cases}$ (6)

that is, they form an orthonormal basis of $L^{2}(\mathcal{T})^{P}$. The terms
$\xi_{im}$ in Equation (4) are called the MFPC scores, or simply scores, and
defined as

$\xi_{im}=\sum_{p=1}^{P}\int_{\mathcal{T}}X_{ip}(t)\psi_{mp}(t)dt\quad
i=1\ldots,n.$ (7)

It can be shown that $\text{E}(\xi_{m})=0$,
$\text{E}(\xi_{m}^{2})=\lambda_{m}$ and $\text{E}(\xi_{m_{1}},\xi_{m_{2}})=0$
when $m_{1}\neq m_{2}$. The decomposition in Equation (4) is optimal in the
sense that, for each finite $M\in N$, $\boldsymbol{X}_{i}(t)$ is best
approximated by
$\mathbf{\hat{X}}_{i}=(\hat{X}_{i1},\ldots,\hat{X}_{ip},\ldots,\hat{X}_{iP})^{T}$,
obtained by truncating the Karhunen-Loéve expansion, as follows

$\hat{X}_{ip}(t)=\sum_{m=1}^{M}\xi_{im}{\psi}_{mp}(t),\quad
t\in\mathcal{T},\quad p=1,\ldots,P,\quad i=1\ldots,n.$ (8)

It can be proved that, upon using $\mathbf{\hat{X}}_{i}$ in place of
$\mathbf{X}_{i}$, the functional coefficient $\boldsymbol{\beta}$ in Equation
(2) can be replaced with
$\hat{\boldsymbol{\beta}}=\left(\hat{\beta}_{1},\ldots,\hat{\beta}_{p},\ldots,\hat{\beta}_{P}\right)^{\top}\in
L^{2}(\mathcal{T})^{P}$ by using the same truncated basis expansion, that is

$\hat{\beta}_{p}(t)=\sum_{m=1}^{M}b_{m}\psi_{mp}(t),\quad
t\in\mathcal{T},\quad p=1,\ldots,P,$ (9)

where $b_{1},\ldots,b_{M}$ are the basis coefficients. Upon using the
approximation in Equations (8) and Equation (9), Equation (2) can be rewritten
as

$y_{i}=\alpha+\sum_{m=1}^{M}\xi_{im}b_{m}+\varepsilon_{i}^{*},\quad
i=1,\ldots,n.$ (10)

Hence, in this form, the model parameters, namely the intercept $\alpha$ and
the coefficient $b_{m}$, $m=1,\ldots,M$, can be estimated by using the least
square approach, as $\hat{\alpha}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$ and
$\hat{b}_{m}=\sum_{i=1}^{n}y_{i}\xi_{im}/\sum_{i=1}^{n}\xi_{im}^{2}$,
respectively. Accordingly, the least square prediction of the scalar response
$y_{i}$ can be obtained as

$\hat{y}_{i}=\hat{\beta}_{0}+\sum_{p=1}^{P}\int_{\mathcal{T}}X_{ip}(t)\hat{\beta}_{p}(t)dt,\quad
i=1,\ldots,n$ (11)

As in the multivariate setting, the number $M$ can be chosen such that the
retained MFPCs
$\boldsymbol{\psi}_{m}=\left(\psi_{m1},\ldots,\psi_{mP}\right)_{m=1,\ldots,M}$
explain at least a given percentage, say 80%, of the total variability. In
this paper, we use a different strategy based on the final model’s predictive
ability. That is, the $M$ MFPCs to be retained are the first $M$ that achieve
a given reduction in the prediction sum of squares statistic, defined as
$\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{[i]}\right)^{2}$, where $\hat{y}_{[i]}$ is
the $i$-th fitted value of the scalar response based on the SOF regression
model with the $i$-th observation removed from the data set used to fit the
linear model. In this way, this strategy ensures that the MFPCs with a small
predictive ability are not retained in the SOF model. More details on this
problem can be found in (Jolliffe and Cadima,, 2016, p. 173-177).

###  2.2 Functional neural network

NNs are computational models inspired by the structure and functioning of the
human brain and consisting of interconnected artificial neurons, also known as
nodes or units. The latter are organized into layers, typically an input
layer, one or more hidden layers, and an output layer. The input layer
receives the initial input data, and the output layer produces the final
output or prediction. The hidden layers are intermediate layers between the
input and output layers and play a crucial role in learning complex
representations of the data. Each neuron takes multiple inputs from the
previous layer, performs a weighted sum, and then applies a nonlinear function
to produce an output. Let $n^{(j)}$ and $\mathbf{h}^{(j)}$ be the number of
neurons in the $j$th hidden layer and the output of the $j$th hidden layer,
respectively. More formally, $\mathbf{h}^{(j)}$ is defined as
$\mathbf{h}^{(j)}=g\left(\mathbf{W}^{(j)}\boldsymbol{h}^{(j-1)}+\boldsymbol{b}^{(j)}\right)$,
where $\mathbf{h}^{(j-1)}\in R^{n^{(j-1)}}$ represents the output of the
previous $(j-1)$th layer, $\mathbf{W}^{(j)}$ is a $n^{(j)}\times n^{(j-1)}$
weight matrix and $\boldsymbol{b}^{(j)}\in R^{n^{(j)}}$ is the intercept,
often referred to as the bias in the machine learning field. The function
$g:\mathbb{R}^{n^{(j)}}\rightarrow\mathbb{R}^{n^{(j)}}$ is called the
activation function and introduces nonlinearity into the output of the neuron
(Hastie et al.,, 2009), which may enable more complex pattern recognition and
more accurate predictions. If $g$ is the identity function, any NN is proved
to specialize into a linear regression model. The choice of the activation
function depends on the type of problem to be solved. Some common ones are (a)
the sigmoid function (or logistic function) (Han and Moraga,, 1995) that maps
the input to a range between 0 and 1, making it useful for binary
classification problems; (b) the rectified linear unit (ReLU) (Hahnloser et
al.,, 2000) that returns the input if it is positive, and zero otherwise and
is widely used due to its simplicity and computational efficiency; (c) the
hyperbolic tangent (Rumelhart et al.,, 1985) that squashes the input values to
the range between -1 and 1, making it useful for classification tasks; and (d)
the softmax function (Ackley et al.,, 1985) that is commonly used in the
output layer of a NN for multi-class classification problems, as it converts
the outputs into a probability distribution.

Traditional NNs accept only finite-dimensional vectors as input and thus, they
cannot easily handle profiles, whereas, the FNN introduced by Thind et al.,
(2023) is instead purposely designed for it.

Given $X_{1}(t),\ldots,X_{P}(t)$ and $z_{1},\ldots,z_{J}$ functional and
scalar covariates, respectively, Thind et al., (2023) introduce the functional
weights
$\boldsymbol{\gamma}(t)=\left(\gamma_{1},\ldots,\gamma_{p},\ldots,\gamma_{P}\right)^{\top}$
to effectively weigh the functional covariates at every point along their
domain $\mathcal{T}$, and define the output $h_{k}^{(1)}$ of the $k$th neuron
in the first hidden layer corresponding to the $i$th observation as follows

$h_{ik}^{(1)}=g\left(\sum_{p=1}^{P}\int_{\mathcal{T}}\gamma_{kp}(t)X_{ip}(t)dt+\sum_{j=1}^{J}w_{kj}^{(1)}z_{ij}+b_{k}^{(1)}\right)\quad
i=1,\ldots,n\quad k=1,\ldots,n^{(1)},$ (12)

where $n^{(1)}$ is the number of neurons in the first hidden layer and
$g(\cdot)$ is the activation function. This layer is referred to as a
functional hidden layer as it consists of neurons capable of handling infinite
dimensional functional weights $\boldsymbol{\gamma}(t)$. Each neuron in the
first layer produces a scalar value that is then fed into a regular NN. It is
worth noting that only the output of the first hidden layer has this
functional structure and this is the reason why $\gamma_{kp}(t)$ in Equation
(12) does not need for any superscript. This implies that the rest of the
hidden layers of the FNN can be of any of the usual forms (e.g., feedforward,
recurrent, residual) (LeCun et al.,, 2015). Similarly to the results in
Equation (9), the functional weights ${\gamma_{kp}}(t)$ in this first hidden
layer can be approximated through a linear combination of basis functions

$\hat{\gamma}_{kp}(t)=\sum_{m=1}^{M_{p}}c_{kpm}\zeta_{kpm}(t)=\boldsymbol{c}_{kp}^{T}\boldsymbol{\zeta}_{kp}(t)\quad
p=1,2,\ldots,P\quad k=1,2,\ldots,n^{(1)},$ (13)

where
$\boldsymbol{\zeta}_{kp}(t)=\left(\zeta_{kp1}(t),\ldots,\zeta_{kpM_{p}}(t)\right)^{T}$
is a vector of basis functions,
$\boldsymbol{c}_{kp}=\left(c_{kp1},\ldots,c_{kpM_{p}}\right)^{T}$ is the
corresponding vector of basis coefficients to be estimated by the NN, and
$M_{p}$ denotes the number of basis functions for each of the $P$ functional
covariates. Using the basis approximation in Equation (13), the general form
of the $k$th neuron in the first hidden layer in Equation (12) can be
rewritten as

$\displaystyle h_{k}^{(1)}$
$\displaystyle=g\left(\sum_{p=1}^{P}\int_{\mathcal{T}}\sum_{m=1}^{M_{P}}c_{kpm}\zeta_{kpm}(t)X_{p}(t)dt+\sum_{j=1}^{J}w_{kj}^{(1)}z_{j}+b_{k}^{(1)}\right)$
(14)
$\displaystyle=g\left(\sum_{p=1}^{P}\sum_{m=1}^{M_{P}}c_{kpm}\int_{\mathcal{T}}\zeta_{kpm}(t)X_{p}(t)dt+\sum_{j=1}^{J}w_{kj}^{(1)}z_{j}+b_{k}^{(1)}\right).$

In Equation (14), the integral can be approximated by using any numerical
integration method, e.g., the Simpson’s rule (Süli and Mayers,, 2003). To
compute the integral, the functional covariate $X_{p}(t)$ can be replaced by
$\hat{X}_{p}(t)$ as in Equation (8). The basis coefficients $c_{kpm}$ are
differently initialized for each functional weight $\gamma_{kp}(t)$ using the
Xavier uniform distribution (Glorot and Bengio,, 2010), and then will be
updated as the FNN learns together with the weights of the other non-
functional layers of the FNN. However, any other weight initialization methods
can be used, e.g., see Rastrigin, (1963); Kim and Ra, (1991); He et al.,
(2015). To train and assess the generalization performance of the FNN, the
mean squared error (MSE) $\sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}$ is
used as the loss function, where $y_{i},i=1,\ldots,N$, is the true scalar
response and $\hat{y}_{i}$ the FNN fitted output. The FNN can be trained with
the usual backpropagation algorithm (Rumelhart et al.,, 1985) and the Adam
optimizer (Kingma and Ba,, 2014). It is worth emphasizing that the FNN often
requires fewer parameters compared to standard NNs (e.g., MLP, CNN, RNN) that
directly process raw data. For instance, let us consider a single functional
covariate, measured $C$ times along its domain. The first hidden layer of a
traditional NN has a number of parameters equal to $(C+1)\times n^{(1)}$,
whereas the FNN requires $(M_{1}+1)\times n^{(1)}$, where $M_{1}$ is typically
lower than $C$ to avoid functional weight to overfitting.

Before training a NN there are many parameters to be specified, sometimes
referred to as hyperparameters. Tuning the values of the hyperparameters plays
a critical role in the generalization of the NN model. Typical hyperparameters
are the learning rate, which determines the step size in the optimization
process and significantly affects the convergence speed and the possibility of
getting stuck in local minima or overshooting global optima; the batch size,
which determines the number of samples processed before updating the model’s
weights, affecting the trade-off between computational efficiency and
generalization accuracy; the activation function, which introduces non-
linearity feature and may affect the NN convergence and ability to handle
vanishing or exploding gradients (Rumelhart et al.,, 1986); the NN
architecture, which includes the number of layers and the number of neurons in
each layer and significantly influence the ability of the NN to capture
complex patterns in the data. The number and the type of the basis functions,
which approximate each functional weight, can be considered hyperparameters as
well. In order to optimize the generalization performance of the FNN, we use a
hyperparameter tuning approach combining 5-fold cross-validation (Hastie et
al.,, 2009) with a grid search, which explores all possible combinations of
hyperparameters within a grid containing a wide range of values for each
hyperparameter of interest. By combining 5-fold cross-validation and grid
search, we ensure comprehensive coverage of the hyperparameter space and
identify the optimal parameter values that yield the smallest 5-fold cross-
validated MSE, defined as $\sum_{b=1}^{5}\sum_{i\in
S_{b}}\left(\hat{y}_{i}-y_{i}\right)^{2}/N$, where $S_{b}$ is the $b$th set of
observations in the held-out fold, and $\hat{y}_{i}$, is the predicted value
for $y_{i}$ by the FNN trained on the rest of the $B-1$ folds. In the
subsequent analysis, we tune all the FNN hyperparameters to some degree and
then use the early stop strategy (Keskar and Socher,, 2017) to prevent
overfitting and improve the generalization performance of the final model.
Instead of training the model for a fixed number of epochs, the early stop
strategy consists of monitoring the MSE on the validation set, which is a
separate subset of unseen observations, and stops the training process when
the latter increases.

The problem of interpretability in NNs has become a significant concern in the
field of DL. While NNs have demonstrated remarkable performance across various
tasks, their inherent black-box nature has hindered their wider adoption and
acceptance, particularly in domains where interpretability and transparency
are paramount. The black box refers to the opacity of NNs, where the internal
decision-making processes are not easily understandable or explainable to
humans. Researchers and practitioners are striving to develop techniques and
methodologies that can shed light on the black box, enabling us to gain deeper
insights into how NNs arrive at their predictions or classifications. The
pursuit of interpretability, whose rigorous definition is still debated in
machine learning literature (Molnar,, 2020; Lepore et al., 2022a, ; Stevens
and De Smedt,, 2023), is crucial for ensuring transparency and enabling domain
experts to understand, validate, and trust the decisions made by NNs. To this
aim, FNNs guarantee hidden semantic interpretability that refers to the human
ability to understand hidden layers (Fan et al.,, 2021). Specifically, the
functional weights introduced in Equation (14) differ from the traditional
weights as they can be easily visualized over the continuum, helping the
interpretation of the relationship between the functional covariates and the
scalar response while preserving the autocorrelation structure associated with
the data. These functional weights coincide with those estimated in the linear
functional regression model reported in Equation (10). If the functional
hidden layer has more than one neuron, the average of the estimated functional
weights $\hat{\gamma}_{p}(t)=\sum_{k=1}^{n^{(1)}}\hat{\gamma}_{kp}(t)/n^{(1)}$
for each $P$ functional covariate can be considered.

###  2.3 The functional neural network control chart

We propose a control charting procedure, referred to as FNNCC, to monitor a
scalar quality characteristic adjusted by the effect, possibly nonlinear, of
one or more functional covariates, which relies on the following main steps:

  1. [i]

  2. 1.

the relationship between the scalar response and the functional covariates is
modeled through an FNN;

  3. 2.

the functional model is estimated through FNN hyperparameter tuning, objective
function, and optimization algorithm definition (Goodfellow et al.,, 2016, p.
294-316);

  4. 3.

the monitoring strategy of residuals $e_{i}=y_{i}-\hat{y}_{i}$,
$i=1,2,\ldots,n$ obtained from the FNN is defined. Residuals act as the scalar
quality characteristic of interest to indirectly monitor the stability of the
functional relationship between the scalar response $y_{i}$ and the
multivariate functional covariates $X_{ip}$. For conciseness of notation, we
will hereinafter denote by $\hat{y}_{i}$ the fitted value of $y_{i}$, even
when an objective function different from MSE (as in Equation (11)) is used.

A current observation of the scalar response variable $y_{i}$, given the
corresponding functional covariate vector
$\mathbf{X}_{i}(t)=(X_{i1},\ldots,X_{iP})$, is monitored by using the FRCC
approach, that is the FNNCC results in a univariate control chart based on the
FNN residual $e=y-\hat{y}$. This paper focuses on the prospective monitoring
of the residuals, referred to as Phase II. That is, a data set of observations
representative of the IC process performance, referred to as Phase I sample or
reference data set, is assumed to be available. It is worth noting that the
retrospective monitoring, referred to as Phase I is crucial to check the
stability of historical functional data and to obtain accurate estimates of
the unknown parameters used for Phase II monitoring (Zhang et al.,, 2015). The
reference data set is randomly split into three non-overlapping sets, referred
to as training, validation, and tuning sets. The first is used to design and
train the FNN; the second to implement the early stopping strategy (Goodfellow
et al.,, 2016), which allows stopping the training when the performance does
not improve and thus, to prevent overfitting and improve the generalization of
the NN; the third is used to estimate the upper and lower control limits
(CLs), respectively, as $\alpha/2$ and $(1-\alpha/2)$ empirical quantiles of
the sampling distribution of the FNN residuals estimated from the tuning set,
being $\alpha$ the type-I error rate (Qiu,, 2013). In Phase II, the residual
of a new observation ($\mathbf{X}_{new},y_{new}$) is calculated as

$e_{new}=y_{new}-\hat{y}_{new},$ (15)

where $\hat{y}_{new}$ is the value fitted by the FNN model identified in Phase
I, according to the FNN model estimated with choices made in steps (i) and
(ii). An alarm is issued if $e_{new}$ is larger than UCL or lower than LCL.
The implementation of the FNNCC is outlined in Figure 1.

Figure 1: Outline of the FNNCC approach.

If the FNNCC issues an alarm, a change in the relationship between the
response and the covariates may have occurred. This could be due to changes in
the regression coefficients associated with one or more covariates, or
potential causes outside the set of covariates included in the model may be
investigated (Shu et al.,, 2004).

##  3 Simulation study

The overall performance of the proposed FNNCC is compared with FRCC and with a
univariate Shewhart control chart, referred to as SCC, which monitors the
scalar response without considering any information on the functional
covariates, in terms of $ARL_{1}$ at a given $ARL_{0}$ (Qiu,, 2013) by means
of an extensive Monte Carlo simulation. Without loss of generality, the
compact domain $\mathcal{T}$ is assumed as $[0,1]$, and the number of
covariates $P$ is set equal to 1. Then, the scalar response observations are
generated from an IC process as in the following scenarios

$\displaystyle\text{Scenario A}\quad y^{*}$ $\displaystyle=$ $\displaystyle
G(\mathbf{X}_{i})+\epsilon^{*}=\alpha+\int_{\mathcal{T}}\beta(t)x(t)\mathrm{d}t+\epsilon^{*}$
(16) $\displaystyle\text{Scenario B}\quad y^{*}$ $\displaystyle=$
$\displaystyle
G(\mathbf{X}_{i})+\epsilon^{*}=\exp\left(\alpha+\int_{\mathcal{T}}\beta(t)x(t)\mathrm{d}t\right)+\epsilon^{*}$
(17) $\displaystyle\text{Scenario C}\quad y^{*}$ $\displaystyle=$
$\displaystyle
G(\mathbf{X}_{i})+\epsilon^{*}=\Bigg{|}\alpha+\int_{\mathcal{T}}\beta(t)x(t)\mathrm{d}t\Bigg{|}+\epsilon^{*}$
(18) $\displaystyle\text{Scenario D}\quad y^{*}$ $\displaystyle=$
$\displaystyle
G(\mathbf{X}_{i})+\epsilon^{*}=\log\left(\left|\alpha+\int_{\mathcal{T}}\beta(t)x(t)\mathrm{d}t\right|+u\right)+\epsilon^{*}$
(19) $\displaystyle\text{Scenario E}\quad y^{*}$ $\displaystyle=$
$\displaystyle
G(\mathbf{X}_{i})+\epsilon^{*}=\left(\alpha+\int_{\mathcal{T}}\beta(t)x(t)\mathrm{d}t\right)^{2}+\epsilon^{*},$
(20)

where the noise $\epsilon^{*}$ is sampled from the Gaussian distribution
$N(0,0.1)$. where in each scenario, apart from Scenario A, $G(\cdot)$, defined
in Equation (1), specializes to a nonlinear SOF mapping, and
$\epsilon_{i}^{*}$ is the error term. It is worth noting that in Scenario D,
$u$ is an arbitrary positive constant, which is set equal to 2 to avoid
numerical problems due to small values. $y_{i}$ and
$\alpha+\int_{\mathcal{T}}\beta(t)x_{i}(t)\mathrm{d}t$ are generated using the
functions simulate_mfd(), which is called by the wrapper function
sim_funcharts() from the funcharts package (Capezza et al., 2023a, ). The mean
and the variance function $\mu^{X}(t)$ and $v^{X}(t)$ of the functional
covariates are generated according to the following model

$f(z)=P(z)+r\sum_{i=1}^{I}h_{i}\left(z;m_{i},s_{i}\right),\quad z\in(0,1),$
(22)

with

$P(z)=az^{2}+bz+c,$ (23)

where $a,b$ and $c$ are real numbers, and the terms
$h_{i}\left(z;m_{i},s_{i}\right)$ are normal probability density functions
with mean $m_{i}$ and standard deviation $s_{i}$. The functional covariate
${X}(t)$ is characterized by the Bessel (Abramowitz et al.,, 1964) correlation
function and is evaluated at 150 equally spaced discrete points of the
functional domain $\mathcal{T}=[0,1]$. As the simulated functional covariate
data are observed at noisy discrete values, each functional observation is
obtained by Equation (8) with $M=30$ cubic B-splines estimated through the
spline smoothing approach. Then, the scalar response is generated through a
SOF linear model so that the determination coefficient $R^{2}$ defined by Yao
et al., (2005), which measures the proportion of the variance in the response
variable explained by the functional covariate in the model, is set equal to
0.97. The mean and the variance of $y$ are set to $\mu^{y}=0$ and $v^{y}=1$,
respectively.

The performance of the proposed control chart is studied under a shift in
$\mu^{y}$ only, and both in the functional covariate mean $\mu^{X}(t)$ and
$\mu^{y}$. In the latter case, to study the unwanted effect of the shift in
$\mu^{X}(t)$ on the FNNCC performance, a translation of the profile pattern is
generated by using the model defined in Equation (22) with $P(z)$ defined as

$P(z)=az^{2}+bz+(c+\delta),\quad z\in(0,1),$ (24)

where $\delta$ is a real number defining the translation magnitude and is set
to be equal to 0.5. In both scenarios, the mean shift in the scalar response
is obtained by adding, to the simulated values, a fixed quantity that defines
the mean shift size. That is, $\mu_{y}$ is shifted by as much as
$\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$, where
$s^{y^{*}}=\sqrt{(v^{y^{*}})^{2}}$ is the standard deviation of the
transformed scalar response.

For each simulated scenario, a set of 4000, 1000, and 10000 IC scalar
responses and functional covariates are randomly generated to form the
training, validation, and tuning sets, respectively, for the reason discussed
in Section 2.3. To evaluate the $ARL_{1}$, an additional set of 20000 OC
patterns is randomly generated. $ARL_{0}$ is set to 20, which corresponds to
$\alpha=0.05$. In all scenarios, based on 5-fold cross-validation with grid
search, we use in step (ii) a two-layer FNN with 8 neurons and with ReLu and
linear activation functions. Each functional coefficient is expressed as a
linear combination of 5 cubic B-spline (see Equation (13)). Then, the FNN is
trained in step (ii) using the backpropagation algorithm with the Adam
optimizer to minimize the MSE chosen as the objective function. For the FRCC,
regression coefficients are estimated using the training set, while the tuning
set is used to compute the CLs as the $\alpha/2$ and $1-\alpha/2$ of empirical
quantiles of the scalar response.

Figure 2: Estimated $ARL_{1}$ achieved by FNNCC, FRCC, and SCC for each
simulated scenario, as a function of the mean shift size of the scalar
response $\Delta\mu_{y}={0.5,1,1.5,2s^{y^{*}}}$.

Figure 3: Estimated $ARL_{1}$ achieved by FNNCC, FRCC, and SCC for each
simulated scenario, as a function of the mean shift size of the scalar
response $\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$ when the functional
covariate is subject to a mean shift

Figures 2 graphically represent the $ARL_{1}$ performance achieved by FNNCC,
FRCC, and SCC for each simulated scenario, as a function of the mean shift
size of the scalar response $\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$, when
no covariate mean shift is considered. As we would expect, the SCC shows the
worst performance for all the considered scenarios and shifts as it is not
capable of adjusting the monitoring of the scalar response by the effect of
the additional information provided by the functional covariates. In Scenario
A, the FRCC performs comparably to the FNNCC as the linear functional
regression model is able to capture the true linear relationship between the
scalar response and the functional covariate. In all other scenarios, the
FNNCC outperforms the FRCC for all the considered shifts. This confirms the
SOF regression model is not capable of modeling the true nonlinear
relationship and, thus, the FRCC performance is the same as the SCC. For
example, for Scenario C, the $ARL_{1}$ of the FNN, FRCC, and SCC are 3.46,
11.83, and 11.89, respectively, at $\Delta\mu_{y}=0.5s^{y^{*}}$. The gain in
efficiency decreases as the shift size increases. As an example, when the
relation between the response and the functional covariate is modeled by an
exponential function (Scenario B), the $ARL_{1}$ of the FNNCC and of the FRCC
are 9.01 and 18.06 at a small shift size $\Delta\mu_{y}=0.5s^{y^{*}}$,
respectively, whereas for a higher shift size, say
$\Delta\mu_{y}=1.5s^{y^{*}}$, their performance decrease to 1.02 and 1.84,
respectively. When the functional covariate is subject to a mean shift, Figure
3 displays the estimated $ARL_{1}$ of the three competing control charting
schemes for all the simulated scenarios, as a function of the size of the
response mean shift. This figure points out that shifts in the covariate mean
function may impact the $ARL_{1}$ of the FNNCC and FRCC. It is trivial to note
that the SCC performance is unaffected by a change in $\mu^{x}(t)$. Simulation
results, displayed in Figure 3, show that the FNNCC still results more
sensitive than the FRCC in detecting OC condition of the scalar response and
that the $ARL_{1}$ performance of the FNNC and FRCC generally increases (or at
least remains the same) in the presence of a covariate mean shift. These
results are consistent with Centofanti et al., (2021); Shu et al., (2004).

The simulation study clearly highlights the superiority of the proposed method
in dealing with nonlinearity. As the true nature of the relation between the
scalar response and the functional covariates is never known in real-world
processes, FNNCC has proven to be a more flexible strategy to be recommended
in situations where the influence of the functional covariates is not
necessarily linear.

Moreover, by means of an extensive Monte Carlo simulation, the FNNCC is
compared with other two NN-based monitoring strategies implemented by using an
MLP in place of the FNN in step (i) of the monitoring strategy described in
Section 2.3. While both methods exhibit similar performance, the FNNCC is the
preferred choice due to its interpretability. The functional coefficients of
the FNN enable visualizing the relationship between functional covariates and
the scalar response. In contrast, the weights and intercepts of MLPs are
challenging to interpret, and the existing literature on interpreting
conventional NNs is primarily limited to computer vision applications. The
details of the numerical analysis are reported in Appendix A.

##  4 Case study: monitoring of HVAC systems on modern passenger trains

The case study mentioned in the Introduction is presented to demonstrate the
applicability of the proposed control chart in real situations. In recent
years, European regulations have been established to set operational standards
for the thermal comfort of passenger rail coaches. These standards, such as
EN, (2006), were developed to meet the different operating requirements of
rail vehicles and to ensure a high-quality air environment for passengers. In
view of those standards, railway companies are increasingly installing sensing
systems to collect data from onboard HVAC systems.

HVAC systems regulate the indoor temperature of each train’s coach through a
combination of ventilation, heating, and cooling operations. Ventilation is
the process of replacing or exchanging indoor air with outdoor air to remove
harmful particles like dust, smoke, and bacteria. Heating and cooling, on the
other hand, provide warmth or cold inside the coaches. An HVAC system has
three main components: the compressor, the condenser, and the thermal
expansion valve. The compressor moves the refrigerant gas to the condenser,
where the gas changes into a liquid. Then, the liquid refrigerant moves
through the evaporator section, where it evaporates into a cold gas, absorbing
heat from the surrounding air and cooling down the coach interior. Eventually,
the thermal expansion valve converts the cold gas back into a liquid, and the
process repeats multiple times until the indoor temperature reaches the
desired level.

The five passenger trains object of this study have six coaches, which are
equipped with a dedicated HVAC system, A central unit is installed to control
the heating and cooling modes of each HVAC system based on temperature sensors
that measure changes in the outdoor ($T_{\text{out}}$) and indoor
($T_{\text{in}}$) temperature signals that are streamed to maintenance
engineers for monitoring purposes and potentially improving reliability and
maintenance programs. At each time instant, on each coach, the HVAC system
activates until $T_{\text{in}}$ matches the target temperature
($T_{\text{set}}$), which is automatically and independently set as a function
of $T_{\text{out}}$, to comply with the current regulation on passengers’
comfort. The raw measurements of these temperatures are contained in the HVAC
data set where, for confidentiality reasons, the train names, acquisition
year, and routes are omitted. Trains and coaches are then identified with a
number from 1 to 5 and 1 to 6, respectively. This case study focuses on one
specific route where temperature measurements are acquired at a regular grid
of points equally spaced by 30 seconds. Railway engineers confirmed that the
30 HVAC systems can be assumed to be equal and working under the same
operating conditions (e.g., outdoor temperature, and passenger load) on the
same route. Thus, the coach effect is not considered hereinafter. Raw
measurements are then grouped to form raw profiles referred to different
voyages, which are identified by a unique voyage number (VN). Based on
experts’ opinion, exceptional voyages that do not represent normal operating
conditions have been removed to define the Phase I sample that is thus formed
by 1853 voyages and, according to Section 2.3, is randomly split into 740,
186, and 927 voyages to form training, validation, and tuning sets,
respectively.

Temperature profiles are re-mapped as a function of the fraction of the total
distance traveled by train at each voyage, and, to avoid the modeling of
transitional regimes, due to the HVAC restart at the beginning of each voyage,
the first $25\%$ of the traveled distance of each profile is discarded. This
operation can be regarded as a landmark registration (Ramsay and Silverman,,
2005, p. 129-132) of the functional data set from the function-specific
temporal domain to the common domain $\mathcal{T}$, which, without loss of
generality, is set as $[0,1]$. To get smooth profiles, i.e., functional data
observations, we choose a B-spline basis system with 70 basis functions and
equally spaced knots estimated by solving a regularization problem with a
roughness penalty on the integrated-squared second derivative and smoothing
parameter chosen through a generalized cross-validation (Ramsay and
Silverman,, 2005, p. 97- 99).

The scalar quality characteristic of interest is the root mean square (RMS) of
the difference between measurements $T_{\text{in}}andT_{\text{set}}$ acquired
during each train voyage, referred to as DevTemp. That is, the scalar
response, DevTemp, is computed once the voyage is ended as engineers cannot
perform maintenance operations until the train has finished its voyage or
reached the terminal station. DevTemp is assumed to be influenced by
$T_{\text{out}}$ and the $T_{\text{set}}$ derivative (with respect to the
fraction of the total distance travel), hereinafter denoted by
$\dot{T}_{\text{set}}$. The former is included in the model to the extent of
accounting for the thermal load under which the system works, whereas the
latter, based on experts’ opinion, is able to account for the thermal inertia
of the process, which is the time the HVAC system needs to allow the
$T_{\text{in}}$ to match $T_{\text{set}}$. $T_{\text{set}}$ may indeed vary as
a function of $T_{\text{out}}$, which, in turn, may be subject to rapid
changes, e.g., in under and over-ground route segments.

For illustrative purposes, Figure 4 displays a random slice of 100
observations from the training sample of the $T_{\text{out}}$ and
$T_{\text{set}}$ derivative profiles.

Figure 4: 100 random samples of the two functional profiles from the HVAC
training set.

The FNN is trained using a backpropagation algorithm with Adam optimizer and
early stopping, with the same set of hyperparameters discussed in the
simulation study in Section 3. Specifically, we first use the tuning data set
to build the FNNCC and estimate the CLs with a type-I error rate
$\alpha=0.05$. 41 additional profiles from coach 4 of train 5, which are known
to contain a fault due to a diagnosed failure in one of the two HVAC
compressors, are used as Phase II observations. Figure 5 shows the Phase II
monitoring of the aforementioned voyages. The $x$-axis label is the VN, while
the FNN residuals, defined in Section (2.3), are reported on the $y$-axis.

Figure 5: Phase II FNNCC. Each point corresponds to a voyage and the values of
the corresponding residual are reported. Horizontal dashed lines are the CLs
and the point above the CLs denotes the OC observation.

Voyage 32 shows an overly large value of the monitoring statistics, thus
correctly signaling the OC state of the corresponding HVAC system, which was
promptly repaired by train maintenance service. Indeed, subsequent voyages
plot inside the CLs, further demonstrating the proposed control chart to
properly track both the IC and OC states in practice.

##  5 Conclusions

A novel profile monitoring charting scheme is proposed in this work and
referred to as functional neural network control chart (FNNCC). Based on a
deep learning (DL) model, the FNNCC can adjust the monitoring of a scalar
quality characteristic of interest for the nonlinear effect of the functional
covariates when are available. Specifically, the FNNCC relies on functional
neural network (FNN), which has recently appeared in the DL literature to
allow a neural network to learn possibly nonlinear relationships based on
covariates in the form of functional data, aka profiles. The residuals
obtained from the FNN are elaborated to build a functional regression control
chart (FRCC), which also appeared in the statistical literature, although
implemented only under the (functional) linearity assumption. The proposed
FNNCC is thus the first DL-based profile monitoring scheme that can
efficiently exploit additional information on functional covariates, in a
possibly nonlinear fashion.

An extensive Monte Carlo simulation is carried out to assess the performance
of the proposed FNNCC in identifying mean shifts in the quality characteristic
of interest, which is in a scalar form, in the presence or absence of
covariate mean shifts. Then, the FNNCC is compared with the FRCC and the
Shewhart control chart (SCC) built on the scalar response. The results show
that the FNNCC is far better than the two competitors when the relation
between the scalar response and the functional covariate is nonlinear.
Additionally, the FNNCC is compared with two NN-based control charting
strategies and, even though they show similar performance, the FNNCC is
preferred due to its interpretability of the functional coefficients. The
practical applicability of the proposed control chart is finally illustrated
through a case study in the monitoring of heating, ventilation and air
conditioning systems installed on board six coach passenger trains, where the
favorable performance of the proposed method in properly tracking the IC and
OC states of the process is shown in practice.

The integration of neural networks into functional data analysis remains an
interesting topic as it efficiently enables nonlinear learning when the
covariates are in the form of profiles. Future research can be addressed to
extend the FNNCC to different and more sophisticated monitoring statistics.

## Acknowledgements

The authors are extremely grateful to the Operation Service and Maintenance
Product Evolution Department of Hitachi Rail Italy S.p.A. and, in particular,
to engineers Giuseppe Giannini, Vincenzo Criscuolo, and Guido Cesaro for their
technological insights in the interpretation of results.

This work was supported by the MOST – Sustainable Mobility National Research
Center and received funding from the European Union Next-GenerationEU (PIANO
NAZIONALE DI RIPRESA E RESILIENZA (PNRR) – MISSIONE 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.

The computing resources and the related technical support used for this work
have been provided by CRESCO/ENEAGRID High-Performance Computing
infrastructure and its staff (Iannone et al.,, 2019). CRESCO/ENEAGRID High-
Performance Computing infrastructure is funded by ENEA, the Italian National
Agency for New Technologies, Energy and Sustainable Economic Development, and
by Italian and European research programs, see
http://www.cresco.enea.it/english for information.

##  Appendix A Additional simulation study

In the simulation study, the analyses were carried out by comparing the
performance of the FNNCC with that of the FRCC and SCC in terms of the
$ARL_{1}$. In this Appendix, additional simulations are run for the nonlinear
scenarios (Scenario B-E) at each shift type and size discussed in Section 3,
to compare the proposed control chart with other two NN-based control charting
schemes, implemented by using an MLP in place of the FNN at step (i) of the
monitoring strategy discussed in Section 2.3 to model the nonlinear
relationships between the scalar response and the functional covariate. In
particular, in the first control charting strategy, referred to as
RawdataMLPCC, the raw discrete values of the functional covariate are directly
fed into the MLP (Rossi et al.,, 2002), whereas the second competitor, namely
BsplineMLPCC, involves a pre-processing step of the functional covariates into
a vector of scores from its B-spline expansion (Rossi et al.,, 2005). The two
MLPs, with the same set of the FNN hyperparameters defined in the simulation
study in Section 3, are trained on the training set using the backpropagation
algorithm with the Adam optimizer, and the CLs of the two competing methods
are estimated on the tuning set. Figures A1-A4 show the estimated $ARL_{1}$
for the FNNCC, RawdataMLPCC, and BsplineMLPCC for the nonlinear scenarios from
B to E, respectively, both in the presence and in the absence of a mean
covariate shift. It is clear from these figures that the RawdataMLPCC shows
worse performance than the functional counterparts as it is not able to
account for the functional nature of the covariate. FNNCC and BsplineMLPCC
achieve similar performance for all OC scenarios. However, the implementation
of the FRCC is to be preferred thanks to its easier interpretability of the
functional coefficients.

Figure 6: Estimated $ARL_{1}$ achieved by FNNCC, RawdataMLPCC, and
BsplineMLPCC in Scenario B both in the presence (a) and absence (b) of a mean
covariate shift, as a function of the mean shift size of the scalar response
$\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$.

Figure 7: Estimated $ARL_{1}$ achieved by FNNCC, RawdataMLPCC, and
BsplineMLPCC in Scenario C both in the presence (a) and absence (b) of a mean
covariate shift, as a function of the mean shift size of the scalar response
$\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$.

Figure 8: Estimated $ARL_{1}$ achieved by FNNCC, RawdataMLPCC, and
BsplineMLPCC in Scenario D both in the presence (a) and absence (b) of a mean
covariate shift, as a function of the mean shift size of the scalar response
$\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$.

Figure 9: Estimated $ARL_{1}$ achieved by FNNCC, RawdataMLPCC, and
BsplineMLPCC in Scenario E both in the presence (a) and absence (b) of a mean
covariate shift, as a function of the mean shift size of the scalar response
$\Delta\mu_{y}=\\{0.5,1,1.5,2\\}s^{y^{*}}$.

## References

  * Abramowitz et al., (1964) Abramowitz, M., Stegun, I. A., et al. (1964).  Handbook of mathematical functions, volume 55.  Dover New York.
  * Ackley et al., (1985) Ackley, D. H., Hinton, G. E., and Sejnowski, T. J. (1985).  A learning algorithm for boltzmann machines.  Cognitive science, 9(1):147–169.
  * (3) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2022a).  Functional regression control chart for monitoring ship co 2 _2 emissions.  Quality and Reliability Engineering International, 38(3):1519–1537.
  * (4) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2023a).  funcharts: control charts for multivariate functional data in r.  Journal of Quality Technology, 0(0):1–18.
  * (5) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2023b).  Functional regression control charts with an application to ship fuel consumption monitoring.  Wiley StatsRef: Statistics Reference Online, pages 1–10.
  * (6) 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.
  * 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.
  * Chen et al., (2020) Chen, S., Yu, J., and Wang, S. (2020).  Monitoring of complex profiles based on deep stacked denoising autoencoders.  Computers & Industrial Engineering, 143:106402.
  * 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.
  * Conan-Guez and Rossi, (2002) Conan-Guez, B. and Rossi, F. (2002).  Multi-layer perceptrons for functional data analysis: a projection based approach.  In Artificial Neural Networks—ICANN 2002: International Conference Madrid, Spain, August 28–30, 2002 Proceedings 12, pages 667–672. Springer.
  * Cybenko, (1989) Cybenko, G. (1989).  Approximation by superpositions of a sigmoidal function.  Mathematics of control, signals and systems, 2(4):303–314.
  * EN, (2006) EN, B. (2006).  Uni-en 14750-1: Railway applications—air conditioning for urban and suburban rolling stock. part 1: Comfort parameters.  British standard. London: British Standards Institution.
  * Fan et al., (2021) Fan, F.-L., Xiong, J., Li, M., and Wang, G. (2021).  On interpretability of artificial neural networks: A survey.  IEEE Transactions on Radiation and Plasma Medical Sciences, 5(6):741–760.
  * Ferraty, (2006) Ferraty, F. (2006).  Nonparametric functional data analysis.  Springer.
  * Glorot and Bengio, (2010) Glorot, X. and Bengio, Y. (2010).  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings.
  * Goodfellow et al., (2016) Goodfellow, I., Bengio, Y., and Courville, A. (2016).  Deep learning. book in preparation for mit press.  URL¡ http://www. deeplearningbook. org, 1.
  * Hahnloser et al., (2000) Hahnloser, R. H., Sarpeshkar, R., Mahowald, M. A., Douglas, R. J., and Seung, H. S. (2000).  Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit.  nature, 405(6789):947–951.
  * Han and Moraga, (1995) Han, J. and Moraga, C. (1995).  The influence of the sigmoid function parameters on the speed of backpropagation learning.  In From Natural to Artificial Neural Computation: International Workshop on Artificial Neural Networks Malaga-Torremolinos, Spain, June 7–9, 1995 Proceedings 3, pages 195–201. Springer.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., Friedman, J. H., and Friedman, J. H. (2009).  The elements of statistical learning: data mining, inference, and prediction, volume 2.  Springer.
  * He et al., (2015) He, K., Zhang, X., Ren, S., and Sun, J. (2015).  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034.
  * Hochreiter and Schmidhuber, (1997) Hochreiter, S. and Schmidhuber, J. (1997).  Long short-term memory.  Neural computation, 9(8):1735–1780.
  * Hornik et al., (1989) Hornik, K., Stinchcombe, M., and White, H. (1989).  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366.
  * Howard et al., (2018) Howard, P., Apley, D. W., and Runger, G. (2018).  Identifying nonlinear variation patterns with deep autoencoders.  IISE Transactions, 50(12):1089–1103.
  * Hsing and Eubank, (2015) Hsing, T. and Eubank, R. (2015).  Theoretical foundations of functional data analysis, with an introduction to linear operators, volume 997.  John Wiley & Sons.
  * Iannone et al., (2019) Iannone, F., Ambrosino, F., Bracco, G., De Rosa, M., Funel, A., Guarnieri, G., Migliori, S., Palombi, F., Ponti, G., Santomauro, G., and Procacci, P. (2019).  Cresco enea hpc clusters: a working example of a multifabric gpfs spectrum scale layout.  In 2019 International Conference on High Performance Computing Simulation (HPCS), pages 1051–1052. IEEE.
  * James et al., (2013) James, G., Witten, D., Hastie, T., Tibshirani, R., et al. (2013).  An introduction to statistical learning, volume 112.  Springer.
  * 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.
  * Keskar and Socher, (2017) Keskar, N. S. and Socher, R. (2017).  Improving generalization performance by switching from adam to sgd.  arXiv preprint arXiv:1712.07628.
  * Kim and Ra, (1991) Kim, Y. and Ra, J. (1991).  Weight value initialization for improving training speed in the backpropagation network.  In [Proceedings] 1991 IEEE International Joint Conference on Neural Networks, pages 2396–2401. IEEE.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Kokoszka and Reimherr, (2017) Kokoszka, P. and Reimherr, M. (2017).  Introduction to functional data analysis.  CRC press.
  * LeCun et al., (2015) LeCun, Y., Bengio, Y., and Hinton, G. (2015).  Deep learning.  nature, 521(7553):436–444.
  * LeCun et al., (1998) LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998).  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324.
  * (38) Lepore, A., Palumbo, B., and Poggi, J.-M. (2022a).  Interpretability for Industry 4.0: Statistical and Machine Learning Approaches.  Springer Nature.
  * (39) Lepore, A., Palumbo, B., and Sposito, G. (2022b).  Neural network based control charting for multiple stream processes with an application to hvac systems in passenger railway vehicles.  Applied Stochastic Models in Business and Industry, 38(5):862–883.
  * Lin et al., (2013) Lin, M., Chen, Q., and Yan, S. (2013).  Network in network.  arXiv preprint arXiv:1312.4400.
  * Mandel, (1969) Mandel, B. (1969).  The regression control chart.  Journal of Quality Technology, 1(1):1–9.
  * Menafoglio and Secchi, (2017) Menafoglio, A. and Secchi, P. (2017).  Statistical analysis of complex and spatially dependent data: a review of object oriented spatial statistics.  European journal of operational research, 258(2):401–410.
  * Molnar, (2020) Molnar, C. (2020).  Interpretable machine learning.  Lulu. com.
  * Montgomery, (2020) Montgomery, D. C. (2020).  Introduction to statistical quality control.  John Wiley & Sons.
  * Morris, (2015) Morris, J. S. (2015).  Functional regression.  Annual Review of Statistics and Its Application, 2:321–359.
  * Pacella and Semeraro, (2011) Pacella, M. and Semeraro, Q. (2011).  Monitoring roundness profiles based on an unsupervised neural network algorithm.  Computers & Industrial Engineering, 60(4):677–689.
  * Qiu, (2013) Qiu, P. (2013).  Introduction to Statistical Process Control.  Chapman & Hall.
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Ramsay and Silverman, (2005) Ramsay, J. O. and Silverman, B. W. (2005).  Functional Data Analysis.  Wiley Online Library.
  * Rastrigin, (1963) Rastrigin, L. (1963).  The convergence of the random search method in the extremal control of a many parameter system.  Automaton & Remote Control, 24:1337–1342.
  * Reiss et al., (2017) Reiss, P. T., Goldsmith, J., Shang, H. L., and Ogden, R. T. (2017).  Methods for scalar-on-function regression.  International Statistical Review, 85(2):228–249.
  * Rossi and Conan-Guez, (2006) Rossi, F. and Conan-Guez, B. (2006).  Theoretical properties of projection based multilayer perceptrons with functional inputs.  Neural Processing Letters, 23(1):55–70.
  * Rossi et al., (2002) Rossi, F., Conan-Guez, B., and Fleuret, F. (2002).  Functional data analysis with multi layer perceptrons.  In Proceedings of the 2002 International Joint Conference on Neural Networks. IJCNN’02 (Cat. No. 02CH37290), volume 3, pages 2843–2848. IEEE.
  * Rossi et al., (2005) Rossi, F., Delannay, N., Conan-Guez, B., and Verleysen, M. (2005).  Representation of functional data in neural networks.  Neurocomputing, 64:183–210.
  * Rumelhart et al., (1985) Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1985).  Learning internal representations by error propagation.  Technical report, California Univ San Diego La Jolla Inst for Cognitive Science.
  * Rumelhart et al., (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986).  Learning representations by back-propagating errors.  nature, 323(6088):533–536.
  * Saghaei et al., (2013) Saghaei, A., Noorossana, R., and Amiri, A. (2013).  Statistical analysis of profile monitoring.  Wiley Online Library.
  * Sergin and Yan, (2021) Sergin, N. D. and Yan, H. (2021).  Toward a better monitoring statistic for profile monitoring via variational autoencoders.  Journal of Quality Technology, 53(5):454–473.
  * Shu et al., (2004) Shu, L., Tsung, F., and Tsui, K.-L. (2004).  Run-length performance of regression control charts with estimated parameters.  Journal of Quality Technology, 36(3):280–292.
  * Stevens and De Smedt, (2023) Stevens, A. and De Smedt, J. (2023).  Explainability in process outcome prediction: Guidelines to obtain interpretable and faithful models.  European Journal of Operational Research.
  * Stuart et al., (1996) Stuart, M., Mullins, E., and Drew, E. (1996).  Statistical quality control and improvement.  European journal of operational research, 88(2):203–214.
  * Süli and Mayers, (2003) Süli, E. and Mayers, D. F. (2003).  An introduction to numerical analysis.  Cambridge university press.
  * Thind et al., (2023) Thind, B., Multani, K., and Cao, J. (2023).  Deep learning with functional inputs.  Journal of Computational and Graphical Statistics, 32(1):171–180.
  * Wang et al., (2016) Wang, J.-L., Chiou, J.-M., and Müller, H.-G. (2016).  Functional data analysis.  Annual Review of Statistics and its application, 3:257–295.
  * (65) Wang, Q., Zheng, S., Farahat, A., Serita, S., and Gupta, C. (2019a).  Remaining useful life estimation using functional data analysis.  In 2019 ieee international conference on prognostics and health management (icphm), pages 1–8. IEEE.
  * (66) Wang, Q., Zheng, S., Farahat, A., Serita, S., Saeki, T., and Gupta, C. (2019b).  Multilayer perceptron for sparse functional data.  In 2019 International Joint Conference on Neural Networks (IJCNN), pages 1–10. IEEE.
  * Williams et al., (2007) Williams, J. D., Woodall, W. H., and Birch, J. B. (2007).  Statistical monitoring of nonlinear product and process quality profiles.  Quality and Reliability Engineering International, 23(8):925–941.
  * Woodall et al., (2004) Woodall, W. H., Spitzner, D. J., Montgomery, D. C., and Gupta, S. (2004).  Using control charts to monitor process and product quality profiles.  Journal of Quality Technology, 36(3):309–320.
  * Yao et al., (2005) Yao, F., Müller, H.-G., and Wang, J.-L. (2005).  Functional linear regression analysis for longitudinal data.
  * Yao et al., (2021) Yao, J., Mueller, J., and Wang, J.-L. (2021).  Deep learning for functional data analysis with adaptive basis layers.  In International Conference on Machine Learning, pages 11898–11908. PMLR.
  * Yeganeh et al., (2022) Yeganeh, A., Abbasi, S. A., Pourpanah, F., Shadman, A., Johannssen, A., and Chukhrova, N. (2022).  An ensemble neural network framework for improving the detection ability of a base control chart in non-parametric profile monitoring.  Expert Systems with Applications, 204:117572.
  * Yeganeh and Shadman, (2021) Yeganeh, A. and Shadman, A. (2021).  Monitoring linear profiles using artificial neural networks with run rules.  Expert Systems with Applications, 168:114237.
  * Zhang et al., (2015) Zhang, J., Ren, H., Yao, R., Zou, C., and Wang, Z. (2015).  Phase i analysis of multivariate profiles based on regression adjustment.  Computers & Industrial Engineering, 85:132–144.
  * Zou et al., (2007) Zou, C., Tsung, F., and Wang, Z. (2007).  Monitoring general linear profiles using multivariate exponentially weighted moving average schemes.  Technometrics, 49(4):395–408.