# Real-World Anomaly Detection by using Digital Twin Systems and Weakly-
Supervised Learning

Castellani Andrea, Schmitt Sebastian, and Stefano Squartini A. Castellani and
S. Schmitt are with the Honda Research Institute Europe GmbH, 63071 Offenbach,
Germany (e-mail<EMAIL_ADDRESS>sebastian.schmitt@honda-ri.de)S. Squartini is with the Department of
Information Engineering, Università Politecnica delle Marche, 60131 Ancona,
Italy (e-mail: s.squartini@univpm.it)Manuscript received XXXXXXX XX, 2020;
revised XXXXXXX XX, 2020. DOI 10.1109/TII.2020.3019788

###### Abstract

The continuously growing amount of monitored data in the Industry 4.0 context
requires strong and reliable anomaly detection techniques. The advancement of
Digital Twin technologies allows for realistic simulations of complex
machinery , therefore , it is ideally suited to generate synthetic datasets
for the use in anomaly detection approaches , when compared to actual
measurement data. In this paper, we present novel weakly-supervised approaches
to anomaly detection for industrial settings, which make use of the Digital
Twin to generate a training dataset which simulates the normal operation of
the machinery, along with a small set of labeled anomalous measurement from
the real machinery. In particular, we introduce a clustering-based approach,
called Cluster Centers (CC), and a neural architecture based on the Siamese
Autoencoders (SAE), which are tailored for weakly-supervised settings with
very few labeled data samples. The performance of the proposed methods is
compared against various state-of-the-art anomaly detection algorithms on an
application to a real-world dataset from a facility monitoring system, by
using a multitude of performance measures. Also, the influence of hyper-
parameters related to feature extraction and network architecture is
investigated. We find that the proposed SAE based solutions outperform state-
of-the-art anomaly detection approaches very robustly for many different
hyper-parameter settings on all performance measures.

###### Index Terms:

Anomaly Detection, Digital Twin Systems, Siamese Neural Networks, Weakly-
Supervised Learning, Industry 4.0.

##  I Introduction

In recent years, there has been a strong tendency to equip technical
machinery, ranging from single machines to complete buildings and
manufacturing plants with sensors to constantly monitor their operation,
especially in the context of Industry 4.0 strategies[1]. Generally, component
failures or complete system failures need to be avoided as this severely
impacts the functioning of the machines and leads to significant increase in
maintenance, overhaul and repair (MRO) costs. In many situations, precursors
of component failures can be observed in the time-series of measured sensor
data and predictive maintenance approaches try to use this to reduce MRO
downtime and cost[2]. An essential part of these approaches are robust and
reliable anomaly detection methods which work in real-world settings.

Anomaly detection [3] refers to the problem of finding patterns in data that
do not conform to expected or normal behavior. It is an active area of
research with a wide range of application areas, such as energy [4],
manufacturing [5], network sensors [6], health care and video surveillance
[7]. Anomaly detection techniques based on machine learning can be separated
into different types of approaches [8]: supervised approaches, where a
sufficiently large set of training samples with labelled data is available;
unsupervised approaches, where only the unlabelled measurement data is
available; and weakly-supervised approaches, where a large amount of
unlabelled data with a very small set of labeled data is available. The
distinction between these cases is not clear-cut, as there could be a
supervised situation where a large amount of labelled data is available, but
this data is only from the normal operation. In that case, the training set
does not reflect the true distribution of data in the real-world, and even
more importantly, the most important information carrying data, i.e. labelled
anomaly samples, are missing.

For a practical application in realistic circumstances, even if there are many
data available, only unlabelled data is usually available, since it requires a
tremendous effort by human experts to manually create a fully labeled dataset
with a large portion of the possible anomalous scenarios. This requires
unsupervised anomaly detection methods to be used. In those methods, the bulk
statistics of the data is learned and the prevailing features of the data are
considered normal. The decision whether a data sample is normal or not is then
based on the comparison to the learned statistics so that infrequent data
samples are more likely to be considered anomalous. This implies a very high
class imbalance for the normal and the abnormal classes. Tuning the
sensitivity of the anomaly detection algorithm is a rather difficult problem
in such a situation. Even a very small false positive rate leads to a very
large absolute amount of wrongly detected anomalies which renders many
approaches useless for practical applications. Additionally, for most complex
machinery, the definition of an anomaly is not so clear. There might exist
normal operation modes, which are very rare and therefore, from a statistical
viewpoint, anomalies. Also real-world machinery exhibits drift of data due to
de-calibration of sensors, wear and tear and degradation of the machinery. So
the performance of an anomaly detection which was trained at a specific
instance in time might degrade over time due to the drift in the data and
regular re-training might be necessary.

In recent years, simulation technology has advanced substantially and even
complex machinery can nowadays be simulated quite accurately. It is even
possible to incorporate input data from actual measurements and to mirror the
operational history of individual machines, which is known as Digital Twin
systems first introduced in 2003 [9] (see also [10]). Depending on the effort
spent to create the Digital Twin, it can be realistic simulation which
captures the qualitatively correct behavior of the machinery, or an almost
perfect digital copy whose output can be directly compared to the real-world
machinery. The former can be at least used to create a large dataset
containing data samples of normal operation conditions which can be utilized
in machine learning approaches. The latter can play a key role in the anomaly
detection problem [11, 12, 13], when it runs in parallel to the physical
system with the same input values and environmental conditions.

In this paper, we present novel weakly-supervised approaches to anomaly
detection which we apply to the data from a company facility monitoring
system. We generate a large dataset of normal operation data covering a
complete year of operation using a Digital Twin simulation of the system. This
enables the accurate unsupervised learning of the statistics of the normal
operation states statistics including the very rare but normal states. In
order to increase the sensitivity for anomalies, we additionally use a very
small dataset of labeled anomalous samples which is obtained from the real-
world measurements. This weakly-supervised approach has the advantage to
produce a substantially lower false positive rate making it suitable for usage
in the real-world. Additionally, the Siamese Networks can be retrained very
efficiently as soon as new labelled data is available, making it ideal for
refinement during usage as well as for adjusting to drift and re-
configurations of the monitored machinery during operation.

The main contributions of this work are: (i) novel approaches to anomaly
detection using data from a Digital Twin simulation for the normal operational
state of a machinery, (ii) a clustering-based algorithm capable to solve the
anomaly detection task in both unsupervised and weakly-supervised settings,
(iii) a Siamese Autoencoder architecture for weakly-supervised anomaly
detection where very few labelled training samples are needed to improve the
performance over unsupervised anomaly detection methods, (iv) a thorough
comparison of experimental results of the proposed methods and their
performance to state-of-the-art algorithms for anomaly detection.

The rest of the paper is organized as follows. Section II presents an overview
of the recent literature on anomaly detection on time-series and the usage of
Digital Twin. Section III shows the proposed approach and presents the
developed algorithms. Section IV presents the used Digital Twin framework.
Section V shows the comparative anomaly detection algorithm investigated in
this work and Section VI presents the dataset and describes the experimental
setup. Section VII discusses the results obtained in our experimentation.
Finally, Section VIII concludes the paper and presents future developments.

##  II Related Works

Anomaly detection can be approached in different ways depending on the kind of
data available and the requirements of the particular use case. Chandola et
al. [3] presented a survey of different anomaly detection techniques applied
to different domains and recently the problem has been addressed from the
perspective of Deep Learning by Chalapathy et al. [8]. One category of deep
anomaly detection methods is based on reducing the dimension of the input by
mapping it to a low dimensional manifold followed by reconstruction. These
networks are usually trained in unsupervised fashion with a dataset containing
only data that reflects the normal state of operation. The reconstruction
error is then used to distinguish anomalies with several different neural
architectures, for example Autoencoders (AE) [14], Long Short Memory Network
(LSTM) [Wu2019a], Adversarial Autoencoder (AAE) [15] or Generative Adversarial
Networks (GAN) [16]. When dealing with multivariate time-series data, Canizo
et al. [17] use a CNN-RNN capable to reach an AP score of 0.994 on a real-
world dataset, but the proposed architecture is trained with a fully
supervised approach. Other approaches to anomaly detection make use of data-
driven statistical methods [18, 4], unsupervised clustering [19] and density-
based clustering [5].

Tao et al. [10] presents a detailed survey of the actual state-of-the-art of
the Digital Twin in industry. Several applications are presented in different
areas of design, production, prognostic and health management. To the best of
our knowledge, there are just few references in literature to the usage of
Digital Twin in the context of anomaly detection [12, 13, 11] and all of those
are more focused to present the Digital Twin system rather than a anomaly
detection algorithm.

While there are several supervised and unsupervised approaches to anomaly
detection described in literature, only few weakly-supervised approaches have
been proposed. In [20] the authors use domain adversarial training to transfer
the knowledge learned from a normal dataset to another with few labeled
samples. The Siamese Networks approach was introduced for signature
verification in [21] in a fully supervised setting. A first weakly-supervised
approach using Siamese Networks was made by Koch et al. [22]. Regarding
anomaly detection, there are some applications of Siamese Autoencoders
including human fall [23] and medical application [24]. Both architectures
make use of Siamese Autoencoder to extract features from the data, but the
anomaly detection process is made by an additional classifier, a KNN in the
former and OCSVM in the latter. Utkin et al. [25] presents a variation of the
Siamese Autoencoder with application to robotics. To the best of our
knowledge, no previous attempts to approach the anomaly detection problem with
weakly-supervised Siamese Networks and the combination of Digital Twin and
real-world data have been proposed.

##  III Proposed Algorithms

Figure 1: Proposed pipeline for Anomaly Detection.

The overall pipeline of the proposed anomaly detection scheme is sketched in
Figure 1. It consists of a training phase of an anomaly detection algorithm
(dotted rectangle) with the combination of two datasets: normal operation
dataset $\mathcal{N}$, generated with the Digital Twin simulation, and a small
set $\mathcal{A}$ of labeled anomalous samples from the real-world measurement
system. During the regular operational phase (blue boxes), the physical values
from the machinery provide a real-world dataset $\mathcal{R}$, which is fed to
the trained anomaly detection algorithm. An Anomaly Score (AS) is calculated
for each data sample and it is considered an anomaly if its score exceeds a
certain threshold.

###  III-A Cluster Centers

