# Two-phase Dual COPOD Method for Anomaly Detection in Industrial Control
System

Emmanuel Aboah Boateng, Student Member, IEEE, and J.W. Bruce, Senior Member,
IEEE The authors are with the Department of Electrical and Computer
Engineering, Tennessee Technological University, Cookeville, TN, USA.

###### Abstract

Critical infrastructures like water treatment facilities and power plants
depend on industrial control systems (ICS) for monitoring and control, making
them vulnerable to cyber attacks and system malfunctions. Traditional ICS
anomaly detection methods lack transparency and interpretability, which make
it difficult for practitioners to understand and trust the results. This paper
proposes a two-phase dual Copula-based Outlier Detection (COPOD) method that
addresses these challenges. The first phase removes unwanted outliers using an
empirical cumulative distribution algorithm, and the second phase develops two
parallel COPOD models based on the output data of phase 1. The method is based
on empirical distribution functions, parameter-free, and provides
interpretability by quantifying each feature’s contribution to an anomaly. The
method is also computationally and memory-efficient, suitable for low- and
high-dimensional datasets. Experimental results demonstrate superior
performance in terms of F1-score and recall on three open-source ICS datasets,
enabling real-time ICS anomaly detection.

###### Index Terms:

Anomaly detection, industrial control systems, machine learning, cyber-
physical systems.

##  I Introduction

Modern critical infrastructures like water treatment facilities, oil
refineries, power grids, and nuclear and thermal power plants all include
industrial control systems (ICS) and are used to control and monitor a
physical process. A system known as an ICS is created by combining
computational, and communication components, sensors, actuators, Programmable
Logic Controllers (PLCs), Human Machine Interfaces (HMIs), and Supervisory
Control and Data Acquisition (SCADA) systems are some of the devices and
subsystems that make up an ICS.

For a given ICS setup, the physical layer’s field devices, or sensors and
actuators, control and manage the underlying industrial process. The
distributed control layer’s PLCs get information about the process’s present
condition via sensors. Pumps, valves, generators, and circuit breakers are
examples of actuators that receive control actions from the PLCs and carry
them out. For the purpose of executing human-assisted control actions, other
devices, such as the SCADA and HMIs at the supervisory control layer, enable
communication between a plant operator and the PLCs [1]. Figure 1 shows a
conceptual view of a typical ICS.

Figure 1: Conceptual View of Industrial Control Systems (ICS).

Because of the deployment of ICSs in critical infrastructure, ICSs are
desirable targets for adversaries. Any of the ICS subsystems may be
compromised by a skilled adversary who would then be able to manipulate the
actuators or sensor readings over time until a nefarious objective is
achieved. ICSs may experience service interruptions or material losses as a
result of malicious cyberattacks, which could have a detrimental impact on the
quality of life in societies. A significant instance is the 2010 Stuxnet-based
cyberattack against the Iranian nuclear program [2]. The worm used pertinent
data from the ICS to harm the centrifuges inside the factory by repeatedly
altering their rotation speed after infecting the Programmable Logic
Controllers (PLCs). Another such is the cyberattack in late 2015 [3] that
targeted the electrical grid in Ukraine. Three local power distribution
businesses had their information systems compromised by the attackers using
the BlackEnergy malware.

Machine learning-based intrusion detection solutions for ICSs have been
developed in response to the rising frequency and complexity of these attacks
[4, 5, 6]. A common and effective method among the several strategies
suggested in the scientific literature is the one-class classification [4, 7].
In the training stage, one-class classification-based solutions create a model
that embodies the ICS’s typical behavior. The detection system employs that
model during the production stage to check whether the behavior of the live
system adheres to a required standard. Anomalies are deviations from the norm
that are typically identified using classification error criteria. These
algorithms may be particularly sensitive to unusual behavior, including zero-
day anomalies or cyberattacks and malfunctioning hardware or sensors that
might avoid detection.

Another key challenge of deploying machine learning techniques for ICS is the
computational demands of the machine learning detection model [8]. ICS
technology has been around for several decades; as a result, deploying complex
machine learning models with high computational requirements can be
challenging. However, existing approaches are expensive to compute. Existing
approaches performances are limited to low-dimensional data due to the curse
of dimensionality. The majority of existing methods require the selection and
tuning of hyperparameters. Sometimes the selection of the hyperparameters is
arbitrary and subjective, producing various results during different
experiments of the same ICS setup.

Model interpretability is a crucial challenge in ICS anomaly detection. The
majority of the existing literature makes use of complex machine learning
techniques [9, 4, 5] or an ensemble of different anomaly detection algorithms
[10, 11, 12] in an effort achieve high anomaly detection performance. However,
these complex anomaly detection techniques suffer from model interpretability.
In ICS setup, it is essential to understand the root cause of anomalies to
make informed decisions during ICS troubleshooting and during ICS forensic
analysis. Because ICS control critical infrastructure in our society, anomaly
detection models must be capable of providing reasons for flagging any system
operation as an anomaly.

False alarms are also a challenge in ICS anomaly detection. To preserve system
availability and operational continuity, every anomaly detection system must
lower the erroneous categorization of ICS routine procedures. A false alarm
could occur due to the selection of an improper machine learning approach for
ICS anomaly detection. For example, supervised learning has been used by the
vast majority of researchers for ICS anomaly detection [13, 6]. Because this
supervised anomaly detection method is unable to detect unidentified
cyberattacks, huge, up-to-date databases, including every possible
cyberattack, are required to train the classification model. Such requirements
are infeasible in practice. Unsupervised machine learning techniques have also
been used in ICS anomaly detection [4, 9] in an effort to account for unknown
anomalies. Unsupervised anomaly detection algorithms assume that all training
data are a true representation of the normal operation of the ICS because
anomalies are unknown ahead of time. The assumption may not hold for complex
ICS with several sensors and actuators. Sensor and actuator readings may be
noisy. If noisy data is used in training anomaly detection models, their
predictions may not be a true reflection of the normal ICS operation. There is
a need to remove the noise and any apparent outliers in the data before
training.

This work tackles the aforementioned challenges by proposing a novel
unsupervised Two-phase Dual (TPD) Copula-Based Outlier Detection (COPOD)
method for ICS anomaly detection. The proposed TPD COPOD is an anomaly
detection method that consists of two sequential modeling stages and a dual
(parallel) modeling stage. The first phase of the method makes use of
Empirical Cumulative Distribution Functions (ECOD) to remove any obvious noise
in the training dataset in order to substantiate the hypothesis that the
training dataset truly represents the normal operation of a given ICS dataset.
ECOD uses information about the data distribution of the training dataset to
identify anomalies based on the assumption that anomalies are mostly isolated
events located in the tails of a distribution [14]. The second phase consists
of a dual COPOD architecture that utilizes the normal process data (output) of
phase 1 to develop two COPOD models. The raw normal process data is used to
develop the first COPOD model. The second COPOD model is developed using
normalized process data (excluding discrete variables). COPOD algorithm uses
the normal dataset to formulate an empirical copula, and then uses the copula
to predict the tail probabilities of each observation in the dataset [15]. The
magnitude of the tail probabilities represents the degree to which an
observation is an outlier. A decision function strategy is also introduced for
assigning final anomaly scores for the second phase. The motivation for the
proposed architecture is to be able to exploit ICS data by examining the data
using two latent representations to extract useful information to minimize
overfitting and improve the model’s anomaly detection capabilities.

The contributions of this work can be summarized below:

  1. 1.

Propose first known two-phase dual COPOD method for ICS anomaly detection;

  2. 2.

Introduce an efficient and scalable anomaly detection method suitable for both
low and high dimensional ICS datasets; and capable of real-time ICS anomaly
detection; and

  3. 3.

Introduce a robust, parameter-free ICS anomaly detection method based on
empirical distribution functions. The deterministic nature of all stages of
the method leads to the mitigation of the challenges involved in
hyperparameter selection.

  4. 4.

Propose a highly interpretable anomaly detection method that quantifies each
feature’s contribution towards an ICS anomaly