We propose the Cluster Centers (CC)algorithm as a weakly supervised
classification approach to the anomaly detection problem. In a first step, the
Digital Twin dataset of normal samples $\mathcal{N}$ is processed with an
unsupervised clustering algorithm, which identifies distinct normal
operational modes of the industrial machinery. Here, we use a $k$-means
approach, but in general any clustering algorithm is suitable, as long as it
is possible to calculate the set of cluster centers,
$\mathcal{C}=\\{\bm{c}_{j}\>\mathrm{for}\>j=1,...,N_{C}\\}$.

The number of clusters $N_{C}$ and the clusters themselves do not need to
represent the real semantically correct operational modes of the machinery,
but rather organize the data into different clusters which cover the complete
statistical variations of data found during normal operation. For this reason,
the number of clusters should be considerably larger than the (guessed or
known) true number of operational states so the clusters can also cover larger
variations, transitions and switching behavior properly.

The cluster centers $\bm{c}_{j}$ from the clustering algorithms of only normal
data samples from the Digital Twin can already be used to predict the
operational clusters on the real-world data set $\mathcal{R}$. The
distribution of cluster centers reflects only normal operation, and only from
the Digital Twin. The expectation is that the normal operations data from the
real-world aligns with this clustering, while anomalous data samples are quite
distinct from normal operation data and cannot be assigned reasonably to these
cluster centers. Therefore, the cluster centers can provide a fully
unsupervised estimate for the AS, which is given by the distance ($d$) from
each data point ($\bm{x}_{r}$) to its nearest cluster center ($\bm{c}$). For
the $k$-means approach we use the Euclidean distance, but any other distance
measure could also be used.

The set of labeled anomalies, $\mathcal{A}$, is used to refine the AS of each
real-world sample $\bm{x}_{r}$ by adding a penalty term,

$AS(\bm{x}_{r})={\underbrace{\vphantom{\frac{1}{\min(d(\bm{x}_{r},\bm{x}_{a})+\zeta)}}\min_{j}(d(\bm{c}_{j},\bm{x}_{r}))}_{\text{Unsupervised
AS}}}+{\underbrace{\eta\,\frac{1}{\min_{a}(d(\bm{x}_{r},\bm{x}_{a}))+\zeta}}_{\text{Penalty
term}}},$ (1)

with the cluster centers $\bm{c}_{j}\in\mathcal{C}$, the anomalous data
samples $\bm{x}_{a}\in\mathcal{A}$, a hyper-parameter $\eta>0$ determining the
influence of the penalty and a regularization factor $\zeta>0$. This penalty
term increases the AS whenever a known anomaly is close to the data sample.
With this, the information about the known anomalies is used to refine the
information provided by the cluster centers, which represents only normal
operation, and allows for a better discrimination between normal an anomalous
data samples.

The proposed CC and a regular $k$-Nearest Neighbors (KNN)clustering approach
differ in the way the AS is calculated. We use the distance to the nearest
cluster center, whereas in a KNN approach the nearest data samples are used.
Also, the penalty term is not used in regular KNN approaches.

In the above algorithm more hyper-parameters could have been introduced, for
example, an exponent different from 1 for the penalty term in (1) which would
allow for asymmetric weighting for small and large distances. This was
purposely not done to keep the approach as simple as possible. Also, instead
of including the labeled anomaly samples via a penalty term after the
clustering, they could have been included from the beginning by using a
constrained clustering algorithm. This was not done in order to deal more
easily with a constantly growing and changing set of labeled anomalies, which
can be used in the AS without redoing the clustering.

###  III-B Siamese Autoencoder

The Siamese Networks consist of two identical networks with shared weights
which can be efficiently used to decide if a pair of input data samples comes
from the same distribution or not.

The structure of the proposed architecture is shown in Figure 2. It comprises
an encoder $E$ which maps the input data sample $\bm{x}\in\mathbb{R}^{N}$ to a
latent representation $\bm{h}(\bm{x})\in\mathbb{R}^{M}$, with $M<N$. The
decoder network $D$ takes the latent representation and reconstructs a data
sample in the original space, $\tilde{\bm{x}}\in\mathbb{R}^{N}$. We refer the
reconstruction of a given input sample, $\bm{x}$, as
$\tilde{\bm{x}}=D(E(\bm{x}))=D(\bm{h}(\bm{x}))$, and its the latent
representation as $\bm{h}(\bm{x})=E(\bm{x})$.

The main goal of the proposed network is to reconstruct normal data samples
with low error and to create a clear separation of the normal and the
anomalous data distribution in the latent space. This can be achieved by the
following behaviour of the network:

  1. 1.

Normal data samples should be reconstructed as good as possible:
$\bm{x}\,\simeq\,\tilde{\bm{x}}$ for $\bm{x}\>\in\>\mathcal{N}$.

  2. 2.

The distance between two latent representations of any pair of normal data
samples should be as small as possible:
$d(\bm{h}(\bm{x}),\bm{h}(\bm{x}^{\prime}))$ small for
$\bm{x},\bm{x}^{\prime}\>\in\>\mathcal{N}$.

  3. 3.

Anomalous data samples should not be reconstructed well, i.e. their
reconstruction error should be much larger than the reconstruction errors for
normal data: $d(\bm{x},\tilde{\bm{x}})$ large for $\bm{x}\>\in\>\mathcal{A}$.

  4. 4.

The distance between the latent representations of a normal and an anomalous
data sample should be large: $d(\bm{h}(\bm{x}),\bm{h}(\bm{x}^{\prime}))$ large
for $\bm{x}\>\in\>\mathcal{N}$ and $\bm{x}^{\prime}\>\in\>\mathcal{A}$.

Figure 2: Siamese Autoencoder structure with the three loss functions.

The desired behaviour can be achieved by implementing an appropriate training
procedure for the neural network. The network always evaluates a pair of data
samples, $(\bm{x}_{DT},\bm{x}_{S})$, where the first sample is always taken
from the normal data set, $\bm{x}_{DT}\>\in\>\mathcal{N}$. In our approach
normal data samples for training are always taken from the Digital Twin
dataset, which is why we include the subscript $DT$ for clarity. The second
sample can be either from the normal or the anomalous dataset,
$\bm{x}_{S}\>\in\>\mathcal{N}\cup\mathcal{A}$. The loss function for training
then has the following three contributions

$\mathcal{L}=\mathcal{L}_{REC}+\mathcal{L}_{CL}+\mathcal{L}_{PCL}$ (2)

where each contribution is calculated as sums over all pairs of data
$(\bm{x}_{DT_{i}},\bm{x}_{S_{i}})$ and are given by the following expressions:

  * •

Reconstruction Loss: the Mean Square Error (MSE)between the input and its
reconstruction for normal operation data, which is just the cost function of a
typical Autoencoder:

$\mathcal{L}_{REC}=\frac{1}{N}\sum_{\mathrm{pairs}\>i}\big{(}\tilde{\bm{x}}_{DT_{i}}-\bm{x}_{DT_{i}}\big{)}^{2}.$
(3)

  * •

Contrastive Loss: the Euclidean distance between the latent vectors of the
input pair with a local modification:

$\mathcal{L}_{CL}=\frac{1}{2N}\sum_{\mathrm{pairs}\>i}\Bigg{(}(1-Y_{i})\,d\left(E(\bm{x}_{DT_{i}}),E(\bm{x}_{S_{i}})\right)^{2}\\\
+Y_{i}\left\\{\max\big{(}0,m-d\left(E(\bm{x}_{DT_{i}}),E(\bm{x}_{S_{i}})\right)\big{)}\right\\}^{2}\Bigg{)},$
(4)

where $Y_{i}=0$ if $\bm{x}_{S_{i}}\in\mathcal{N}$ and $Y_{i}=1$ if
$\bm{x}_{S_{i}}\in\mathcal{A}$. This contribution minimizes the difference for
input from the same class ($Y=0$) while maximize the differences for input of
different class ($Y=1$). The parameter $m>0$ allows only samples whose
distance is less than the radius defined by $m$ to contribute to the loss
function in order to avoid the domination of individual samples.

  * •

Partial Contrastive Loss: enforces a large reconstruction error for anomalous
data sample, which can also be viewed as a gradient inversion in the
Autoencoder training process:

$\mathcal{L}_{PCL}=\frac{1}{2N}\sum_{\mathrm{pairs}\>i}Y_{i}\>\max\big{(}0,m-d\left(E(\bm{x}_{DT_{i}}),E(\bm{x}_{S_{i}})\right)\big{)}$
(5)

Due to the possibility to construct a very large training dataset of distinct
pairs from the large normal dataset and very small anomalous dataset, this
approach is able to deal with extremely unbalanced datasets,
$|\mathcal{N}|\gg|\mathcal{A}|$.

After training, the AS for a new real-world data sample
$\bm{x}_{r}\in\mathcal{R}$ is calculated as:

$AS(\bm{x}_{r})\\!=\\!{\underbrace{\vphantom{\frac{1}{N^{\prime}}\sum_{i=1}^{N^{\prime}}}\left(\tilde{\bm{x}}_{r}-\bm{x}_{r}\right)^{2}}_{\text{Reconstruction
error}}}+{\underbrace{\frac{1}{N^{\prime}}\\!\sum_{i=1}^{N^{\prime}}\left\lVert
E(\bm{x}_{DT_{i}})-E(\bm{x}_{r})\right\rVert_{2}}_{\text{Embedding
distance}}},$ (6)

where $\bm{x}_{DT_{i}}$ are from a subset of the normal operation dataset from
the training phase with $N^{\prime}\leq N$ elements.

In this work we employ the Siamese architecture for two types of input data.
In the so-called Siamese Autoencoder (SAE)we extract feature vectors from the
measurement time series data and use standard feed forward autoencoder
architectures for encoder and decoder which are chosen to be symmetric. The
second variant called Convolutional Siamese Autoencoder (CNN-SAE)directly
operates on raw time series data and employs 1-D Convolutions, where encoder
and decoder are also chosen to be symmetric.

##  IV Digital Twin

The Digital Twin model used in this work simulates the electrical power
system, the Heating, Ventilation and Air Conditioning (HVAC), Combined Heat
and Power (CHP) systems for a medium sized company facility. The simulation is
realized with the Green City library111http://ea-energie.de/en/products/green-
city-simulationsbibliothek-2-2/ and the SimulationX software based on Modelica
programming language.The detailed process of the calibration of the Digital
Twin with the measurement data from the machinery is presented in a separate
publication [26]. This work focuses on the CHP model which takes the actual
measurements form weather data and total power demand time-series as input.

Due the high complexity of modelling realistic machines and the heating and
cooling of a company facility, the Digital Twin is only able to simulate with
very high fidelity power and energy consumption of the machines in a
reasonably large window of time, i.e. not less than few hours. The exact
transient behavior of the machines cannot be simulated exactly so comparing
raw time series between Digital Twin and the measurement data is not suitable
on the sampling rate of one data point per minute. However, we have validated
the Digital Twin model with several periods of correct operation of the
machine in 2018 and there is a discrepancy with the real-world electrical and
thermal energy production of $2.18\%$.

##  V Comparative Methods

The simplest method for deriving an AS is given by directly comparing the raw
time-series of the measurement with the simulation and taking the Mean
Absolute Error (MAE). As state-of-the-art algorithms for anomaly detection we
employ an Isolation Forest (IF), which builds an ensemble of isolation trees
for a given dataset, where anomalies are those data samples with short average
path lengths in the trees. In the $k$-Nearest Neighbors (KNN)approach the
measure of outlierness of an observation is obtained based on the distance to
its neighbors. We also compare the results of supervised Support Vector
Machine (SVM), its unsupervised variant One-Class SVM (OCSVM), as well as to
the Local Outlier Factor (LOF), which measures the outlierness of data points
by its local deviation of densities with respect to its neighbors. The
unsupervised dimensional reduction method Principal Component Analysis (PCA)is
used as an anomaly detection technique by taking the reconstruction error of
each sample as AS. Finally we also compare to a Deep Learning-based Feed
Forward Autoencoder (FF-AE)and a supervised Multi Layer Perceptron (MLP).

All unsupervised algorithms are trained on the Digital Twin dataset
$\mathcal{N}$ and then evaluated against the real-world data $\mathcal{R}$. In
addition, one OCSVM is also trained and tested on real-world data only. The
weakly-supervised and supervised algorithms are trained with the Digital Twin
dataset and few randomly sampled anomalies from $\mathcal{A}$.

##  VI Experimental Setup

###  VI-A Dataset Description

In this study, we use a dataset from a medium sized company, where an
infrastructure monitoring system recording various energy-related modules is
installed. The recorded data consists of a large number of sensors for heat,
cold and electricity consumption and production, as well as a local weather
station to monitor the ambient environment. This work focuses on data recorded
from the CHP module, where natural gas is burned to produce heat and
electrical power. The recorded time-series are: consumed and produced energy,
produced heat, heating fluid volume flow, and flow and return temperatures for
the heating fluid, as shown in Table I.

The dataset consists of data from November 2017 to February 2019. The sampling
rate is one data point per minute, leading to a total of 658081 sample points.

The synthetic Digital Twin dataset consists of the same time-series and the
same processing as for the real-world dataset, but only the year 2018 was
simulated.

The anomalous real-world dataset has been labeled manually. A total of 100
failure instances have been identified which translates to 24.2% of the total
number of samples being anomalous. There are two kinds of anomalies. One is
produced by sensor failures which result in a flat time series, and which are
rather easy to detect. More interesting and harder to detect are those
situations in which each time-series appears to provide a valid measurement,
but the combination multiple sensor reading does not reflect normal operation.
For the CHP, this is the case when, for example, the ambient temperature is
rather low but still no electricity and heat are produced. An example of a
normal and anomalous data sample is shown in Fig. 5.

###  VI-B Data Pre-Processing and Evaluation Metrics

TABLE I: Data Hyperparameters explored in the partial grid search. Values in bold have been used for the final configuration. Hyperparameter | Value
---|---
Variables |  Ambient Temperature ($T\text{a}$)
Produced Thermal Power ($P_{\text{Th}}$)
Produced Thermal Energy ($E_{\text{Th}}$)
Produced Electrical Power ($P_{\text{El}}$)
Produced Electrical Energy ($E_{\text{El}}$)
Water flux ($Flux$)
Water flow temperature ($T_{\text{Flow}}$)
Water return temperature ($T_{\text{Return}}$)
Flow/return temperature difference ($T_{\text{Diff}}$)
Statistical Features |  Mean ($\mu$)
Standard deviation ($\sigma$)
Skewness
Kurtosis
Sum
Root Mean Square
Time Features | True, False
Contextual Feature | True, False
Window Length | 240, 360, 480, 720, 1440
Window Stride | 30, 60, 120

For all approaches except the CNN-SAE, we use features derived from the raw
data as input to the models. The implemented statistical features are listed
in Table I and are extracted with a sliding window approach with parameters
also shown in the aforementioned table. The contextual feature refers to the
number of machine shutdowns (which can be easily detected from the data) and
the total working time within the time window. The time features include the
information related to working days and season of the year.

The feature vectors are composed of data from different domains and scales,
and therefore need to be standardized with the z-score
$:Z=(X-\bar{x})/\sigma$, where $\bar{x}$ is the mean of $X$ and $\sigma$ its
standard deviation. Regarding time and contextual values, One-Hot Encoding is
used to represent them as linearly uncorrelated vectors.

We evaluate the performance of the discussed approaches by using multiple
metrics. We use the $F_{2}$ Score, the Area Under Receiver Operating
Characteristic curve (AUC ROC) and the Average Precision (AP) as a measure of
the Area Under the Precision-Recall Curve (PRC). In the situation with highly
imbalanced class sizes, it was shown that the PRC is more informative than ROC
[27], since it better reflects the correct prediction of the minority class.

###  VI-C Implementation and Training Details

The state-of-the-art anomaly detection algorithms used for the comparison are
taken from the open source library PyOD [28]. The neural networks are
implemented with the Keras library and TensorFlow 1.13 back-end. The data
processing and training of the algorithms is done on two workstations, one
with a quad-core CPU (Intel Core i5-7300HQ) and 16 GB of DDR4 RAM and another
one with an Intel Xeon X5680 (6 cores and 12 threads), 128 GB of DDR4 RAM and
a NVIDIA TITAN X GPU.

TABLE II: Autoencoders Networks hyperparameters explored in the Grid Search.

Hyperparameter | Range | FF-AE | SAE | CNN-SAE
---|---|---|---|---
Batch size | $\\{16,32,64,128\\}$ | 64 | 128 | 32
Learning rate | $\left[10^{-2},10^{-5}\right]$ | $10^{-4}$ | $10^{-3}$ | $10^{-4}$
Dropout | $\left[0,0.5\right]$ | 0 | 0 | 0.1
Epoch | $\left[50,200\right]$ | 150 | 150 | 100
Activation function | $\\{\text{tanh},\text{ReLu},\text{sigmoid}\\}$ | tanh | ReLu | ReLu
Nr. of Encoder/Decoder layers | $\left[1,5\right]$ | 3 | 4 | 5
Encoder/Decoder neurons / kernel$*$ | $\left[8,512\right]$ | 64 32 16 | 128 64 32 16 | 32 64 128 128 256
Encoder/Decoder kernel length | $\left[1,10\right]$ | / | / | 10 7 5 4 3
Encoder/Decoder kernel stride | $\left[1,3\right]$ | / | / | 2 2 2 3 1
Nr. of fully-connected layers | $\left[1,3\right]$ | / | / | 1
Fully-connected neurons | $\left[10,100\right]$ | / | / | 60
Embedding space dimension | $\left[1,3\right]$ | 3 | 2 | 2

  * $*$

Only for Convolutional Networks

We use 10-fold Cross Validation, where the Digital Twin dataset and the subset
of labeled anomalies for training are divided in 10 folds to train and
validate the algorithm’s performance. The real-world dataset is divided in two
parts, $20\%$ for validation and $80\%$ for the final test evaluation.

In order to determine the best performing combination of hyper-parameters for
the feature extraction and the architectures of the Neural Networks, a Partial
Grid Search [29] has been used. Instead of searching in the complete hyper-
parameter space by computing all the possible combination of the parameters,
the search is organized as an iterated line-search. All parameters except one
are fixed and the minimization is performed only along the dimension of the
free parameter. Then, the best performing value for this parameter is selected
and the search is done along another parameter dimension. The resulting hyper-
parameters related to data and feature extraction are shown in Table I. The
values in bold constitute the best configuration of parameters, determined by
the search approach with Cross Validation. An example of this approach is
reported in Fig. 3, which illustrates the performance variations for some
different settings related to the feature selection. The left panel shows the
performance as function of the statistical quantities included in the feature
vector. The performance for different input time-series is instead shown in
the right panel. The qualitative changes and trends in the performance are the
same for all algorithms, only the SAE sometimes has a slightly different
trend. The variation of the statistical features (left panel) shows that the
best performance is achieved when only mean ($\mu$) and standard deviation
($\sigma$) are used. The inclusion of higher statistical moments like skewness
and kurtosis leads to a degradation of the performance most likely due to the
fact that the short time-scales, switching and state transitions, are not
captured well by the Digital Twin simulation and therefore cannot contribute
to a successful classification on the real-world data. By inspecting the curve
in the right panel we can state that the ambient temperature $T_{a}$ is a
rather important measurement, since not including it in the feature vector
leads to a performance drop from around 0.8 to around 0.3. This was expected
since the CHP is controlled on that variable. More sensors do not improve the
performance, also the inclusion of the flow and return temperatures do not
seem to be effective in this sense, for most algorithms. Only the SAE shows
high performance when the flow temperature is used.

Figure 3: Comparison of the AP performance for different settings for the
selected statistical features (left) and input time series (right). The
symbols and names shown on the $x$-axis are explained in Table I. Other
features are those shown in bold in Table I.

The resulting final dataset of feature vectors includes 8737 samples from the
Digital Twin and 10945 samples from the real-world, of which 2162 samples are
labeled anomalies. Each sample consists of a feature vector of dimension 6.

The hyper-parameters of the Neural Network architectures are shown in Table
II, with ranges used in the search approach reported in the first column. We
use a symmetrical structure for each Autoencoder. In particular, in the CNN
case, the encoder is composed of five convolutional layers followed by a
global max-pooling layer and then a fully-connected layer. The decoder has
five de-convolutional layers (up-sampling layer followed by a convolutional
layer). In order to revert the loss of information in the global max-pooling
layer, a skip-connection with linear activation is used between the last
convolutional layer of the encoder and the first de-convolutional layer of the
decoder. This is crucial to reconstruct the input data.