The rest of this article is organized as follows. Section II presents an
overview of the related works, followed by Sections III, which discusses the
details of the experimental setup and various datasets used for developing the
proposed anomaly detection method. Section IV presents the proposed method,
and after that, section V presents the results and discussions. Finally,
Section VI presents conclusions and provides recommendations for future work.

##  II Related Work

Anomaly and intrusion detection in ICSs has been the subject of extensive
research. System state prediction is the first stage in physics-based
detection, as mentioned in [16]. For instance, some research employed linear
dynamical system modeling [17, 18] or autoregressive models [19]. The linear
dynamical system proposed in [17] is based on a secure state estimation
algorithm that involves Kalman filters for attack detection and state
estimation against sensor attacks in a noiseless dynamical system. The work in
[19] makes use of Autoregressive Moving Average models for meter management
data anomaly detection. The ARMA model in [19] is a linear process that
assumes that the meter management data is stationary and when the data
fluctuates, it does so uniformly around a particular time. Unfortunately, the
linearity of the modeled system, which is generally not satisfied in ICSs, is
one of the assumptions made by both techniques.

Invariant-based techniques [20, 21] and specification-based system modeling
[22] have been shown to be efficient anomaly detection approaches. Specificity
or the requirement that the solution is suited to the system and its operating
circumstances is one of the key downsides of invariant-based techniques and
specification-based systems.

Approaches based on statistical models in the context of anomaly detection
typically start by fitting probability distributions to data points. Using the
fitted models, the statistical models decide which points are outliers.
Parametric and nonparametric methods are two main classes into which these
techniques are typically divided. The main distinction is that parametric
methods take the data to come from a parametric distribution under the
assumption; hence fitting such a distribution entails learning the parameters
of the assumed parametric distribution. Both linear regression and Gaussian
mixture models (GMM) are used parametric techniques for anomaly detection [23,
24]. In contrast, nonparametric techniques do not rely on a parametric model
of the data. Examples include histogram-based techniques (HBOS) [25], Kernel
Density Estimation (KDE) [26]. After the model fitting process, parametric
models for anomaly detection are typically quick to utilize; however,
nonparametric models for anomaly detection can be more expensive to deal with.

To obtain reliable anomaly detection results, ensemble-based methods integrate
the output from different base outlier detectors. Notable works of ensemble
based anomaly detection methods include feature bagging [27], isolation forest
[28], locally selective methods [29], and scalable unsupervised outlier
detection [30]. Feature bagging method employs a variety of sub-feature
spaces. Isolation forests method combines data from multiple base trees.
Locally selective combination in parallel outlier ensembles method dynamically
selects the best base estimator for each data point. Scalable unsupervised
outlier detection method employs numerous heterogeneous estimators. In
general, ensemble-based approaches for anomaly detection frequently perform
well in practice, even for high-dimensional datasets [31, 32]. However,
ensemble-based methods can require non-trivial tweaking, such as choosing the
appropriate meta-detectors [33]. Ensemble anomaly detection methods are
frequently difficult to interpret [34].

The COPOD anomaly detection algorithm is modeled after copulas for
multivariate data distribution [15]. Copulas are mathematical functions that
let COPOD distinguish between the dependency structure of a given multivariate
distribution and the marginal distributions of a dataset. Empirical Cumulative
Distribution Functions (ECDF) are first calculated by the authors using a
specified dataset. The empirical copula functions are then created using the
(ECDF). Finally, the COPOD anomaly detection algorithm uses the empirical
copula to approximate the tail probabilities and quantifies the probabilities
as the anomaly scores of the data records. The authors in [15] evaluated their
proposed algorithm on 30 public outlier detection benchmark datasets. In each
instance, 60% of the data was used for training, and the remaining 40% of the
data was used for testing. Through extensive experiments conducted, the
authors claim the detection algorithm is a state-of-the-art outlier detection
algorithm in terms of detection accuracy and computational cost [15].

An unsupervised outlier detection algorithm known as ECOD, which is inspired
by the fact that outliers are often the rare events that occur at the tails of
a distribution, is proposed in [14]. The authors’ method models each dimension
in the dataset using a nonparametric statistical method. The proposed method
calculates the ECDF per dimension for each data point in order to first
estimate the underlying distribution of the input dataset. The algorithm then
estimates the dataset’s tail probabilities across all dimensions using the
obtained ECDF. Finally, quantification of the tail probability is utilized to
score the data records for outliers in the dataset. The authors conducted in-
depth tests using 30 benchmark datasets and claimed that the ECOD algorithm
outperforms 11 state-of-the-art baselines in terms of detection accuracy,
model efficiency, and scalability. Despite their model’s robust performance,
the work in [14] made a simplistic unimodal assumption about the input
dataset, and so the approach is not appropriate for multimodal distributions
in which an outlier could be in neither left nor right tails.

The detection of anomalies in PLCs has also been done using NN-based
prediction algorithms [35]. For the purpose of identifying PLC anomalies
employing long short-term memory (LSTM), the authors in [35] presented a non-
intrusive power-based anomaly detection approach. With the proposed models,
accuracy of up to 98% was obtained. However, identifying malicious code by
monitoring the power consumption of the PLC is insufficient as a faulty power
supply and a power electronics breakdown, can result in false positive
readings [36].

##  III Experiment Setup

This section provides the details of the threat model for this study and the
various experimental setups and datasets considered for developing and
evaluating the proposed anomaly detection method.

###  III-A Threat Model

In this work, it is assumed that an attacker is capable of physically
compromising sensors and actuators inside the ICS, gaining remote access to
the SCADA workstation, and gaining access to the networked control system.
This work further assumes that the attacker is familiar with the targeted ICS,
including the physical characteristics measured by each sensor and the effects
of actuation commands. The objective of the attacker is to harm or alter the
ICS activities using the aforementioned capabilities and prior information.
This includes; manipulating actuator states through an attack akin to Man-In-
the-Middle (MITM) attack in which the attacker sends orders to the actuator
rather than a PLC; delivering forged sensor values to the PLC to influence the
PLC to make poor judgments; and altering the PLC firmware with the intention
to disable the ICS or alter the programmed logic.

###  III-B Datasets Description

The proposed method is trained, evaluated, and compared with previous work
using three open source ICS datasets, namely Secure Water Treatment (SWaT)
[37], Water Distribution (WADI) [38], and the Traffic Light (TLIGHT) datasets
[4, 39]. This section provides a summary of the main properties of the
datasets utilized in this work.

The SWaT testbed is a typical representation of the dynamic nature of Cyber-
physical Systems (CPSs) used in our societies [40]. The dataset is subdivided
into a 7-day fraction of regular operations, which serve as the training set,
and a 4-day chunk of anomalous operations and 36 attacks produced using the
attack model in [38]. The WADI dataset represents a scaled-down version of a
large water distribution network in a city. The data records in the WADI
dataset, each with 123 attributes broken down into 69 sensor readings and 54
actuator states, were gathered over the course of 16 days. For the first 14
days, normal operating circumstances were recorded and divided into training
(95%) and validation (5%) sets. The test set consists of the final two days
containing 15 attacks interspersed with normal operating conditions. A total
of 784,568 data records were used for training the proposed model and 172,800
data records were used for testing the proposed model. The experimental setup
for emulating the behavior of the traffic light (TLIGHT) system for the
purposes of data collection is described in [4, 39]. The TLIGHT dataset
consisting of normal and abnormal data records, is recorded for a period of
four days. Overall, five test sets are presented in the TLIGHT dataset, with
each set containing normal and anomalous TLIGHT system operations.

##  IV Proposed Method

This section describes the fundamental background of the two-phase dual COPOD
anomaly detection technique.

###  IV-A Two-phase Dual COPOD Method

Anomaly detection algorithms assume that all training data are true
representations of the normal operation of a system because the anomalies are
unknown ahead of time. The assumption may not hold for complex ICS with
several sensors and actuators. Sensor and actuator readings may be noisy. If
noisy data is used in training anomaly detection models, anomaly predictions
may not be accurate. Therefore, noise and apparent outliers need to be removed
in the data before training. To this end, a two-phase dual (TPD) COPOD method
is proposed. The proposed TPD COPOD is an anomaly detection method that
consists of two sequential modeling stages and a dual (parallel) modeling
stage. Figure 2 shows the high-level overview of the TPD COPOD. The first
phase is responsible for data cleaning and noise removal using ECOD. The
second stage of the method is trained using the cleaned data from phase 1. The
details of the individual stages of the method are discussed in the subsequent
sections.