For training all neural architectures, we use the Adam optimizer with “LeCun
Uniform” weight initialization [30] and early stopping. We also apply gradient
clipping with maximum $L_{2}$ norm of 4 to improve training stability. For
weakly-supervised methods, the set of labeled anomalies for training is
composed of only 10 samples, randomly selected from the available labeled
data. For the CC algorithm, we set $\zeta=0.001$ and $\eta=0.15$, since those
are the best performing combination. For the OCSVM and SVC, we use a RBF
kernel with $C=1$. For the MLP we use a network with 3 hidden layer of 50
neurons each. Each experiment has been repeated 10 times and, as results, we
present the mean and standard deviation of the performance measures calculated
on the test set. The AS threshold, $\theta_{AS}$ is the value of AS in
correspondence of which the $25\%$ of the test-data has higher AS and it is
considered anomalous. This choice of $\theta_{AS}$ is based on an a-priori
assumption of the expected anomaly rate of the machinery under investigation.
Furthermore, the $\theta_{AS}$ needs to calculate confusion matrix of the
classifier, and then the $F_{2}$ score. In our experiments we use the AP score
to compare the performance of different algorithm and for the hyperparameter
search, because it is threshold-free and it is showed in [27] that it is the
best metric to use when dealing with an imbalanced classification task.

##  VII Experimental Results

TABLE III: Performance values of all evaluated algorithms , using the Digital Twin to generate the training-set. The globally best scores are highlighted in bold, while the best unsupervised scores are in italic font. For the weakly supervised methods only 10 additional labeled anomalies were used during training. Training Approach | Algorithm | Train Time | AP | AUC ROC | $\bm{F_{2}}$
---|---|---|---|---|---
Unsupervised |  Mean Absolute Error (MAE) | 0.41 s | $0.692\pm 0.002$ | $0.792\pm 0.001$ | $0.512\pm 0.001$
Isolation Forest (IF) | 0.54 s | $0.603\pm 0.024$ | $0.888\pm 0.011$ | $0.697\pm 0.026$
One-Class Support Vector Machine (OCSVM) | 3.54 s | $0.778\pm 0.001$ | $0.843\pm 0.001$ | $0.738\pm 0.001$
$k$-Nearest Neighbors (KNN) | 0.56 s | $0.807\pm 0.008$ | $0.893\pm 0.002$ | $0.720\pm 0.002$
Principal Component Analysis (PCA) | 0.09 s | $0.519\pm 0.001$ | $0.649\pm 0.001$ | $0.478\pm 0.001$
Local Outlier Factor (LOF) | 0.18 s | $0.811\pm 0.002$ | $0.869\pm 0.001$ | $0.711\pm 0.002$
Feed Forward Autoencoder (FF-AE) | 32.4 s | $\mathit{0.842\pm 0.007}$ | $\mathit{0.914\pm 0.006}$ | $\mathit{0.746\pm 0.013}$
Cluster Centers (CC) | 0.56 s | $0.831\pm 0.001$ | $0.909\pm 0.001$ | $0.734\pm 0.001$
Weakly-supervised |  Support Vector Classifier (SVC) | 7.8 s | $0.798\pm 0.028$ | $0.864\pm 0.024$ | $0.490\pm 0.040$
Multi Layer Perceptron (MLP) | 12.7 s | $0.732\pm 0.051$ | $0.813\pm 0.054$ | $0.017\pm 0.051$
Cluster Centers (CC) | 0.83 s | $0.848\pm 0.005$ | $0.921\pm 0.003$ | $0.762\pm 0.005$
Siamese Autoencoder (SAE) | 15.8 s | $\bm{0.872\pm 0.015}$ | $\bm{0.935\pm 0.013}$ | $\bm{0.823\pm 0.017}$
Convolutional Siamese Autoencoder (CNN-SAE) | 132.6 s | $0.866\pm 0.016$ | $0.931\pm 0.014$ | $0.798\pm 0.022$

(a)  FF-AE.  (b)  Unsupervised CC.  (c)  Weakly-supervised CC.  (d)  SAE.  (e)
CNN-SAE.

Figure 4: Confusion matrices of the proposed algorithms trained with Digital
Twin data: FF-AE (a), unsupervised CC (b), weakly-supervised CC (c), SAE (d),
CNN-SAE (e). The numbers are the absolute detections while the numbers given
in parenthesis are the in-class percentages.

First, we report the performance of the proposed unsupervised and weakly-
supervised approaches against other state-of-the-art algorithm for anomaly
detection, presented in section V, applied to a real-world use-case data,
introduced in section VI, by using the Digital Twin to generate the simulation
train data. Also, for the weakly-supervised approaches, we analyze how the
labeled anomalies used in the training affects the performance. Then, for
completeness, we show the usefulness of the Digital Twin simulation data, by
reporting the performance of the anomaly detection algorithms by using only
the real-world data.

###  VII-A Results obtained with Digital Twin data simulation

The results obtained in our experiments, by using the Digital Twin to generate
the training set, are reported in Table III. For the weakly-supervised
methods, 10 randomly selected anomalous data samples were included in the
training.

The naive MAE approach of directly comparing Digital Twin data to real
measurement data gives a medium performance of $AP=0.69$. It is worth noting
that simple unsupervised clustering-based approaches like KNN and the proposed
CC already give rather good scores of about $AP=0.81$ and $AP=0.83$,
respectively. The best performing unsupervised method is the deep learning
based FF-AE with a slightly better $AP=0.84$ while the training time increases
substantially.

All weakly-supervised algorithms show better performance in all metrics than
the unsupervised methods, which is to be expected since additional valuable
information on actual anomalies is used during training. The best anomaly
detection algorithm in our studies is the proposed SAE which reaches an AP
score of 0.872, AUC ROC of 0.935 and $F_{2}$ of 0.823, but with a fairly high
variance among our experiments. The larger variance is a direct result of the
generation of the training data, since randomly selected samples of normal
data are paired with only 10 randomly selected samples from anomalous dataset
to produce the input pairs. This produces a large variance in the selected
pairs from run to run, which explains the variance in the performance of
Siamese approaches. The SAE outperforms the weakly-supervised CC by $2.8\%$,
SVC by $19.1\%$ and MLP by $9.2\%$, in terms of AP score, when trained with
the same number of labeled samples. The CNN-SAE also reaches performance
levels above the state-of-the-art methods and comparable to the SAE but
requires a much longer training time.

In Figure 4 we show the confusion matrices of the best performing unsupervised
algorithm FF-AE (4a) and the proposed methods: unsupervised CC (4b) and
weakly-supervised CC (4c), SAE (4d) and CNN-SAE (4e). The false positive rate
goes from 12% to less than 10% with the change of the training approach from
unsupervised to weakly-supervised. And it is is about 9% for the best
performing algorithm in our studies, SAE, which is rather large but just about
to be acceptable for the application to the real world.

Figure 5:  Input, top row, and reconstructed sample, bottom row, of CNN-SAE.
The data have been standardized to have zero mean and unitary standard
deviation.  Figure 6: Overall performances of our experiments based on varying
the number of labeled anomalies.

An example data from normal and anomaly class and their respective
reconstruction, made with the CNN-SAE, is shown in Figure 5. The upper row
shows a normal and anomalous data sequence which is given as input to the
model. The bottom row shows the reconstruction of the convolutional network.
As it can be seen, the initial targets are achieved, since the normal time
series is accurately reconstructed, while the anomalous sample produces an
uncorrelated output.

For the previous results, we used only 10 anomalous data samples in the
weakly-supervised setting. In order to get some insights into the dependency
on the number of labeled anomalous data samples, we trained the algorithms
with differently sized sets of anomalies,
$|\mathcal{A}|=\\{5,10,50,100,500,1000\\}$, which corresponds to $0.05\%$,
$0.11\%$, $0.57\%$, $1.14\%$, $5.7\%$ and $11.44\%$ of the normal training
samples. As it is expected, there is a general trend that the performance
increases with more labeled anomaly samples as observable in Figure 6.
Remarkably, the proposed Siamese approaches already give a rather high AP
score above 0.85 with just 5 labeled anomaly samples. For the largest number
of labeled anomalies, the proposed CC algorithm show the best performance.
This is understandable, as in the clustering approach with the penalty term,
accordingly to Eq. 1, the information of each anomaly is directly used.
However, this is only beneficial as long as there is no noise in the set of
labeled anomalies. In a real-world setting with noisy labels, the performance
might not improve as much with more labeled data samples. The proposed Siamese
approaches are always clearly above the state-of-the-art methods. The increase
with additional labeled anomaly samples is rather slow which we take as a hint
that the Siamese approaches can be operated rather robustly in a real-world
setting with noisy anomaly data.

###  VII-B Results obtained with only Real-World data

TABLE IV:  Performance comparison of the investigated algorithm trained only with real-world data. The globally best scores are in bold, while the best unsupervised scores are in italic font. The $\Delta_{\%}$ is with respect to the use of Digital Twin data.  Training Approach | Algorithm | AP | $\bm{\Delta_{AP,\%}}$ | AUC ROC | $\bm{\Delta_{ROC,\%}}$
---|---|---|---|---|---
Unsupervised |  Isolation Forest (IF) | $0.479\pm 0.008$ | $-20.6$ | $0.800\pm 0.009$ | $-9.9$
One-Class Support Vector Machine (OCSVM) | $0.526\pm 0.004$ | $-32.4$ | $0.775\pm 0.005$ | $-8.1$
$k$-Nearest Neighbors (KNN) | $0.260\pm 0.002$ | $-67.8$ | $0.491\pm 0.002$ | $-45.0$
Principal Component Analysis (PCA) | $0.507\pm 0.002$ | $-2.3$ | $0.614\pm 0.001$ | $-5.4$
Local Outlier Factor (LOF) | $0.295\pm 0.004$ | $-63.6$ | $0.496\pm 0.006$ | $-42.9$
Feed Forward Autoencoder (FF-AE) | $\mathit{0.612\pm 0.012}$ | $-27.3$ | $\mathit{0.822\pm 0.012}$ | $-10.1$
Cluster Centers (CC) | $0.468\pm 0.017$ | $-43.7$ | $0.686\pm 0.008$ | $-24.5$
Weakly-supervised |  Support Vector Classifier (SVC) | $0.554\pm 0.085$ | $-30.6$ | $0.737\pm 0.036$ | $-14.7$
Multi Layer Perceptron (MLP) | $0.711\pm 0.031$ | $-2.9$ | $0.767\pm 0.026$ | $-5.7$
Cluster Centers (CC) | $0.584\pm 0.033$ | $-31.1$ | $0.774\pm 0.023$ | $-16.0$
Siamese Autoencoder (SAE) | $\bm{0.756\pm 0.019}$ | $-13.3$ | $\bm{0.876\pm 0.012}$ | $-6.3$
Convolutional Siamese Autoencoder (CNN-SAE) | $0.589\pm 0.023$ | $-32.0$ | $0.789\pm 0.020$ | $-15.3$

For the sake of completeness, we show the usefulness of the Digital Twin
simulation data by reporting the results obtained by using only the real-world
unlabeled data to train the algorithms. For the weakly-supervised training
approach, 10 randomly selected labeled anomalous samples were used during the
training. Thus, the anomalous samples were removed from the train-set. The
same networks structure and hyperparameters set of the previous results, were
used in those analysis.

The results obtained by using only the real-world data are reported in Table
IV, in terms of AP, ROC and their difference in percentage with respect to the
use of Digital Twin data ($\Delta_{AP,\%}$ and $\Delta_{ROC,\%}$). They show
the same trend of those with the Digital Twin simulation data, the best
performing unsupervised algorithm is the FF-AE, and the best algorithm is the
weakly-supervised SAE. But the performance is lower for all the investigated
algorithm for both the metrics reported, AP and ROC. This is due the fact that
the algorithms are not trained with a dataset that reflects only the normal
operation modes of the machinery, but the training data is noisy and with some
anomalous samples.

##  VIII Final Remarks and Conclusion

In this work, we present novel approaches to multi-variate time-series anomaly
detection. We demonstrate the approaches at an application to real-world data
from a facility monitoring system of a medium-sized company. We focus on the
Combined Heat and Power (CHP) module and use a Digital Twin simulation of the
facility to generate normal operation data for training. A small set of
labeled anomalies from the actual monitoring system is also available to train
the proposed weakly-supervised algorithms. We propose a simple clustering-
based approach where the anomaly score includes a penalty term which accounts
for the labeled data samples. We also present two approaches realizing Siamese
Neural Networks architectures, one taking features derived from the time-
series as input and the other directly operating with raw time-series data.
These architectures implement autoenconder neural networks, the loss function
targets the perfect reconstruction of normal operation data, while anomalous
data samples should not be reconstructed well. In order to enhance the
discriminatory power of the networks, the latent representations of normal and
anomalous samples is forced to have large distances.

We evaluate many statistical feature types and parameters for the time-series
features and observe that simple mean and standard deviation of a time
interval of one day, with a stride of one hour, give the best results. We
compare the proposed approaches to many state-of-the-art approaches and
evaluate a multitude of performance measures for comparison.

We validate the usefulness of the Digital Twin to generate the normal
operation data for training, all the investigated algorithms degrade their
performance by a remarkable amount when just the real-world data is used in
the training phase.

All the proposed weakly-supervised algorithms showed better performance that
the state-of-the-art approaches according to all performance measures. The
overall best performing algorithm is the Siamese Autoencoder operating on
time-series features, while the Feed Forward Autoencoder exhibits best
performance of all unsupervised algorithms. When varying the number of labeled
anomalous samples in the training set, the performance values of the proposed
weakly-supervised approaches change mildly and still have very good
performance for only five anomalous samples. In contrast, the performance
strongly degrades for smaller number of samples in the training set for the
Multi Layer Perceptron and the Support Vector Classifier.

The false positive and false negative detection rates of the best performing
method, SAE, are around 9% and 12% in each class, which is acceptable but
still rather large for a real-world application. However, the current approach
constitutes a generic proof of concept for using Siamese Networks for anomaly
detection tasks in real-world settings. Multiple options are available to
improve the performance for a specific application. Including real-world
normal states into the labeled dataset should enhance the performance due to
improving the alignment and transfer from purely simulated Digital Twin data
to real measurement data. Instead of just using a fixed predefined threshold
for the Anomaly Score, which is determined in a way to get 25% anomalies, a
more appropriate way should be used where, for example, the threshold should
be learned from the available data and maybe even adjusted with feedback from
an expert.

In the present study, the set of labeled anomalies is assumed to contain no
noise, i.e. no mislabeled samples. This assumption is not valid in a real-
world environment, as labelling errors are to be expected, but also the
definition of what constitutes an anomaly is not clear and might change over
time. As a consequence, a human-in-the-loop approach with expert feedback is
highly desirable where the set of labeled anomalies is re-inspected and
possibly extended over time. The proposed algorithms are designed to easily
incorporate the expert-feedback therein. However, their evaluation is
intentionally left for future work.

## References

  * [1] D. Preuveneers and E. Ilie-Zudor, “The intelligent industry of the future: A survey on emerging trends, research challenges and opportunities in industry 4.0,” _Journal of Ambient Intelligence and Smart Environments_ , vol. 9, no. 3, pp. 287–298, 2017.
  * [2] T. P. Carvalho, F. A. A. M. N. Soares, R. Vita, R. da P. Francisco, J. P. Basto, and S. G. S. Alcalá, “A systematic literature review of machine learning methods applied to predictive maintenance,” _Computers & Industrial Engineering_, vol. 137, p. 106024, 2019.
  * [3] V. Chandola, A. Banerjee, and V. Kumar, “Survey of Anomaly Detection,” _ACM Computing Survey (CSUR)_ , vol. 41, no. 3, pp. 1–72, 2009.
  * [4] Y. Zhang, Z. Y. Dong, W. Kong, and K. Meng, “A composite anomaly detection system for data-driven power plant condition monitoring,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 7, pp. 4390–4402, 2020.
  * [5] M. H. Chang, C. Chen, D. Das, and M. Pecht, “Anomaly detection of light-emitting diodes using the similarity-based metric test,” _IEEE Transactions on Industrial Informatics_ , vol. 10, no. 3, pp. 1852–1863, 2014\.
  * [6] D. Wu, Z. Jiang, X. Xie, X. Wei, W. Yu, and R. Li, “LSTM Learning With Bayesian and Gaussian Processing for Anomaly Detection in Industrial IoT,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 8, pp. 5244–5253, 2020.
  * [7] R. Nawaratne, D. Alahakoon, D. De Silva, and X. Yu, “Spatiotemporal anomaly detection using deep learning for real-time video surveillance,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 1, pp. 393–402, 2020.
  * [8] R. Chalapathy and S. Chawla, “Deep learning for anomaly detection: A survey,” _CoRR_ , vol. abs/1901.03407, 2019. [Online]. Available: http://arxiv.org/abs/1901.03407
  * [9] M. Grieves, “Digital twin: Manufacturing excellence through virtual factory replication,” White paper, 2014.
  * [10] F. Tao, H. Zhang, A. Liu, and A. Y. Nee, “Digital Twin in Industry: State-of-the-Art,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 4, pp. 2405–2415, 2019.
  * [11] G. Sugumar and A. Mathur, “Assessment of a Method for Detecting Process Anomalies Using Digital-Twinning,” _Proceedings - 2019 15th European Dependable Computing Conference, EDCC 2019_ , pp. 119–126, 2019.
  * [12] E. C. Balta, D. M. Tilbury, and K. Barton, “A digital twin framework for performance monitoring and anomaly detection in fused deposition modeling,” _IEEE International Conference on Automation Science and Engineering_ , vol. 2019-Augus, no. 70, pp. 823–829, 2019.
  * [13] P. Pileggi, J. Verriet, J. Broekhuijsen, C. Van Leeuwen, W. Wijbrandi, and M. Konsman, “A digital twin for cyber-physical energy systems,” _7th Workshop on Modeling and Simulation of Cyber-Physical Energy Systems, MSCPES 2019 - Held as part of CPS Week, Proceedings_ , pp. 1–6, 2019.
  * [14] E. Principi, D. Rossetti, S. Squartini, and F. Piazza, “Unsupervised electric motor fault detection by using deep autoencoders,” _IEEE/CAA Journal of Automatica Sinica_ , vol. 6, no. 2, pp. 441–451, 2019.
  * [15] J. Wu, Z. Zhao, C. Sun, R. Yan, and X. Chen, “Fault-attention generative probabilistic adversarial autoencoder for machine anomaly detection,” _IEEE Transactions on Industrial Informatics_ , pp. 1–10, 2020\.
  * [16] W. Jiang, Y. Hong, B. Zhou, X. He, and C. Cheng, “A GAN-Based Anomaly Detection Approach for Imbalanced Industrial Time Series,” _IEEE Access_ , vol. 7, pp. 143 608–143 619, 2019.
  * [17] M. Canizo, I. Triguero, A. Conde, and E. Onieva, “Multi-head CNN-RNN for multi-time series anomaly detection: An industrial case study,” _Neurocomputing_ , vol. 363, pp. 246–260, 2019.
  * [18] B. Wang, Z. Li, Z. Dai, N. Lawrence, and X. Yan, “Data-driven Mode Identification and Unsupervised Fault Detection for Nonlinear Multimode Processes,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 6, pp. 3651–3661, 2019.
  * [19] P. Li, O. Niggemann, and B. Hammer, “A Geometric Approach to Clustering Based Anomaly Detection for Industrial Applications,” _IECON 2018 - 44th Annual Conference of the IEEE Industrial Electronics Society_ , vol. 1, pp. 5345–5352, 2019.
  * [20] X. Li, W. Zhang, Q. Ding, and X. Li, “Diagnosing rotating machines with weakly supervised data using deep transfer learning,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 3, pp. 1688–1697, 2020\.
  * [21] J. Bromley, I. Guyon, Y. LeCun, E. Säckinger, and R. Shah, “Signature verification using a” siamese” time delay neural network,” in _Advances in neural information processing systems_ , 1994, pp. 737–744.
  * [22] G. Koch, R. Zemel, and R. Salakhutdinov, “Siamese neural networks for one-shot image recognition,” in _ICML deep learning workshop_ , vol. 2, 2015.
  * [23] D. Droghini, S. Squartini, E. Principi, L. Gabrielli, and F. Piazza, “Audio metric learning by using siamese autoencoders for one-shot human fall detection,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , pp. 1–11, 2019.
  * [24] Z. Alaverdyan, J. Jung, R. Bouet, and C. Lartizien, “Regularized siamese neural network for unsupervised outlier detection on brain multiparametric magnetic resonance imaging: Application to epilepsy lesion screening,” _Medical Image Analysis_ , vol. 60, p. 101618, 2020.
  * [25] L. V. Utkin, Y. A. Zhuk, and V. S. Zaborovsky, “An anomalous behavior detection of a robot system by using a hierarchical Siamese neural network,” _Proceedings of 2017 20th IEEE International Conference on Soft Computing and Measurements, SCM 2017_ , pp. 630–634, 2017.
  * [26] T. Schwan, S. Schmitt, and A. Castellani, “Calibration of HVAC system models with monitoring data - Digital Twin meets measurement data,” in _ESI Forum in Deutschland_. ESI, 2019.
  * [27] T. Saito and M. Rehmsmeier, “The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets,” _PLoS ONE_ , vol. 10, no. 3, pp. 1–21, 2015.
  * [28] Y. Zhao, Z. Nasrullah, and Z. Li, “Pyod: A python toolbox for scalable outlier detection,” _Journal of Machine Learning Research_ , vol. 20, no. 96, pp. 1–7, 2019.
  * [29] J. Bergstra and Y. Bengio, “Random search for hyper-parameter optimization,” _Journal of Machine Learning Research_ , vol. 13, no. Feb, pp. 281–305, 2012\.
  * [30] Y. A. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller, “Efficient backprop,” in _Neural networks: Tricks of the trade_. Springer, 2012, pp. 9–48.