Figure 2: Two-phase Dual COPOD Method.

###  IV-B Phase 1: ECOD Model

Phase 1 of the method uses the ECOD algorithm [14] to remove anomalies and
obvious noise in the training dataset so that the assumption about the
training dataset being an authentic representation of the normal ICS
operations holds. Anomalies are rare data points that occur in low-density
areas of a probability distribution [27, 41]. If the distribution is unimodal,
these rare events are found in the distribution’s tails. Determining the
likelihood of finding a data point at least as ”extreme” as $X_{i}$ in terms
of tail probabilities forms the basis of the approach ECOD algorithm of phase
1. In this work, it is assumed that for a given ICS dataset, $n$ data points
$X_{1},X_{2},...,X_{n}\in R^{d}$ are sampled independently and identically
distributed. Figure 3 shows the architecture of phase 1 in the TPD COPOD
method.

Figure 3: Phase 1 Architecture of the Proposed Method.

The joint cumulative distribution function over all $d$ dimensions is
represented by $F:R^{d}\to[0,1]$. Then, for a vector $z\in R^{d}$, let
$z^{(j)}$ represent the $j-th$ entry. Hence by definition of a joint CDF, for
any $x\in R^{d}$,

$F(x)=P(X^{(1)}\leq x^{(1)},X^{(2)}\leq x^{(2)},...,X^{(d)}\leq x^{(d)})$ (1)

Regarding the left tail, (1) determines how extreme $X_{i}$ is. The smaller
$F(X_{i})$ is, the less likely a point $X$ sampled from the same distribution
as $X_{i}$ will satisfy the inequality $X\leq X_{i}$. Similarly, $1-F(X_{i})$
measures the extremeness of $X_{i}$ by focusing on the right tails of each
dimension rather than the left tails. As a result, if either $F_{X}(X_{i})$ or
$1-F(X_{i})$ is extremely small, then this suggests that $X_{i}$ corresponds
to a rare event and is, therefore, likely to be an anomaly. To simplify,
assume that a dataset’s various dimensions are independent so that the
combined CDF can be factored

$F(x)=\prod_{j=1}^{d}F^{(J)}(x^{(j)})\quad\text{for}\quad x\in R^{d},$

where $F^{(j)}:R\to[0,1]$ represents the univariate CDF of the j-th dimension.
Now it is sufficient to notice that univariate CDFs can be accurately
estimated by utilizing the ECDF. The left tail ECDF is

$\hat{F}_{left}^{(j)}(z)=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{X_{i}^{(j)}\leq
z\\}\quad\text{for}\quad z\in R$ (2)

where $\mathbbm{1}\\{\cdot\\}$ is the indicator function that is 1 when its
argument is true and is 0 otherwise. Similarly, the right tail ECDF is

$\hat{F}_{right}^{(j)}(z)=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{X_{i}^{(j)}\geq
z\\}\quad\text{for}\quad z\in R.$ (3)

For every point $X_{i}^{(j)}$, the tail probabilities
$\hat{F}_{left}^{(j)}(X_{i}^{(j)})$ and $\hat{F}_{right}^{(j)}(X_{i}^{(j)})$
are aggregated by multiplying them together to achieve the final anomaly score
$O_{i}\in[0,\infty)$. Higher $O_{i}$ means more likely to be an anomaly.

The skewness of a dimension’s distribution is used in aggregating the tail
probabilities to determine whether the left or right tail probability should
be used for a given dimension. The sample skewness coefficient,
$\tilde{\mu_{3}}$ of dimension $j$, for a given dataset can be derived as [42]

$\tilde{\mu_{3}}=\frac{\sum_{i=1}^{n}(X_{i}^{(j)}-\bar{X^{(j)}})^{3}}{(n-1)\times\sigma^{3}}$
(4)

where $\sigma$, the standard deviation is

$\sigma=\sqrt{\frac{\sum_{i=1}^{n}(X_{i}^{(j)}-\bar{X}^{(j)})^{2}}{(n-1)}}$

and $\bar{X}^{(j)}$ is the mean of dimension $j$’s distribution, and
$\bar{X}^{(j)}=\sum_{i=1}^{n}X_{i}^{(j)}/n$. When $\tilde{\mu_{3}}<0$, points
on the left tail can be considered outliers, whereas when $\tilde{\mu_{3}}>0$,
points on the right tail can be considered outliers.

The tail probabilities are transformed to log negative probabilities to
calculate the final anomaly score per data point for a given training dataset.
The log negative probability to the left and right tails are

$O_{left}(X_{i})=-\sum_{j=1}^{d}\log(\hat{F}_{left}^{(j)}(X_{i}^{(j)})),$ (5)

and

$O_{right}(X_{i})=-\sum_{j=1}^{d}\log(\hat{F}_{right}^{(j)}(X_{i}^{(j)})),$
(6)

respectively.The automated form of selecting the left or right tail of the
$j-th$ dimension based on whether $\tilde{\mu_{3}}<0$ or $\tilde{\mu_{3}}>0$
is

$\begin{split}O_{auto}(X_{i})&=\sum_{j=1}^{d}[\mathbbm{1}\\{\tilde{\mu_{3}}<0\\}\log(\hat{F}_{left}^{(j)}(X_{i}^{(j)}))\\\
&+\mathbbm{1}\\{\tilde{\mu_{3}}\geq
0\\}\log(\hat{F}_{right}^{(j)}(X_{i}^{(j)}))].\end{split}$ (7)

The final anomaly score is calculated in the space of negative log
probabilities, where a lower probability translates into a larger negative log
probability and, as a result, a higher likelihood of being an anomaly. The
largest magnitude of the three computed anomaly scores is chosen as the final
output anomaly score $O_{i}$ for the data point $X_{i}$ with

$O_{i}=\max\\{O_{left}(X_{i}),O_{right}(X_{i}),O_{auto}(X_{i})\\}.$ (8)

Algorithm 1 contains the ECOD algorithm’s pseudocode, which is utilized in
phase 1 of the proposed method. Also, $O_{i}^{(j)}$ is the dimensional anomaly
score for dimension j of $X_{i}$. Since the log function is monotonic,
$O_{i}^{(j)}$ represents the degree of outlierliness and the anomaly
contribution by $X_{i}^{(j)}$. This representation creates model
interpretability for the ECOD model in phase 1 of the method.

Input: Training dataset $\\{X_{i}^{(J)}|i=1,2,3...,n\\}$ where $X_{i}^{(J)}$
refers to the $j-th$ feature (dimension) of the $i-th$ data point

Output: Set of anomaly scores $O_{i}\in R^{n}$

for  _each dimension j in 1,…,d_ do

Estimate left and right tail ECDFs (using equation 2 and IV-B)

Compute the sample skewness coefficient {$\tilde{\mu_{3}}$} for the j-th
feature’s distribution using (4)

end for

for  _each sample i in 1,…,n_ do

Aggregate tail probabilities of $X_{i}$ to obtain anomaly score $O_{i}$:

$O_{left}(X_{i})=-\sum_{j=1}^{d}\log(\hat{F}_{left}^{(j)}(X_{i}^{(j)}))$

$O_{right}(X_{i})=-\sum_{j=1}^{d}\log(\hat{F}_{right}^{(j)}(X_{i}^{(j)}))$

$O_{auto}(X_{i})=\sum_{j=1}^{d}[\mathbbm{1}\\{\tilde{\mu_{3}}<0\\}\log(\hat{F}_{left}^{(j)}(X_{i}^{(j)}))+\mathbbm{1}\\{\tilde{\mu_{3}}\geq
0\\}\log(\hat{F}_{right}^{(j)}(X_{i}^{(j)}))].$

Set the final anomaly score for point $X_{i}$ as:

$O_{i}=\max\\{O_{left}(X_{i}),O_{right}(X_{i}),O_{auto}(X_{i})\\}$.

end for

return _Anomaly scores, $O=(O_{1},...,O_{n})$_