|  Castellani Andrea was born in Civitanova Marche, Italy, in 1993. He
received the B.Sc. and M.Sc. degree (cum laude) in Electronic Engineering
respectively in 2016 and 2019 from Universitá Politecnica delle Marche
(Italy). In 2020 he started his Ph.D. perdiod at Uni Bielefeld (Germany) in
collaboration with Honda Research Institute (Germany). His current research
interests include Machine Learning and time-series analysis.
---|---
|  Schmitt Sebastian Sebastian Schmitt was born in 1975 in Germany. He
received his PhD in theoretical solid state physics in 2008 from Technical
University of Darmstadt, Germany. After a post-doc period from 2009 to 2011 at
the Technical University of Dortmund, Germany, he joined the Honda Research
Institute Europe in 2012. In his role as Senior Scientist he evaluates recent
technologies in optimization, data analytics and machine learning for their
use in real-world engineering settings relevant to R&D departments within
Honda.
---|---
|  Squartini Stefano (IEEE Senior Member, IEEE CIS Member) was born in Ancona,
Italy, on March 1976. He got the Italian Laurea with honors in electronic
engineering from Polytechnic University of Marche, UnivPM, in 2002. He
obtained his PhD at the same university in 2005. He joined the Department of
Information Engineering as Assistant Professor in Circuit Theory in 2007. He
is now Full Professor at UnivPM since 2020. His current research interests are
in the area of computational intelligence and digital signal processing, with
special focus on audio processing and energy management. He is author and
coauthor of more than 210 international scientific papers. He is Associate
Editor of the IEEE Transactions on Neural Networks and Learning Systems, IEEE
Transactions on Cybernetics and IEEE Transactions on Emerging Topics in
Computational Intelligence. He joined the Organizing and the Technical Program
Committees of more than 80 International Conferences and Workshops.
---|---


# Three Subtyping Algorithms for Binary Session Types and their Complexity
Analyses

Thien Udomsrirungruang University of Oxford, Oxford, UK
<EMAIL_ADDRESS>University of Oxford, Oxford, UK Nobuko
Yoshida University of Oxford, Oxford, UK<EMAIL_ADDRESS>
###### Abstract

Session types are a type discipline for describing and specifying
communication behaviours of concurrent processes. _Session subtyping_ ,
firstly introduced by Gay and Hole [4], is widely used for enlarging
typability of session programs. This paper gives the complexity analysis of
three algorithms for subtyping of synchronous binary session types. First, we
analyse the complexity of the algorithm from the original paper, which is
based on an inductive tree search. We then introduce its optimised version,
which improves the complexity, but is still exponential against the size of
the two types. Finally, we propose a new quadratic algorithm based on a graph
search using the concept of $\mathcal{XYZW}$-simulation, recently introduced
by Silva et al. [13].

##  1 Introduction

Session types [6, 14] are a type discipline for describing and specifying
communication behaviours of concurrent processes. They stem from the
observation that many real-world processes communicate in a highly-structured
manner; these communications may include sending and receiving messages,
selection and branching from a set of labels, and recursion. These session
types allow programs to be validated for type safety using typechecking
algorithms.

With regards to the synchronous session type system, Gay and Hole [4] proposed
a subtyping scheme, and proved the soundness of the type system; Chen et al.
[3] later proved that the type system was complete (i.e. there is no larger
subtyping relation between session types that is sound). Thus this notion of
subtyping is of interest as it is most general. Our presentation of session
types will use this system, as described in Gay and Hole’s paper [4]. For
simplicity, we omit the channel type $\hskip
3.00003pt\widehat{}\left[S\right]$: they can be adopted into the algorithms
with little issue.

As an example, consider the following type, which represents an interface to a
server which is able to respond to pings, as well as terminate itself:

$\displaystyle T^{\mathsf{interface}}_{\\_}1$ $\displaystyle=\mu
X.\oplus\negmedspace\langle\mathsf{respond}\\!:\>?[\mathsf{end}];X,\mathsf{exit}\\!:\mathsf{end}\rangle\negmedspace$
(1)

Here, the channel offers two choices of label to the process using it: (1)
choosing $\mathsf{respond}$ then receiving a message of unit type
$\mathsf{end}$, or (2) choosing $\mathsf{exit}$ and terminating the channel.
After performing (1), the $\mu$-recursion means that the channel reverts to
its original state.

Next, consider the scenario where we upgrade the interface to the server such
that it can be replicated.

$\displaystyle T^{\mathsf{interface}}_{\\_}2$ $\displaystyle=\mu
X.\oplus\negmedspace\langle\mathsf{respond}\\!:\>?[\mathsf{end}];X,\mathsf{exit}\\!:\mathsf{end},\mathsf{replicate}\\!:\>?[X];X\rangle\negmedspace$
(2)

The interface offers another choice: (3) choosing $\mathsf{replicate}$ and
receiving a new channel with the same type $T^{\mathsf{interface}}_{\\_}2$. We
can immediately see that any program that expects a channel of type
$T^{\mathsf{interface}}_{\\_}1$ will also work on a channel of type
$T^{\mathsf{interface}}_{\\_}2$. This is because the functionality of the new
channel type is identical to the old type, as long as the $\mathsf{replicate}$
label is never chosen. Thus, it would make sense to accept
$T^{\mathsf{interface}}_{\\_}2$ in type checking. It turns out that
$T^{\mathsf{interface}}_{\\_}2$ is a subtype of
$T^{\mathsf{interface}}_{\\_}1$ (denoted
$T^{\mathsf{interface}}_{\\_}2\leq_{c}T^{\mathsf{interface}}_{\\_}1$), which
characterises this property.

This subtyping relation is easy to see with a syntactic approach: the shapes
of the types are almost identical, except for one additional branch in
$T^{\mathsf{interface}}_{\\_}2$; we could then argue inductively that the
subtyping relation holds. However, for some types, this is not possible:

$\displaystyle T^{\mathsf{interface}}_{\\_}3$ $\displaystyle=\mu
Y.\oplus\negmedspace\langle\mathsf{respond}\\!:\>?[\mathsf{end}];Y,\mathsf{exit}\\!:\mathsf{end},\mathsf{replicate}\\!:\>?[T^{\mathsf{interface}}_{\\_}1];Y\rangle\negmedspace$
(3)

In essence, $T^{\mathsf{interface}}_{\\_}3$ is an interface that can be
replicated to get copies, but those copies cannot themselves be replicated. We
have that
$T^{\mathsf{interface}}_{\\_}2\leq_{c}T^{\mathsf{interface}}_{\\_}3$.
Syntactic subtyping does not work because the types do not have the same
shapes (namely, the messages sent in the $\mathsf{replicate}$ branch are $X$
and $T^{\mathsf{interface}}_{\\_}1$ respectively). The difficulty of checking
subtyping comes from the ability to "branch off" into multiple subterms: in
the above example, checking the subterm $\>?[T^{\mathsf{interface}}_{\\_}1];Y$
would require us to check both $T^{\mathsf{interface}}_{\\_}1$ and $Y$. In
certain cases this can lead to exponential complexity in the standard
inductive algorithm given in [4] (see Example LABEL:ex:exponential_example).
Contributions. We first analyse the algorithm given in Gay and Hole’s original
paper on session type subtyping [4, §5.1], improving the bound given in [9,
§5.1]. Second, we propose its optimised version, in the style of [12, Fig.
21-4]. Both of these algorithms have worst-case exponential complexity, with
the second algorithm having better complexity than the first. Third, we
represent types as labeled transition systems, and formulate the subtyping
problem as checking an $\mathcal{XYZW}$-simulation, similarly to [13]; we then
give a quadratic algorithm for subtyping by checking the validity of the
simulation.

##  2 Preliminaries

We restate the definitions given in [4], which are used in this paper.

######  Definition 2.1 (Session Types).

Session types (denoted $S,S^{\prime},T,U,V,W,\dots$) are defined by the
following grammar:

$\displaystyle S::=$ $\displaystyle\>\mathsf{end}$ (inaction) $\displaystyle|$
$\displaystyle\>\>?[S_{\\_}1,\dots,S_{\\_}n];S$ (input) $\displaystyle|$
$\displaystyle\>\>![S_{\\_}1,\dots,S_{\\_}n];S$ (output) $\displaystyle|$
$\displaystyle\>\oplus\negmedspace\langle
l_{\\_}i:S_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq n}$ (selection)
$\displaystyle|$ $\displaystyle\>\&\negmedspace\langle
l_{\\_}i:S_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq n}$ (branching)
$\displaystyle|$ $\displaystyle\>\mu X.S$ (recursive type) $\displaystyle|$
$\displaystyle\>X$ (type variable)

$X,X_{\\_}i,Y,Z,\dots$ range over a countable set of type variables. We
require that all terms are contractive, i.e. $\mu X_{\\_}1.\mu
X_{\\_}2.\dots\mu X_{\\_}n.X_{\\_}1$ is not allowed as a subterm for any
$n\geq 1$. As shorthand, we use $\>?[\tilde{T}];S$ instead of
$\>?[T_{\\_}1,\dots,T_{\\_}n];S$ when there is no ambiguity on $n$, and
similarly for $\>![\tilde{T}];S$.

######  Definition 2.2 (Unfolding).

$\text{unfold}\negmedspace\left(\mu
X.T\right)=\text{unfold}\negmedspace\left(T[\mu X.T/X]\right)$, and
$\text{unfold}\negmedspace\left(T\right)=T$ otherwise.

Note that as all terms are contractive, this is well-defined.

######  Definition 2.3 (Coinductive subtyping).

A relation $\mathcal{R}$ is a subtyping relation if the following rules hold,
for all $T\mathcal{R}U$:

  * •

If $\text{unfold}\negmedspace\left(T\right)=\>?[\tilde{T}^{\prime}];S_{\\_}1$
then
$\text{unfold}\negmedspace\left(U\right)=\>?[\tilde{U}^{\prime}];S_{\\_}2$,
and $\tilde{T}^{\prime}\mathcal{R}\tilde{U}^{\prime}$ and
$S_{\\_}1\mathcal{R}S_{\\_}2$.

  * •

If $\text{unfold}\negmedspace\left(T\right)=\>![\tilde{T}^{\prime}];S_{\\_}1$
then
$\text{unfold}\negmedspace\left(U\right)=\>![\tilde{U}^{\prime}];S_{\\_}2$,
and $\tilde{U}^{\prime}\mathcal{R}\tilde{T}^{\prime}$ and
$S_{\\_}1\mathcal{R}S_{\\_}2$.

  * •

If $\text{unfold}\negmedspace\left(T\right)=\&\negmedspace\langle
l_{\\_}i:T_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq m}$ then
$\text{unfold}\negmedspace\left(U\right)=\&\negmedspace\langle
l_{\\_}i:U_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq n}$, and $m\leq n$ and
$\forall i\in\\{1,\dots,m\\}.\>T_{\\_}i\mathcal{R}U_{\\_}i$.

  * •

If $\text{unfold}\negmedspace\left(T\right)=\oplus\negmedspace\langle
l_{\\_}i:T_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq m}$ then
$\text{unfold}\negmedspace\left(U\right)=\oplus\negmedspace\langle
l_{\\_}i:U_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq n}$, and $n\leq m$ and
$\forall i\in\\{1,\dots,n\\}.\>T_{\\_}i\mathcal{R}U_{\\_}i$.

  * •

If $\text{unfold}\negmedspace\left(T\right)=\mathsf{end}$ then
$\text{unfold}\negmedspace\left(U\right)=\mathsf{end}$.

Subtyping $\leq_{c}$ is defined by $S\leq_{c}T$ if $(S,T)\in\mathcal{R}$ in
some type simulation $\mathcal{R}$. It immediately follows that subtyping is
the largest type simulation.

######  Definition 2.4 (Coinductive equality).

$T=_{\\_}cT^{\prime}$ if $T\leq_{\\_}cT^{\prime}$ and
$T^{\prime}\leq_{\\_}cT$.

######  Example 2.5 (Interface).

We can now formally prove the subtyping relations mentioned in Section 1. Let
us define
$\mathcal{R}=\\{(T^{\mathsf{interface}}_{\\_}2,T^{\mathsf{interface}}_{\\_}3),(\>?[\mathsf{end}];T^{\mathsf{interface}}_{\\_}2,\>?[\mathsf{end}];T^{\mathsf{interface}}_{\\_}3),(\mathsf{end},\mathsf{end}),(T^{\mathsf{interface}}_{\\_}2,T^{\mathsf{interface}}_{\\_}1),\break(\>?[T^{\mathsf{interface}}_{\\_}2];{T^{\mathsf{interface}}_{\\_}2},\>?[T^{\mathsf{interface}}_{\\_}1];,{T^{\mathsf{interface}}_{\\_}3}),(\>?[\mathsf{end}];{T^{\mathsf{interface}}_{\\_}2},\>?[\mathsf{end}];{T^{\mathsf{interface}}_{\\_}1})\\}$.
We verify that the rules in Definition 2.3 hold. Thus $T\mathcal{R}U$ implies
$T\leq_{c}U$. In particular,
$T^{\mathsf{interface}}_{\\_}2\leq_{c}T^{\mathsf{interface}}_{\\_}3$ and
$T^{\mathsf{interface}}_{\\_}2\leq_{c}T^{\mathsf{interface}}_{\\_}1$.

To reason about the running time of the subtyping algorithms, we define the
size of a session type.

######  Definition 2.6 (Size).

$|\mathsf{end}|=|X|=1$, $|\mu X.T|=|T|+1$;
$|\>![T_{\\_}1,\dots,T_{\\_}n];U|=|\>?[T_{\\_}1,\dots,T_{\\_}n];U|=\sum_{\\_}{i=1}^{n}|T_{\\_}i|+|U|+1$;
$|\oplus\negmedspace\langle
l_{\\_}1:T_{\\_}1,\dots,l_{\\_}n:T_{\\_}n\rangle\negmedspace|=|\&\negmedspace\langle
l_{\\_}1:T_{\\_}1,\dots,l_{\\_}n:T_{\\_}n\rangle\negmedspace|=\sum_{\\_}{i=1}^{n}|T_{\\_}i|+1$.

We will also need a notion of subterms. The following function is defined in
[4].

######  Definition 2.7 (Bottom-up subterms).

The set of bottom-up subterms of $T$ is defined inductively as:
$\text{Sub}(\mathsf{end})=\\{\mathsf{end}\\}$; $\text{Sub}(X)=\\{X\\}$;
$\text{Sub}(\mu X.T)=\\{\mu X.T\\}\cup\\{S[\mu X.T/X]\mid
S\in\text{Sub}(T)\\}$;
$\text{Sub}(\>![T_{\\_}1,\dots,T_{\\_}n];S)=\\{\>![T_{\\_}1,\dots,T_{\\_}n];S\\}\cup\text{Sub}(S)\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}(T_{\\_}i)$;
$\text{Sub}(\>?[T_{\\_}1,\dots,T_{\\_}n];S)=\\{\>?[T_{\\_}1,\dots,T_{\\_}n];S\\}\cup\text{Sub}(S)\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}(T_{\\_}i)$;
$\text{Sub}(\oplus\negmedspace\langle
l_{\\_}1:T_{\\_}1,\dots,l_{\\_}n:T_{\\_}n\rangle\negmedspace)=\\{\oplus\negmedspace\langle
l_{\\_}1:T_{\\_}1,\dots,l_{\\_}n:T_{\\_}n\rangle\negmedspace
S\\}\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}(T_{\\_}i)$; and
$\text{Sub}(\&\negmedspace\langle
l_{\\_}1:T_{\\_}1,\dots,l_{\\_}n:T_{\\_}n\rangle\negmedspace)=\\{\&\negmedspace\langle
l_{\\_}1:T_{\\_}1,\dots,l_{\\_}n:T_{\\_}n\rangle\negmedspace
S\\}\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}(T_{\\_}i)$.

Then we define $\text{Sub}(T,U)=\text{Sub}(T)\cup\text{Sub}(U)$.

Due the use of unfolds (Definition 2.2) in our proof rules, it will be easier
to reason with a definition of subterms that use an operation similar to
unfolding in the treatment of recursive types. Thus, we give the analogue of
[12, Definition 21.9.1] for session types. The full definition is Appendix
A.1.

######  Definition 2.8 (Top-down subterms).

The set of top-down subterms of $T$ is defined: $\text{Sub}_{\text{TD}}(\mu
X.T)=\\{\mu X.T\\}\cup\text{Sub}_{\text{TD}}(T[\mu X.T/X])$, with all other
rules from Definition 2.7 with Sub replaced by $\text{Sub}_{\text{TD}}$.

##  3 Related work

The first algorithm for subtyping of recursive function types dates back to
Amadio and Cardelli [2], where they give an exponential algorithm for
recursive function type subtyping. Their algorithm inductively checks for
$\alpha\rightarrow\beta\leq\gamma\rightarrow\delta$ by recursively checking
that $\alpha\leq\gamma$ and $\beta\leq\delta$, unwrapping $\mu$-recursions,
and keeping track of which pairs have been checked before.

Pierce [12, Chapter 21.12] surveys the development of quadratic subtyping
algorithms for these types. Notably, Kozen et al. [8] represent types as
automata, then does a linear-time check on the product automaton of two types
to check for subtyping, a method similar to our third algorithm
(§LABEL:section:alg_3).

The subtyping relation for synchronous session types was introduced by Gay and
Hole [4], in which they showed that subtyping is sound and decidable. Their
algorithm for subtyping (as presented here in §LABEL:section:alg_1), is
similar to Amadio and Cardelli’s, adapted for session types. Later, Chen et
al.[3] proved that this relation is precise: no strictly larger relation
respects type safety.

Lange and Yoshida [9] investigate subtyping for session types by converting
terms to a modal $\mu$-calculus formula which represents its subtypes, then
using a model checker to check for subtyping. A short analysis of two other
subtyping algorithms is also provided, along with an empirical evaluation. The
first is Gay and Hole’s algorithm [4], in which a doubly-exponential upper
bound is given, which we improve to a singly-exponential bound. The second is
an adaptation of Kozen’s algorithm [8], which is similar to our LTS-based
construction (§LABEL:section:alg_3) in that it builds a product automaton and
checks for reachability. A quadratic algorithm is given; however, the type
system provided in [9] does not allow sending sessions as messages; instead,
sorts are used as payloads. Our algorithm is adapted to remove this
restriction. A comparison of complexity bounds in [9] and this paper is in
Table 1.