Algorithm 1 Phase 1: Noise Reduction using ECOD.

###  IV-C Phase 2: Dual COPOD Model

Phase 2 of the TPD COPOD consists of two parallel COPOD models. This stage
receives the input data from phase 1, which is separated into discrete and
continuous data points. One COPOD model receives the discrete input, whereas
the other receives the continuous input. Phase 2 of the TPD COPOD model is
responsible for the actual anomaly detection after training. The details of
the phase 2 architecture are discussed in this section. Figure 4 shows phase 2
architecture of the proposed method.

Figure 4: Phase 2 Architecture of the Proposed Method.

Copulas are mathematical operations that allow us to decouple marginal
distributions from a given multivariate distribution’s dependence structure.
Formally, a d-variate copula $C:[0,1]^{d}\in[0,1]$, is the CDF of a random
vector $(U^{(1)},U^{(2)},...,U^{(d)})$ with uniform marginals given by

$C_{U}(u)=P(U^{(1)}\leq u^{(1)},...,U^{(d)}\leq u^{(d)})$ (9)

where $P(U^{(j)}\leq u^{(j)})=u^{(j)}$ for $j\in 1,...,d$ and
$u^{(j)}\in[0,1]$. By using inverse sampling, uniform distributions can be
transformed into any desired distributions, such as

$X_{j}=F_{j}^{-1}(U_{j})\sim F_{j}$ (10)

It has been shown in [43] that for any random variables $X_{1},\cdots,X_{d}$
with joint distribution function $F(x_{1},\cdots,x_{d})$ and marginal
distributions $F_{1},\cdots,F_{d}$, there exists a copula such that

$F(x)=C(F_{1}(x_{1}),\cdots,F_{d}(x_{d}))$ (11)

With copula, a joint distribution of $(X_{1},\cdots,X_{d})$ can be described
using only their marginals. This provides flexibility when modeling high-
dimensional ICS datasets because the dataset’s dimensions can be modeled
separately, and a linkage of marginal distributions together to form the joint
distribution is gauranteed. Furthermore [43] shows that if $F$ has univariate
marginal distributions $F_{1},\cdots,F_{d}$, then there exists a copula
$C(\cdot)$ such that (11) holds. Additionally, if the marginal distributions
are continuous, then $C(\cdot)$ can be uniquely determined. The expression for
the copula equation in terms of the joint CDF and inverse CDFs is derived by
replacing (9) with the inverse (10) to yield

$C(u)=F_{x}(F_{x_{d}}^{-1}(u_{1}),\cdots,F_{x_{d}}^{-1}(u_{1})).$ (12)

Together, (11) and (12) are referred to as Sklar’s Theorem, which ensures the
existence of a copula for any multivariate CDF with continuous marginals and
offers a closed form equation for constructing the copula.

The dual COPOD algorithm of phase 2 uses a nonparametric approach based on
fitting ECDFs called empirical copula [15]. Let $X_{i}^{(j)}$ represent the
i-th observation of the j-th feature/dimension for the given d-dimensional ICS
training dataset $X$ with $n$ observations. The empirical CDF is