Algorithm | [9] | This paper
---|---|---
Gay and Hole [4] | $O\negmedspace\left(n^{2^{n}}\right)$ | $O\negmedspace\left(n^{n^{3}}\right),\Omega\negmedspace\left((\sqrt{n})!\right)$
[4] with memoization | – | $2^{O\negmedspace\left(n^{2}\right)},\Omega\negmedspace\left((\sqrt{n})!\right)$
Coinductive subtyping | $O\negmedspace\left(n^{2}\right)$ | $O\negmedspace\left(n^{2}\right)$
Table 1: Comparison between upper and lower bounds for the worst-case
complexity given in [9] and this paper to check the subtyping relation
$T\leavevmode\nobreak\ \leq_{c}\leavevmode\nobreak\ U$ where $n=|T|+|U|$.

Subtyping for non-regular session types, in which there are infinitely many
subterms, are explored by Silva et al. [13], in which the concept of
$\mathcal{XYZW}$-simulation for session type subtyping is introduced; it is
also used in this paper. A sound algorithm is introduced to semi-decide the
subtyping problem for context-free session types, which is accompanied by an
empirical evaluation. In general, subtyping for context-free session types is
undecidable, as shown by Padovani [10].

We believe that the results from this paper could be applied to other session
typing schemes with little difficulty: for example, the subtyping relation for
local synchronous multiparty session types [7], which is also sound and
complete [5].

### Acknowledgements.

The authors thank PLACES’24 reviewers for their careful reading and comments.
The first author is supported by the Keble Association Grant. The second
author is supported by EPSRC EP/T006544/2, EP/Y005244/1, EP/K011715/1,
EP/K034413/1, EP/L00058X/1, EP/N027833/2, EP/T014709/2, EP/V000462/1,
EP/X015955/1 and Horizon EU TaRDIS 101093006.

## References

  * [1]
  * [2] Roberto Amadio & Luca Cardelli (1993): _Subtyping Recursive Types_.  ACM transactions on programming languages and systems 15(4), pp. 575–631, 10.1145/155183.155231.
  * [3] Tzu Chun Chen, Mariangiola Dezani-Ciancaglini, Alceste Scalas & Nobuko Yoshida (2017): _On the Preciseness of Subtyping in Session Types_.  Logical Methods in Computer Science 13(2), 10.23638/LMCS-13(2:12)2017.
  * [4] Simon Gay & Malcolm Hole (2005): _Subtyping for Session Types in the Pi Calculus_.  Acta Informatica 42(2-3), pp. 191–225, 10.1007/s00236-005-0177-z.
  * [5] Silvia Ghilezan, Svetlana Jakšić, Jovanka Pantović, Alceste Scalas & Nobuko Yoshida (2019): _Precise Subtyping for Synchronous Multiparty Sessions_.  Journal of Logical and Algebraic Methods in Programming 104, pp. 127–173, 10.1016/j.jlamp.2018.12.002.
  * [6] Kohei Honda, Vasco Thudichum Vasconcelos & Makoto Kubo (1998): _Language Primitives and Type Discipline for Structured Communication-Based Programming_.  In Chris Hankin, editor: Programming Languages and Systems - ESOP’98, 7th European Symposium on Programming, Held as Part of the European Joint Conferences on the Theory and Practice of Software, ETAPS’98, Lisbon, Portugal, March 28 - April 4, 1998, Proceedings, Lecture Notes in Computer Science 1381, Springer, pp. 122–138, 10.1007/BFB0053567.
  * [7] Kohei Honda, Nobuko Yoshida & Marco Carbone (2008): _Multiparty Asynchronous Session Types_.  In: Proceedings of the 35th Annual ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages, POPL ’08, Association for Computing Machinery, New York, NY, USA, pp. 273–284, 10.1145/1328438.1328472.
  * [8] Dexter Kozen, Jens Palsberg & Michael I. Schwartzbach (1993): _Efficient Recursive Subtyping_.  In: Proceedings of the 20th ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages, POPL ’93, Association for Computing Machinery, New York, NY, USA, pp. 419–428, 10.1145/158511.158700.
  * [9] Julien Lange & Nobuko Yoshida (2016): _Characteristic Formulae for Session Types_.  In Marsha Chechik & Jean-François Raskin, editors: Tools and Algorithms for the Construction and Analysis of Systems, Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, pp. 833–850, 10.1007/978-3-662-49674-952.
  * [10] Luca Padovani (2019): _Context-Free Session Type Inference_.  ACM Transactions on Programming Languages and Systems 41(2), pp. 9:1–9:37, 10.1145/3229062.
  * [11] Robert Paige & Robert E. Tarjan (1987): _Three Partition Refinement Algorithms_.  SIAM Journal on Computing 16(6), pp. 973–989, 10.1137/0216062.
  * [12] Benjamin C. Pierce (2002): _Types and Programming Languages_ , 1st edition.  The MIT Press.
  * [13] Gil Silva, Andreia Mordido & Vasco T. Vasconcelos (2023): _Subtyping Context-Free Session Types_.  279, Schloss Dagstuhl – Leibniz-Zentrum für Informatik, pp. 11:1–11:19, 10.4230/LIPICS.CONCUR.2023.11.
  * [14] Kaku Takeuchi, Kohei Honda & Makoto Kubo (1994): _An Interaction-based Language and Its Typing System_.  In Constantine Halatsis, Dimitris G. Maritsas, George Philokyprou & Sergios Theodoridis, editors: PARLE ’94: Parallel Architectures and Languages Europe, 6th International PARLE Conference, Athens, Greece, July 4-8, 1994, Proceedings, Lecture Notes in Computer Science 817, Springer, pp. 398–413, 10.1007/3-540-58184-7118.

##  Appendix A Appendix

######  Definition A.1 (Full version of Definition 2.8).

The set of top-down subterms of $T$ is defined:

$\displaystyle\text{Sub}_{\text{TD}}(\mathsf{end})$
$\displaystyle=\\{\mathsf{end}\\}$ $\displaystyle\text{Sub}_{\text{TD}}(X)$
$\displaystyle=\\{X\\}$ $\displaystyle\text{Sub}_{\text{TD}}(\mu X.T)$
$\displaystyle=\\{\mu X.T\\}\cup\text{Sub}_{\text{TD}}(T[\mu X.T/X])$
$\displaystyle\text{Sub}_{\text{TD}}(\>![T_{\\_}1,\dots,T_{\\_}n];S)$
$\displaystyle=\\{\>![T_{\\_}1,\dots,T_{\\_}n];S\\}\cup\text{Sub}_{\text{TD}}(S)\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}_{\text{TD}}(T_{\\_}i)$
$\displaystyle\text{Sub}_{\text{TD}}(\>?[T_{\\_}1,\dots,T_{\\_}n];S)$
$\displaystyle=\\{\>?[T_{\\_}1,\dots,T_{\\_}n];S\\}\cup\text{Sub}_{\text{TD}}(S)\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}_{\text{TD}}(T_{\\_}i)$
$\displaystyle\text{Sub}_{\text{TD}}(\oplus\negmedspace\langle
T_{\\_}1,\dots,T_{\\_}n\rangle\negmedspace)$
$\displaystyle=\\{\>![T_{\\_}1,\dots,T_{\\_}n];S\\}\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}_{\text{TD}}(T_{\\_}i)$
$\displaystyle\text{Sub}_{\text{TD}}(\&\negmedspace\langle
T_{\\_}1,\dots,T_{\\_}n\rangle\negmedspace)$
$\displaystyle=\\{\>?[T_{\\_}1,\dots,T_{\\_}n];S\\}\cup\bigcup_{\\_}{i=1}^{n}\text{Sub}_{\text{TD}}(T_{\\_}i)$

1:function SubtypeQuad($\Gamma,T,U$)

2: if $\Gamma=\text{false}$ then

3: return false

4: end if

5: if $T\leq U\in\Gamma$ then

6: return $\Gamma$

7: end if

8: $\Gamma\leftarrow\Gamma\cup\\{T\leq U\\}$

9: if $T=\mathsf{end}$ and $U=\mathsf{end}$ then

10: return $\Gamma$

11: else if $\text{unfold}\negmedspace\left(T\right)=\>?[\tilde{T}];V$ and
$\text{unfold}\negmedspace\left(U\right)=\>?[\tilde{U}];W$ then

12: for $(T_{\\_}i,U_{\\_}i)\leftarrow(\tilde{T},\tilde{U})$ do

13: $\Gamma\leftarrow\textsc{SubtypeQuad}(\Gamma,T_{\\_}i,U_{\\_}i)$

14: end for

15: return SubtypeQuad($\Gamma,V,W$)

16: else if $\text{unfold}\negmedspace\left(T\right)=\>![\tilde{T}];V$ and
$\text{unfold}\negmedspace\left(U\right)=\>![\tilde{U}];W$ then

17: for $(T_{\\_}i,U_{\\_}i)\leftarrow(\tilde{T},\tilde{U})$ do

18: $\Gamma\leftarrow\textsc{SubtypeQuad}(\Gamma,U_{\\_}i,T_{\\_}i)$

19: end for

20: return SubtypeQuad($\Gamma,V,W$)

21: else if $\text{unfold}\negmedspace\left(T\right)=\&\negmedspace\langle
l_{\\_}i:T_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq m}$ and
$\text{unfold}\negmedspace\left(U\right)=\&\negmedspace\langle
l_{\\_}i:U_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq n}$ and $m\leq n$ then

22: for $i\leftarrow 1..m$ do

23: $\Gamma\leftarrow\textsc{SubtypeQuad}(\Gamma,T_{\\_}i,U_{\\_}i)$

24: end for

25: return $\Gamma$

26: else if $\text{unfold}\negmedspace\left(T\right)=\oplus\negmedspace\langle
l_{\\_}i:T_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq n}$ and
$\text{unfold}\negmedspace\left(U\right)=\oplus\negmedspace\langle
l_{\\_}i:U_{\\_}i\rangle\negmedspace_{\\_}{1\leq i\leq m}$ and $m\leq n$ then

27: for $i\leftarrow 1..m$ do

28: $\Gamma\leftarrow\textsc{SubtypeQuad}(\Gamma,T_{\\_}i,U_{\\_}i)$

29: end for

30: return $\Gamma$

31: else

32: return false

33: end if

34:end function

Figure A.1: The algorithm in Theorem LABEL:thm:quadratic_binary_subtyping
written in procedural form.

######  Remark A.2.

The above results hold for the subtyping relation $\leq_{c}$. However, if
coinductive equality $=_{\\_}c$ is required, from Definition
LABEL:def:subtyping_lts it is a bisimulation and so can be checked in $O(n\log
n)$ time using a bisimulation checking algorithm[11].