$\hat{F}(x)=P((-\infty,x])=\frac{1}{n}\sum_{i=1}^{n}[\mathbbm{1}\\{X_{i}\leq
x_{i}\\}.$ (13)

It is possible to determine $\hat{U_{i}}$ by using the empirical copula
observations as the inverse of (10), such that

$(\hat{U}_{i}^{(1)},\cdots,\hat{U}_{i}^{d})=(\hat{F}^{(1)}(X_{i}^{(1)}),\cdots,\hat{F}^{(d)}(X_{i}^{(d)})).$
(14)

Finally, by substituting the empirical copula observations of (14) into the
first equality of (12), the empirical copulas are

$\hat{C}(u^{(1)},\cdots,u^{(d)})=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{\tilde{U}_{i}^{(1)}\leq
u^{(1)},\cdots,\tilde{U}_{i}^{(d)}\leq u^{(d)}\\}.$ (15)

An empirical copula, $\hat{C}(u)$ with multivariate CDFs supported on $n$
points in the grid $\\{1/n,2/n,\cdots,1\\}^{d}$, has discrete uniform
marginals on $\\{1/n,2/n,\cdots,1\\}$, and asymptotically converges to $C(u)$
as a result of the central limit theorem [44].

The COPOD algorithm adopts a three-stage process in detecting ICS anomalies.
First, the COPOD algorithm computes the ECDFs based on the input dataset
received from phase 1. Secondly, the COPOD algorithm uses the ECDFs to produce
the empirical copula function. Finally, the COPOD algorithm uses the empirical
copula in (15) to approximate the tail probability described in (11).

The COPOD algorithm’s objective is to determine the probability of detecting a
point at least as extreme as each observation in the input dataset, $x_{i}$.
That is, assume that $x_{i}$ is distributed according to d-distribution
function $F_{X}(X_{i})$, the COPOD algorithm needs to calculate the tail
probabilities, $F_{X}(x_{i})=P(X\leq x_{i})$ and $1-F(x_{i})=P(X\geq x_{i})$.
If $x_{i}$ is an anomaly, the probability of observing a point at least as
extreme as $x_{i}$ should be small. Therefore, if either $F_{X}(x_{i})$ or
$1-F(x_{i})$ is extremely small, point $x_{i}$ occurs infrequently and is
likely to be an anomaly. In the COPOD algorithm, $F_{X}(x_{i})$ is known as
the left tail probability of $x_{i}$, and $1-F(x_{i})$ is known as the right
tail probability of $x_{i}$. Therefore, an anomaly is considered as an
observation that has a small tail probability if either of the two quantities
($F_{X}(x_{i})$ and $1-F(x_{i})$) are small.

High-dimensional data spaces have challenges not found in low-dimensional
spaces. These challenges are known as the curse of dimensionality, and are
common in anomaly detection domain [15]. In order to prevent diminishing tail
probabilities and to exploit the monotonicity property of the $\log()$
function, the COPOD algorithm uses the sum of the negative log probabilities
similar to the ECOD algorithm of phase 1 in (5) and (6).

###  IV-D The Dual COPOD Approach to Anomaly Detection

The cleaned dataset from phase 1 goes to phase 2 for training the dual COPOD
models. In phase 2, the input data is separated into discrete and continuous
data points. The discrete data points go to the upper COPOD model (COPOD 1 of
Figure 4), whereas the continuous data is input to the lower COPOD model
(COPOD 2 of Figure 4). The reason for the dual architecture model is to allow
phase 2 of the architecture to exploit ICS data by examining the use of two
latent representations to extract useful information to minimize overfitting
and improve the model’s anomaly detection capability.

Each of the COPOD models in the dual architecture in Figure 4 requires a
d-dimensional input dataset from phase 1;
$X=(X_{i}^{(1)},X_{i}^{(2)},\cdots,X_{i}^{(d)})$, where $i=1,\cdots,n$, and
produces an anomaly score vector $O(X)=[X_{1},X_{2},\cdots,X_{n}]$. The
anomaly scores are between $(0,\infty)$, and are to be used comparatively. The
anomaly score does not indicate the probability of $X_{i}$ being an anomaly
but rather the relative measure of how likely $X_{i}$ is when compared to
other points in the dataset. Larger $O(X_{i})$ signifies $X_{i}$ is more
likely anomalous.

Each COPOD model fits d-dimensional left tail CDFs, using (13) and
d-dimensional right tail CDFs by replacing $X$ in (13) with $-X$. Also,
d-dimensional skewness vector $\tilde{\mu_{3}}$ is computed using (4). Next,
the empirical copula observations for each $X_{i}$ are computed using (14) to
obtain the left tail copulas $\hat{U}_{left}^{j}$ and right tail copulas
$\hat{U}_{right}^{j}$. Then, the skewness corrected empirical copula
observations are calculated as

$\hat{W}_{i}^{(j)}=\begin{cases}\hat{U}_{left}^{j}&\text{if}\quad\tilde{\mu_{3}}<0\\\
\hat{U}_{right}^{j}&\text{otherwise.}\\\ \end{cases}$

Finally, the probability of observing a point at least as extreme as each
$x_{i}$ along each dimension is computed. Similar to (8), the maximum of the
negative log of the probabilities generated by the left tail empirical copula,
right tail empirical copula, and skewness corrected empirical copula is
selected as the final anomaly score. That is, the smaller the tail probability
is, the bigger its negative log, and so a data point is considered an outlier
if it has a small left tail probability, a small right tail probability, or a
small skewness corrected tail probability.

The final anomaly scores of the TPD COPOD are generated as a combination of
the anomaly scores of the two COPOD models using a window function. The window
function defines the time intervals based on which anomaly scores can be
divided using a moving average [45]. Let $O_{c_{1}}(X_{i})$ and
$O_{c_{2}}(X_{i})$ be the anomaly scores of COPOD model 1 and COPOD model 2,
respectively. For a given data point $X_{i}$, its decision label is negative
(anomalous) if and only if the predictions of $O_{c_{1}}(X_{i})$ and
$O_{c_{2}}(X_{i})$ are negative; otherwise, the decision label of $X_{i}$ is
positive. The final decision score/label $O_{f}(X_{i})$ is made by observing
the predictions of the dual COPOD models for a time frame of $t_{w}$ time
instants. In this work, a $t_{w}$ of 30s is used for inference. Algorithm 2
summarizes the decision function of the TPD COPOD method. Combining the TPD
COPOD architecture with the algorithm definition of phase 2 results in the
complete TPD COPOD architecture as shown in Figure 4.

Input 1: Discrete input data $\\{X_{i}^{(J)}|i=1,2,3...,n\\}$ where
$X_{i}^{(J)}$ refers to the $j-th$ feature (dimension) of the $i-th$ discrete
data point

Input 2: Continuous input data $\\{X_{i}^{(k)}|i=1,2,3...,n\\}$ where
$X_{i}^{(k)}$ refers to the $k-th$ feature (dimension) of the $i-th$
continuous data point

Output: Set of decision labels $O_{f}$

for  _each dimension j in 1,…,d_ do

Compute anomaly scores $O_{c_{1}}(X_{i}^{(j)})$ for COPOD 1

end for

for  _each dimension k in 1,…,d_ do

Compute anomaly scores $O_{c_{2}}(X_{i}^{(k)})$ for COPOD 2

end for

for  _each data record i in 1,…,n_ do

if  _$O_{c_{1}}(X_{i})\parallel O_{c_{2}}(X_{i})==-1(anomalous)$_ then

Output $O_{i}=1$; ;

/* Anomalous data point */

else

Output $O_{i}=0$ ;

/* Normal data point */

end if

end for

for  _each data record i in 1,…,n_ do

Using a time frame, $t_{w}$

if  _$\frac{1}{t_{w}}\sum_{i=1}^{t_{w}}O_{i}\geq 80\%$_ then

Final decision label $O_{f_{i}}=-1$ ;

/* Anomaly detected */

else

Final decision label $O_{f_{i}}=1$ ;

/* Normal operation */

end if

end for

return _Final decision labels $\\{O_{f_{i}}\\}$_

Algorithm 2 TPD COPOD Method’s Decision Function.

##  V Results and Discussions

The TPD COPOD method’s performance evaluation is based on performance metrics,
results from predictions on the test data, and comparison with prior work
trained on similar datasets. The proposed method was developed using the
Python programming language, and PyOD which is the most comprehensive and
scalable Python library for detecting outlying objects in multivariate data
[46]. Evaluation results and performance metrics calculations are performed by
using the Scikit-learn library [47]. The three different datasets, namely;
SWaT [40], WADI [38], and TLIGHT [4] are used for evaluating the model
performance of the proposed method.

###  V-A Results Evaluation on the SWaT Dataset

This subsection describes the proposed method’s performance on the SWaT
dataset as compared with previous work evaluated on the same dataset. During
training, the SWaT training dataset is passed through the first stage of the
TPD COPOD method for any noise and unwanted signals to be removed. Figure 5
shows the first sample of noise detected by the ECOD model of phase 1 TPD
COPOD. Figure 5 shows the feature-level outlier scores explaining the reason
for detecting the first sample as an outlier. The x-axis indicates the
features (sensors and actuators) of the input dataset represented as numerical
integers. The blue dashed line represents the 90th percentile band, and the
orange dashed line represents the 99th percentile band for the given features.
Sample 1 of Figure 5 is flagged as an anomaly because several dimensions
(features), such as dimensions 2, 6-8, 19, 24, 28, 29, 31, 34, 35 39-43, and
45-48, have outlier scores that exceed the 99th percentile band. The detected
anomalous training sample 1 shown in Figure 5 is selected at random, and the
explanations provided about the figure is similar to the rest of the plots,
and are not shown in this paper. Figure 5 shows that about half of the
dimensions have outlier scores that exceed the 99th percentile band, and this
explains why the sample is detected as an anomaly. Overall, the computer used
for analysis required 16.70s for phase 1 to clean the training data consisting
of 496,800 samples. A total of 49,680 training samples were detected as
outliers in phase 1 of TPD COPOD.

Figure 5: Detected Noise in the SWaT Training Dataset by Phase 1 shown in
Figure 3 of the Proposed Method: Sample 1.

Figure 6 shows the first normal training sample in the SWaT training dataset
detected by the ECOD model of phase 1 TPD COPOD model. The training sample
shown in Figure 6 is detected as a normal sample because none of the
dimensional outlier scores exceeds the 99th percentile band. Similarly, in all
cases where the ECOD model of phase 1 TPD COPOD model detected a training
sample as normal, none of the dimensional anomaly scores of those samples
exceeded the 99th percentile band.

Figure 6: Detected Normal Data Sample in the SWaT Training Dataset by Phase 1
shown in Figure 3 of the Proposed Method: Sample 1.

The clean SWaT dataset consisting of 447,120 samples were separated into
discrete and continuous samples for phase 2 of the TPD COPOD. Training
required 2.35s for phase 2 models, and predictions required 6.02s over test
data consisting of 450,819 samples. After training, only phase 2 of the
proposed method is used for performing inference. Figure 7 shows the first
discrete sample of a detected anomaly in the test data by phase 2 COPOD 1.
Figure 7 shows that because of the discrete nature of the input dataset to
phase 2 COPOD 1, the 90th percentile band for most of the dimensions are 0.
Dimension 9 produced an outlier score that is equal to the 99th percentile
band, and as a result, dimension 9 is a major contributing factor to
classifying the first discrete test sample as an anomaly.

Figure 7: Detected Anomaly in the SWaT Test Dataset by Phase 2 shown in Figure
4 of the Proposed Method: Discrete Sample 1.

The results presented in Figure 8 use the visualization approach proposed in
[4, 48] to better understand the TPD COPOD model performance. The histogram-
based visualization approach normalizes the histogram frequency (y-axis) to a
range between $0\%$ and $100\%$, whereas the x-axis represents the model’s
normalized decision scores indicating the prediction confidence. Figure 8
shows the results of the normalized TP, TN, FP, and FN values on the SWaT test
set by phase 2 COPOD 1 shown in Figure 4. Figure 8 shows that phase 2 COPOD 1
correctly detected about 60% of the SWaT anomalies with over 60% confidence
level. Phase 2 COOPD 1 misclassified less than 10% of the normal instance as
anomalies, whereas about 35% of the anomalies were misclassified as normal
data points. Also, phase 2 COPOD 1 correctly predicted about 45% of the normal
data with prediction confidence of less than 10%. Figure 9 shows the results
of the normalized TP, TN, FP, and FN values on the SWaT test set by phase 2
COPOD 2 shown in Figure 4. Figure 9 shows that about 80% of the anomalies were
correctly detected by phase 2 COPOD 2 with prediction confidence between 95%
and 50%. The TP and FP prediction by phase 2 COPOD are normally distributed
with average prediction confidence of about 50%.

Figure 8: Results of the Normalized TP, TN, FP, and FN Values on the SWaT Test
Set by Phase 2 COPOD 1 shown in Figure 4.

Figure 9: Results of the Normalized TP, TN, FP, and FN Values on the SWaT Test
Set by Phase 2 COPOD 2 shown in Figure 4.

The SWaT dataset has been utilized for anomaly detection by K-Nearest
Neighbors (KNN) [49], Feature Bagging (FB) [49], Support Vector Machines (SVM)
[49], Auto-encoders (AE) [49] and Dual Isolation Forest (DIF) [9]. Table I
shows a comparison of evaluation results between previously published methods
and the proposed method. Table I shows that TPD COPOD has an F1-score of
$93\%$, precision of $93\%$, and recall of $93\%$. TPD COPOD has the highest
F1-score and recall values as compared with previous work. In terms of
precision, TPD COPOD performed at par with NN-one class, DIF, and SVM. The
high recall and F1-score values of TPD COPOD reflect its robustness and
ability to confidently detect SWaT dataset anomalies compared to previous
work. When compared to other approaches with similar computational
complexities, such as NN, the proposed anomaly detection method significantly
outperforms NN in terms of F1-score and recall. The outstanding performance of
TPD COPOD across all performance metrics could be attributed to the noise
removal from the training data in phase 1 during model training. Phase 2 of
the method in Figure 4 trains on clean training data, arguably a true
representation of the normal operations of the SWaT testbed.

TABLE I: Comparison between Detection Methods Evaluated on the SWaT Dataset Method |  F1-Score |  Precision |  Recall |  Complexity
---|---|---|---|---
NN [50] |  0.812 |  0.976 |  0.696 |  Low
RNN [51] |  0.802 |  0.982 |  0.678 |  High
SVM [51] |  0.796 |  0.925 |  0.699 |  High
TABOR(TB) [52] |  0.823 |  0.862 |  0.788 |  Average
ID-CNN [13] |  0.860 |  0.867 |  0.854 |  High
AE [49] |  0.520 |  0.516 |  0.516 |  Average
FB [49] |  0.360 |  0.358 |  0.358 |  Average
KNN [49] |  0.350 |  0.348 |  0.348 |  Average
GAN [49] |  0.510 |  0.406 |  0.677 |  High
DIF [9] |  0.882 |  0.935 |  0.835 |  Average
NN-one class [53] |  0.870 |  0.940 |  0.820 |  Average
TPD COPOD |  0.930 |  0.930 |  0.930 |  Very low

###  V-B Results Evaluation on the WADI Dataset

This subsection describes the proposed method’s performance on the WADI
dataset compared to previous work evaluated on the same dataset. During
training, the WADI training dataset is passed through the first stage of the
TPD COPOD method for any noise and unwanted signals to be removed. The total
number of dimensions (features) consisting of both discrete and continuous
values of the training dataset is 119. The clean WADI data consisting of
706,114 samples were separated into discrete and continuous samples for phase
2 of the TPD COPOD. The computer used for analysis required 13.80s to train
the phase 2 models and 19.64s to make predictions on the test data consisting
of 172800 samples.

The WADI dataset has been utilized for anomaly detection by SVM [49] KNN [49],
AE [49], FB [49], EGAN [49], GAN [49], Deep Autoencoding Gaussian Model
(DAGMM) [54], Long-short Term Memory Variational Autoencoder (LSTM-VAE) [55]
and DIF [9], and Graph Deviation Network (GDN) [56]. Table II shows a
comparison of evaluation results between previous work and the proposed
method. Table II shows that TPD COPOD has superior performance to all previous
work in terms of F1-score and recall. TPD COPOD achieved an F1-score of
$92\%$, precision of $82\%$, and recall of $86\%$. The high recall and
F1-score values of TPD COPOD reflect its robustness and ability to confidently
detect WADI dataset anomalies as compared to previous work. Although LSTM-VAE
and GDN achieved higher precision values as compared to the proposed method in
this study, LSTM-VAE and GDN models have significantly low recall values of
$27\%$ and $40.2\%$, respectively. The low recall values of LSTM-VAE and GDN
resulted in poor F1-scores of $25\%$ and $57\%$, respectively. The outstanding
performance of TPD COPOD across all performance metrics could be attributed to
the noise removal from the training data in phase 1 during model training.
Phase 2 of the method in Figure 4 trains on clean training data, arguably a
true representation of the normal operation of the WADI testbed.

TABLE II: Comparison between Detection Methods Evaluated on the WADI Dataset Method |  F1-Score |  Precision |  Recall |  Complexity
---|---|---|---|---
PCA [50] |  0.250 |  0.504 |  0.166 |  Average
SVM [51] |  0.510 |  0.512 |  0.512 |  High
AE [49] |  0.520 |  0.520 |  0.520 |  Average
FB [49] |  0.340 |  0.336 |  0.336 |  Average
KNN [49] |  0.300 |  0.299 |  0.299 |  Average
EGAN [49] |  0.340 |  0.345 |  0.345 |  High
GAN [49] |  0.620 |  0.538 |  0.749 |  High
DIF [9] |  0.656 |  0.765 |  0.574 |  Average
DAGMM [54] |  0.360 |  0.544 |  0.270 |  Average
LSTM-VAE [55] |  0.250 |  0.878 |  0.145 |  High
GDN [56] |  0.570 |  0.975 |  0.402 |  High
TPD COPOD |  0.920 |  0.820 |  0.860 |  Very low

###  V-C Results Evaluation on the TLIGHT Dataset

This subsection describes the proposed method’s performance on the TLIGHT
dataset as compared with previous work evaluated on the same dataset. During
training, the TLIGHT training dataset is passed through the first stage of the
TPD COPOD method for any noise and unwanted signals to be removed. The total
number of dimensions (features) consisting of both discrete and continuous
values of the training dataset is 36. There are five different test sets of
TLIGHT dataset. The computer used for analysis required 0.67s for phase 1 to
clean the training data consisting of 41,580 samples. A total of 4,144
training samples were detected as outliers in phase 1.

Table III shows the comparison between TPD COPOD and previous methods
evaluated on TLIGHT dataset. Test set 1 of the TLIGHT dataset consisted of
5,000 samples, the computer required 0.328s to make predictions. Table III
shows that TPD COPOD achieves superior precision, recall, and F1-score of
$96\%$, $96\%$, and $95\%$ respectively. OCSVM [4], OCNN [4] and IF [4] have
similar performance in terms of precision, recall, and F1-score. The test set
2 of the TLIGHT dataset consisted of 7,000 samples. The computer used for
analysis required 0.0.437s to make predictions on the test set 2. TPD COPOD
achieves superior performance by correctly predicting all the data points in
test set 2 by achieving $100\%$ precision, recall, and F1-score. OCSVM, OCNN
and IF have similar performance in terms of precision, recall, and F1-score.
IF achieves high precision, recall, and F1-score values of $98\%$, $97\%$, and
$97\%$ respectively, on test set 2. However, OCNN and OCSVM achieved similar
performance on test set 2. The test set 3 of the TLIGHT dataset consisted of
13,130 samples. The computer used for analysis required 0.324s to make
predictions on the test set 3. TPD COPOD achieves superior precsion, recall,
and F1-score of $91\%$, $90\%$, and $88\%$ respectively, on test set 3. OCSVM,
OCNN and IF achieve similar performance in terms of precision, recall, and
F1-score on test set 3.

The test set 4 of the TLIGHT dataset consisted of 15,000 samples. The computer
used for analysis required 0.477s to make predictions on the test set 4. TPD
COPOD achieved relatively low performance on test set 4 as compared to its
performance on test set 1 and 2 by having precision, recall, and F1-score
values of $88\%$, $85\%$, and $83\%$ respectively. OCNN, IF, and TPD COPOD
achieve similar performance whereas OCSVM has the worst performance by having
precision, recall, and F1-score values of $81\%$, $81\%$, and $71\%$
respectively. TPD COPOD has low performance on test set 4 because of the large
proportion of anomalies in test set 4 consisting of timing bits anomalies
which are hard to detect [4]. The test set 5 of the TLIGHT dataset consisted
of 18,269 samples. The computer used for analysis required 0.517 to make
predictions on the test set 5. TPD COPOD achieved its lowest performance on
test set 5 by having precision, recall, and F1-score values of $83\%$, $75\%$,
and $73\%$ respectively. OCNN and OCSVM achieve similar performance whereas IF
has the best performance by having precision, recall, and F1-score values of
$82\%$, $78\%$, and $77\%$ respectively. Again, TPD COPOD has low performance
on test set 5 because of the large proportion of anomalies in test set 5
consisting of timing bits anomalies. Therefore, the TPD COPOD appears to be
ineffective at detecting TLIGHT system errors consisting of system timing
bits.

TABLE III: Comparison between Detection Methods Evaluated on TLIGHT Dataset Dataset |  Method |  Precision |  Recall |  F1-Score |  Complexity
---|---|---|---|---|---
Test |  OCSVM |  0.89 |  0.90 |  0.98 |  High
Set |  OCNN |  0.90 |  0.91 |  0.90 |  Average
1 |  IF |  0.90 |  0.91 |  0.91 |  Average
|  TPD COPOD |  0.96 |  0.96 |  0.95 |  Very low
Test |  OCSVM |  0.81 |  0.87 |  0.89 |  High
Set |  OCNN |  0.81 |  0.88 |  0.84 |  Average
2 |  IF |  0.98 |  0.97 |  0.97 |  Average
|  TPD COPOD |  1.00 |  1.00 |  1.00 |  Very low
Test |  OCSVM |  0.85 |  0.86 |  0.84 |  High
Set |  OCNN |  0.87 |  0.88 |  0.86 |  Average
3 |  IF |  0.87 |  0.87 |  0.87 |  Average
|  TPD COPOD |  0.91 |  0.90 |  0.88 |  Very low
Test |  OCSVM |  0.81 |  0.81 |  0.71 |  High
Set |  OCNN |  0.87 |  0.88 |  0.86 |  Average
4 |  IF |  0.86 |  0.86 |  0.85 |  Average
|  TPD COPOD |  0.88 |  0.85 |  0.83 |  Very low
Test |  OCSVM |  0.78 |  0.70 |  0.68 |  High
Set |  OCNN |  0.79 |  0.70 |  0.68 |  Average
5 |  IF |  0.82 |  0.78 |  0.77 |  Average
|  TPD COPOD |  0.83 |  0.75 |  0.73 |  Very low

##  VI Conclusion

This work proposes the first known two-phase dual (TPD) COPOD anomaly
detection method, which is an unsupervised anomaly detection technique
consisting of two sequential stages and a dual (parallel) modeling stage.
Phase 1 of the method makes use of ECOD to remove any obvious noise or outlier
data records in a given training dataset. Phase 2 of the method consists of a
dual COPOD architecture that utilizes the output data of phase 1 to develop
two COPOD models. The algorithms implemented in the proposed method are
parameter-free and based on empirical distribution functions. The
deterministic nature of all stages of the method results in the mitigation of
the challenges associated with hyperparameter selection in unsupervised
anomaly detection. Furthermore, the proposed method is highly interpretable
and quantifies each feature’s contribution toward an ICS anomaly. The proposed
anomaly detection method is computationally and memory efficient, scalable,
and suitable for low- and high-dimensional ICS datasets. The proposed method
is trained, evaluated, and compared with previous work using three open-source
ICS datasets. The proposed method outperformed previous work in terms of
F1-score and recall on the SWaT, WADI and TLIGHT datasets. The robust
performance of the TPD COPOD method coupled with its speed of anomaly
detection makes the TPD COPOD capable of real-time ICS anomaly detection.
Future work should focus on finding a means by which the ECOD algorithm of
phase 1 may be extended to multimodal training dataset.

## References

  * [1] G. R. MR, C. M. Ahmed, and A. Mathur, “Machine learning for intrusion detection in industrial control systems: challenges and lessons from experimental evaluation,” _Cybersecurity_ , vol. 4, no. 1, pp. 1–12, 2021\.
  * [2] S. Karnouskos, “Stuxnet worm impact on industrial cyber-physical system security,” in _IECON 2011-37th Annual Conference of the IEEE Industrial Electronics Society_. IEEE, 2011, pp. 4490–4494.
  * [3] D. Alert, “Cyber-attack against ukrainian critical infrastructure,” _Cybersecurity Infrastruct. Secur. Agency, Washington, DC, USA, Tech. Rep. ICS Alert (IR-ALERT-H-16-056-01)_ , 2016.
  * [4] E. Aboah Boateng and J. Bruce, “Unsupervised machine learning techniques for detecting plc process control anomalies,” _Journal of Cybersecurity and Privacy_ , vol. 2, pp. 220–244, 2022.
  * [5] E. Aboah Boateng and J. W. Bruce, “Anomaly detection for industrial control systems based on neural networks with one-class objective function,” _Proceedings of Student Research and Creative Inquiry Day_ , vol. 5, 2021\.
  * [6] J. Liu, L. Yin, Y. Hu, S. Lv, and L. Sun, “A novel intrusion detection algorithm for industrial control systems based on cnn and process state transition,” in _2018 IEEE 37th International Performance Computing and Communications Conference (IPCCC)_. IEEE, 2018, pp. 1–8.
  * [7] V. L. Cao, M. Nicolau, and J. McDermott, “One-class classification for anomaly detection with kernel density estimation and genetic programming,” in _European Conference on Genetic Programming_. Springer, 2016, pp. 3–18.
  * [8] J. Vavra, M. Hromada, L. Lukas, and J. Dworzecki, “Adaptive anomaly detection system based on machine learning algorithms in an industrial control environment,” _International Journal of Critical Infrastructure Protection_ , vol. 34, p. 100446, 2021.
  * [9] M. Elnour, N. Meskin, K. Khan, and R. Jain, “A dual-isolation-forests-based attack detection framework for industrial control systems,” _IEEE Access_ , vol. 8, pp. 36 639–36 651, 2020.
  * [10] Z. Zhao, K. G. Mehrotra, and C. K. Mohan, “Ensemble algorithms for unsupervised anomaly detection,” in _International Conference on Industrial, Engineering and Other Applications of Applied Intelligent Systems_. Springer, 2015, pp. 514–525.
  * [11] M. Ganaie, M. Hu _et al._ , “Ensemble deep learning: A review,” _arXiv preprint arXiv:2104.02395_ , 2021.
  * [12] E. Aboah Boateng and J. Bruce, “Unsupervised ensemble methods for anomaly detection in plc-based process control,” _arXiv e-prints_ , pp. arXiv–2302, 2023.
  * [13] M. Kravchik and A. Shabtai, “Detecting cyber attacks in industrial control systems using convolutional neural networks,” in _Proceedings of the 2018 workshop on cyber-physical systems security and privacy_ , 2018, pp. 72–83.
  * [14] Z. Li, Y. Zhao, X. Hu, N. Botta, C. Ionescu, and G. Chen, “Ecod: Unsupervised outlier detection using empirical cumulative distribution functions,” _IEEE Transactions on Knowledge and Data Engineering_ , 2022.
  * [15] Z. Li, Y. Zhao, N. Botta, C. Ionescu, and X. Hu, “Copod: copula-based outlier detection,” in _2020 IEEE International Conference on Data Mining (ICDM)_. IEEE, 2020, pp. 1118–1123.
  * [16] J. Giraldo, D. Urbina, A. Cardenas, J. Valente, M. Faisal, J. Ruths, N. O. Tippenhauer, H. Sandberg, and R. Candell, “A survey of physics-based attack detection in cyber-physical systems,” _ACM Computing Surveys (CSUR)_ , vol. 51, no. 4, pp. 1–36, 2018.
  * [17] S. Mishra, Y. Shoukry, N. Karamchandani, S. N. Diggavi, and P. Tabuada, “Secure state estimation against sensor attacks in the presence of noise,” _IEEE Transactions on Control of Network Systems_ , vol. 4, no. 1, pp. 49–59, 2016.
  * [18] C. Murguia and J. Ruths, “Characterization of a cusum model-based sensor attack detector,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016, pp. 1303–1309.
  * [19] D. Mashima and A. A. Cárdenas, “Evaluating electricity theft detectors in smart grid networks,” in _International Workshop on Recent Advances in Intrusion Detection_. Springer, 2012, pp. 210–229.
  * [20] M. S. Rahman, M. A. Mahmud, A. M. T. Oo, and H. R. Pota, “Multi-agent approach for enhancing security of protection schemes in cyber-physical energy systems,” _IEEE transactions on industrial informatics_ , vol. 13, no. 2, pp. 436–447, 2016.
  * [21] T. Roth and B. McMillin, “Physical attestation in the smart grid for distributed state verification,” _IEEE Transactions on Dependable and Secure Computing_ , vol. 15, no. 2, pp. 275–288, 2016.
  * [22] V. K. Mishra, V. R. Palleti, and A. Mathur, “A modeling framework for critical infrastructure and its application in detecting cyber-attacks on a water distribution system,” _International Journal of Critical Infrastructure Protection_ , vol. 26, p. 100298, 2019.
  * [23] X. Yang, L. J. Latecki, and D. Pokrajac, “Outlier detection with globally optimal exemplar-based gmm,” in _Proceedings of the 2009 SIAM International Conference on Data Mining_. SIAM, 2009, pp. 145–154.
  * [24] M. H. Satman, “A new algorithm for detecting outliers in linear regression,” _International Journal of statistics and Probability_ , vol. 2, no. 3, p. 101, 2013.
  * [25] M. Goldstein and A. Dengel, “Histogram-based outlier score (hbos): A fast unsupervised anomaly detection algorithm,” _KI-2012: poster and demo track_ , vol. 9, 2012.
  * [26] M. Pavlidou and G. Zioutas, “Kernel density outlier detector,” in _Topics in Nonparametric Statistics_. Springer, 2014, pp. 241–250.
  * [27] A. Lazarevic and V. Kumar, “Feature bagging for outlier detection,” in _Proceedings of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining_ , 2005, pp. 157–166.
  * [28] F. T. Liu, K. M. Ting, and Z.-H. Zhou, “Isolation forest,” in _2008 eighth ieee international conference on data mining_. IEEE, 2008, pp. 413–422.
  * [29] Y. Zhao, Z. Nasrullah, M. K. Hryniewicki, and Z. Li, “Lscp: Locally selective combination in parallel outlier ensembles,” in _Proceedings of the 2019 SIAM International Conference on Data Mining_. SIAM, 2019, pp. 585–593.
  * [30] Y. Zhao, X. Hu, C. Cheng, C. Wang, C. Wan, W. Wang, J. Yang, H. Bai, Z. Li, C. Xiao _et al._ , “Suod: Accelerating large-scale unsupervised heterogeneous outlier detection,” _Proceedings of Machine Learning and Systems_ , vol. 3, pp. 463–478, 2021.
  * [31] L. Sun, S. Versteeg, S. Boztas, and A. Rao, “Detecting anomalous user behavior using an extended isolation forest algorithm: an enterprise case study,” _arXiv preprint arXiv:1609.06676_ , 2016.
  * [32] P. Parveen, N. Mcdaniel, Z. Weger, J. Evans, B. Thuraisingham, K. Hamlen, and L. Khan, “Evolving insider threat detection stream mining perspective,” _International Journal on Artificial Intelligence Tools_ , vol. 22, no. 05, p. 1360013, 2013.
  * [33] Y. Zhao, R. Rossi, and L. Akoglu, “Automatic unsupervised outlier model selection,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 4489–4502, 2021.
  * [34] M. Zhou, Y. Wang, A. K. Srivastava, Y. Wu, and P. Banerjee, “Ensemble-based algorithm for synchrophasor data anomaly detection,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 3, pp. 2979–2988, 2018.
  * [35] Y.-j. Xiao, W.-y. Xu, Z.-h. Jia, Z.-r. Ma, and D.-l. Qi, “Nipad: a non-invasive power-based anomaly detection scheme for programmable logic controllers,” _Frontiers of Information Technology & Electronic Engineering_, vol. 18, no. 4, pp. 519–534, 2017.
  * [36] A. N. Sokolov, A. N. Ragozin, I. A. Pyatnitsky, and S. K. Alabugin, “Applying of digital signal processing techniques to improve the performance of machine learning-based cyber attack detection in industrial control system,” in _Proceedings of the 12th International Conference on Security of Information and Networks_ , 2019, pp. 1–4.
  * [37] A. P. Mathur and N. O. Tippenhauer, “Swat: A water treatment testbed for research and training on ics security,” in _2016 international workshop on cyber-physical systems for smart water networks (CySWater)_. IEEE, 2016, pp. 31–36.
  * [38] S. Adepu and A. Mathur, “Generalized attacker and attack models for cyber physical systems,” in _2016 IEEE 40th annual computer software and applications conference (COMPSAC)_ , vol. 1. IEEE, 2016, pp. 283–292.
  * [39] S. Siemens, “S7-300 programmable controller quick start, primer, preface,” C79000-G7076-C500-01, Nuremberg, Germany, Tech. Rep., 1996.
  * [40] J. Goh, S. Adepu, M. Tan, and Z. S. Lee, “Anomaly detection in cyber physical systems using recurrent neural networks,” in _2017 IEEE 18th International Symposium on High Assurance Systems Engineering (HASE)_. IEEE, 2017, pp. 140–145.
  * [41] D. Pokrajac, A. Lazarevic, and L. J. Latecki, “Incremental local outlier detection for data streams,” in _2007 IEEE symposium on computational intelligence and data mining_. IEEE, 2007, pp. 504–515.
  * [42] R. A. Groeneveld and G. Meeden, “Measuring skewness and kurtosis,” _Journal of the Royal Statistical Society: Series D (The Statistician)_ , vol. 33, no. 4, pp. 391–399, 1984.
  * [43] M. Sklar, “Fonctions de repartition an dimensions et leurs marges,” _Publ. inst. statist. univ. Paris_ , vol. 8, pp. 229–231, 1959.
  * [44] R. B. Nelsen, _An introduction to copulas_. Springer Science & Business Media, 2007.
  * [45] R. J. Hyndman, “Moving averages.” 2011.
  * [46] Y. Zhao, Z. Nasrullah, and Z. Li, “Pyod: A python toolbox for scalable outlier detection,” _arXiv preprint arXiv:1901.01588_ , 2019.
  * [47] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * [48] E. Aboah Boateng and J. Bruce, “A new histogram-based visualization tool for analyzing anomaly detection algorithm performance,” _Proceedings of Student Research and Creative Inquiry Day_ , vol. 6, 2022.
  * [49] D. Li, D. Chen, J. Goh, and S.-k. Ng, “Anomaly detection with generative adversarial networks for multivariate time series,” _arXiv preprint arXiv:1809.04758_ , 2018.
  * [50] D. Shalyga, P. Filonov, and A. Lavrentyev, “Anomaly detection for water treatment system based on neural network with automatic architecture optimization,” _arXiv preprint arXiv:1807.07282_ , 2018.
  * [51] J. Inoue, Y. Yamagata, Y. Chen, C. M. Poskitt, and J. Sun, “Anomaly detection for a water treatment system using unsupervised machine learning,” in _2017 IEEE International Conference on Data Mining Workshops (ICDMW)_. IEEE, 2017, pp. 1058–1065.
  * [52] Q. Lin, S. Adepu, S. Verwer, and A. Mathur, “Tabor: A graphical model-based approach for anomaly detection in industrial control systems,” in _Proceedings of the 2018 on asia conference on computer and communications security_ , 2018, pp. 525–536.
  * [53] E. A. Boateng, J. Bruce, and D. A. Talbert, “Anomaly detection for a water treatment system based on one-class neural network,” _IEEE Access_ , vol. 10, pp. 115 179–115 191, 2022.
  * [54] B. Zong, Q. Song, M. R. Min, W. Cheng, C. Lumezanu, D. Cho, and H. Chen, “Deep autoencoding gaussian mixture model for unsupervised anomaly detection,” in _International conference on learning representations_ , 2018.
  * [55] D. Park, Y. Hoshi, and C. C. Kemp, “A multimodal anomaly detector for robot-assisted feeding using an lstm-based variational autoencoder,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1544–1551, 2018\.
  * [56] A. Deng and B. Hooi, “Graph neural network-based anomaly detection in multivariate time series,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 5, 2021, pp. 4027–4035.