# Active Attack Detection and Control in Constrained Cyber-Physical Systems
Under Prevented Actuation Attack

Mehdi Hosseinzadeh, , and Bruno Sinopoli This research has been supported by
National Science Foundation under award numbers ECCS-1932530.M. Hosseinzadeh
and B. Sinopoli are with the Department of Electrical and Systems Engineering,
Washington University in St. Louis, St. Louis, Missouri, USA (email:
<EMAIL_ADDRESS>bsinopoli@wustl.edu).

###### Abstract

This paper proposes an active attack detection scheme for constrained cyber-
physical systems. Despite passive approaches where the detection is based on
the analysis of the input-output data, active approaches interact with the
system by designing the control input so to improve detection. This paper
focuses on the prevented actuation attack, where the attacker prevents the
exchange of information between the controller and actuators. The proposed
scheme consists of two units: 1) detection, and 2) control. The detection unit
includes a set of parallel detectors, which are designed based on the
multiple-model adaptive estimation approach to detect the attack and to
identify the attacked actuator(s). For what regards the control unit, a
constrained optimization approach is developed to determine the control input
such that the control and detection aims are achieved. In the formulation of
the detection and control objective functions, a probabilistic approach is
used to reap the benefits of the a priori information availability. The
effectiveness of the proposed scheme is demonstrated through a simulation
study on an irrigation channel.

##  I Introduction

Cyber-Physical Systems (CPSs) often employ distributed networks of embedded
sensors and actuators that interact with the physical environment. The
availability of cheap communication technologies (e.g., internet) has
certainly improved scalability and functionality features in several
applications. However, they have made CPSs susceptible to cyber security
threats. This makes the cyber security to be of primary importance in safe
operation of CPSs.

By assuming that sensors-to-controller and controller-to-actuators
communication channels are the only ones in CPSs executed via internet and
malicious agents can alter data flows in these channels, in general two
classes of cyber attacks can be considered: (i) False Data Injection (FDI),
and (ii) Denials of Service (DoS). A FDI (a.k.a. deception attack) affects the
data integrity of packets by modifying their payloads [1, 2, 3]. A DoS is the
one that the attacker needs only to disrupt the system by preventing
communication between the components. In this paper, we focus on a specific
type of DoS attack, so-called Prevented Actuation Attack (PA2) [4, 5], where
the attacker prevents the exchange of information between the controller and
the actuators. An attacker can launch such attacks on the physical layer or
cyber layer. Examples of real-world PA2 are: sleep deprivation torture attack
[6] (a.k.a. battery exhaustion attack) that exhausts the battery of a
surveillance robot or a medical implant until it can no longer function; door
lock attack [7] that suppresses the operation of a smart door by injecting
‘close’ command every time an ‘open’ command is received; and fatigue bearing
attack [8] that restrains the operation of the lubricant system in wind
turbines to damage gearboxess.

Regardless of the type of attack, attack detection approaches presented in the
literature can be classified as: (i) passive approaches, and (ii) active
approaches. Note that we use the same terminology of the fault literature [9]
to classify attack detection approaches, as faults and attacks usually
manifest themselves similarly in control systems despite their natural
differences. In passive approaches, the input-output data of the system are
measured (remotely or on-site), analyzed for any possible stealthy behavior,
and then a decision about an attack is made. The passive approaches are widely
studied and commonly used in many today’s applications, e.g., [10, 11, 12, 13,
14]. However, they might not be able to recognize an attack when the input-
output data are not informative enough. Also, they do not address
stability/safety of the system during detection horizon, a time interval from
the instant an attack occurs to the instant when it is detected.

The active approaches interact with the system during the detection horizon by
means of a suitably designed input signal that is injected into the system to
increase the quality of detection, shorten the detection horizon, and enforce
stability/safety during the detection horizon. Contrary to the passive
approaches, the active approaches are historically younger and still under
development. To the best of the authors’ knowledge, the only existing active
attack detection approach in the literature is the physical authentication
(a.k.a. digital watermarking) [15, 16, 17]. The core idea of this method is to
inject a known noisy input to the system and observe its effect on the output
of the system. Thus, if an attacker is unaware of this physical watermark, the
system cannot be adequately emulated, as the attacker is unable to
consistently generate the component of the output associated with this known
noisy input. The physical authentication, which is mainly used in detection of
replay attack (a.k.a. playback attack) [18], can be effective if the noise
injected at the system input is large enough to achieve good detection
performance, which may degrade the control performance. Moreover, this method
injects the noisy input irrespective of the probability of attack occurrence,
which leads to unneeded loss in control performance. Furthermore, in the case
of constrained systems [19, 20, 21], as shown in [22], the extra uncertainty
injected to the system due to the noisy input should be taken into account in
the design procedure, which leads to tighter constraints, and consequently
more conservative behavior.

This paper answers the following question: How to determine the control input
sequence for a constrained CPS such as to improve the detection performance
without degrading the control performance? Inspired by [23], this paper
answers this question in the case of PA2. The proposed structure consists of
two units: (i) detection unit, and (ii) control unit. The detection unit uses
a priori information and input-output data of the system over the detection
horizon with a certain length to generate a decision variable which represents
the situation of the system. More precisely, the detection unit recognizes the
existence/inexistence of PA2 and distinguishes attacked actuators. The control
unit generates the control input which is optimal according to a cost function
and guarantees constraint satisfaction at all times. Both control and
detection aims are defined in the form of stochastic objective functions,
i.e., they are uncertain due to noises and initial condition. The open-loop
information processing strategy [24] is then used to express the stochastic
objective functions as deterministic functions. Finally, in order to evaluate
the quality of the control input sequence in terms of detection and control
aims, a compromise between the two aims is defined in the form of a multi-
objective optimization problem whose solution can be computed by means of
available optimization tools.

##  II Problem Statement

Consider the following discrete-time LTI system:

$\displaystyle x_{k+1}=$ $\displaystyle Ax_{k}+Bu_{k}+w_{k},$ (1)
$\displaystyle y_{k}=$ $\displaystyle Cx_{k}+v_{k},$ (2)

where $x_{k}\in\mathbb{R}^{n}$ is the state vector at time $k$,
$u_{k}\in\mathbb{R}^{p}$ is the control input at time $k$,
$y_{k}\in\mathbb{R}^{m}$ is the vector of measurements from the sensors at
time $k$, and the process noise $w_{k}\in\mathbb{R}^{n}$ and the measurement
noise $v_{k}\in\mathbb{R}^{m}$ are mutually independent white Gaussian noises
with zero mean and covariance matrices $H_{w}\in\mathbb{R}^{n\times n}$ and
$H_{v}\in\mathbb{R}^{m\times m}$, respectively. We assume that the initial
state $x_{0}$ is independent of $w_{k}$ and $v_{k}$, and has a Gaussian
distribution with the known mean $\bar{x}_{0}$ and covariance matrix
$H_{x,(0,0)}$.

In mathematical terms, the PA2 on the $i$-th actuator is equivalent to zeroing
the $i$-th column in the matrix $B$. Thus, the dynamics of attack-free and
under attack systems can be expressed using a single difference equation in
the following form:

$\displaystyle x_{k+1|\mu}=$ $\displaystyle Ax_{k|\mu}+B_{\mu}u_{k}+w_{k},$
(3) $\displaystyle y_{k|\mu}=$ $\displaystyle Cx_{k|\mu}+v_{k},$ (4)

where $\mu\in\\{\mu_{1},\cdots,\mu_{2^{p}}\\}$ is the index of the mode of the
system, each $\mu_{i}$ having a known distribution $P(\mu_{i})$, $B_{\mu}$ is
the corresponding input matrix, $x_{k|\mu}$ is the state of the system
operating in mode $\mu$ with $x_{0|\mu}=x_{0}$, and $y_{k|\mu}$ is the output
of the system operating in mode $\mu$.

Suppose that the system is subject to the following expectational linear
constraints:

$\text{E}\left[G_{x}x_{k|\mu_{i}}+G_{u}u_{k}\right]\leq g,~{}\forall k\geq
0,~{}i\in\\{1,\cdots,2^{p}\\}$ (5)

where $\text{E}[\cdot]$ is the expectation function, and
$G_{x}\in\mathbb{R}^{n_{c}\times n}$, $G_{u}\in\mathbb{R}^{n_{c}\times p}$,
and $g\in\mathbb{R}^{n_{c}}$, with $n_{c}$ as the number of constraints.

###### Problem 1

Consider system (3)-(4) which is subject to constraints (5). Suppose $N>0$ as
the detection horizon, chosen by the designer. Find a control sequence
$u_{k},~{}k=0,\cdots,N-1$ such that111Without loss of generality and for the
sake of simplicity we assume the detection horizon starts from 0. at time $N$
the mode of the system is identified with high probability of correctness,
while optimal control performance and constraint satisfaction are guaranteed
during the detection horizon.

Before starting with the solution of Problem 1, let us compute the conditional
probability density functions of the state and output. According to (3), for
the control sequence $u_{0},\cdots,u_{N-1}$, the mean value of the state at
time $k$ is

$\displaystyle\bar{x}_{k|\mu}=A^{k}\bar{x}_{0|\mu}+A^{k-1}B_{\mu}u_{0}+\cdots+B_{\mu}u_{k-1},$
(6)

and the covariance matrix of the state at times $k$ and $l$ ($k\geq l$) can be
computed as

$\displaystyle H_{x,(k,l)|\mu}:=$
$\displaystyle\text{E}\left\\{(x_{k|\mu}-\bar{x}_{k|\mu})(x_{l|\mu}-\bar{x}_{l|\mu})^{T}\right\\}$
$\displaystyle=$ $\displaystyle
A^{k}H_{x,(0,0)}\left(A^{l}\right)^{T}+A^{k-1}H_{w}\left(A^{l-1}\right)^{T}$
$\displaystyle+A^{k-2}H_{w}\left(A^{l-2}\right)^{T}+\cdots+A^{k-l}H_{w},$ (7)

where $H_{x,(k,l)|\mu}=\left(H_{x,(l,k)|\mu}\right)^{T}\in\mathbb{R}^{n\times
n}$. Therefore, the conditional probability density function of the state for
the interval $[0,N]$ can be expressed as:

$\displaystyle
P\left(x_{0:N}\big{|}\mu,u_{0:N-1}\right)\sim\mathcal{N}\left(\bar{x}_{0:N|\mu},H_{x|\mu}\right),$
(8)

where $x_{0:N}:=[(x_{0})^{T},\cdots,(x_{N})^{T}]^{T}\in\mathbb{R}^{n(N+1)}$,
$u_{0:N-1}:=[(u_{0})^{T},\cdots,(u_{N-1})^{T}]^{T}\in\mathbb{R}^{pN}$,
$\bar{x}_{0:N|\mu}:=[(\bar{x}_{0|\mu})^{T},\cdots,(\bar{x}_{N|\mu})^{T}]^{T}\in\mathbb{R}^{n(N+1)}$,
and $H_{x,(i-1,j-1)|\mu}$ is the element $(i.j)$ of
$H_{x|\mu}\in\mathbb{R}^{n(N+1)\times n(N+1)}$.

Similarly, the conditional probability density function of the output for the
the interval $[0,N]$ can be expressed as

$\displaystyle
P\left(y_{0:N}\big{|}\mu,u_{0:N-1}\right)\sim\mathcal{N}\left(\bar{y}_{0:N|\mu},H_{y|\mu}\right),$
(9)

where $y_{0:N}:=[(y_{0})^{T},\cdots,(y_{N})^{T}]^{T}\in\mathbb{R}^{m(N+1)}$,
and
$\bar{y}_{0:N|\mu}:=[(\bar{y}_{0|\mu})^{T},\cdots,(\bar{y}_{N|\mu})^{T}]^{T}\in\mathbb{R}^{m(N+1)}$
with $\bar{y}_{k|\mu}=C\bar{x}_{k|\mu},~{}k=0,\cdots,N$ as the mean value of
the output at time $k$. Also, $H_{y|\mu}\in\mathbb{R}^{m(N+1)\times m(N+1)}$
is the covariance matrix of the output, where the $(i,j)$ element is

$\displaystyle
H_{y,(i,j)|\mu}=\left\\{\begin{array}[]{rl}CH_{x,(i-1,j-1)|\mu}C^{T},&i>j\\\
CH_{x,(i-1,j-1)|\mu}C^{T}+H_{v}&i=j\end{array}\right..$ (12)

##  III Detection Unit

The system (3)-(4) can be seen as a $2^{p}$-model system, where each model
corresponds to one mode. Thus, in order to detect existence/inexistence of
PA2, and to identify which actuator(s) is under attack, it is only needed to
identify the true $\mu$ at the end of detection horizon. The fact that one of
the $2^{p}$ models is the true one can be modeled by a hypothesis random
variable that must belong to a discrete set of hypothesis
$\\{\mu_{1},\cdots,\mu_{2^{p}}\\}$, where the event $\mu_{i}$ means that the
$i$-th model is the one that is generating the data.

One Bayesian approach to hypothesis testing is to base decisions on the
posterior probabilities, i.e., the probability of the mode $\mu_{i}$
conditioned by the input-output data. In mathematical terms, the posterior
probabilities at time $k$ are denoted as
$P\left(\mu_{i}\big{|}y_{0:k},u_{0:k-1}\right)$, where at time $k=0$ the
posterior probabilities are equal to the prior probabilities, i.e.,
$P\left(\mu_{i}\big{|}u_{0}^{T}\right)=P(\mu_{i})$.

The conditioned posterior probabilities can be computed using the Multiple-
Model Adaptive Estimator (MMAE) structure [25, 26, 27]. The MMAE (a.k.a.
partitioned algorithm) involves the parallel operation of $2^{p}$ Kalman
filters (each matched to one of the postulated models), where the residuals of
the Kalman filters are used to compute the conditional posterior
probabilities. The rationale is that the highest posterior probability
corresponds to the true model of the system. It is shown that the correct
model can be identified “almost surely” [28, 29].

It is easy to show that when the attack happens sometime within a detection
horizon, it might remain undetected until the end of the following detection
horizon. It can be also shown that in the case of a smart attack (i.e., the
attack happens sometime within a detection horizon and lasts for a wisely
selected period of time), a single detector might not be adequate to detect
the attack. Therefore, we propose to deploy $N$ parallel detectors, where the
$d$-th detector identifies the mode $\hat{\mu}^{d}$ via

$\hat{\mu}^{d}=\arg\max\limits_{i\in\\{1,2,\cdots,2^{p}\\}}P\left(\mu_{i}\big{|}y_{0:N},u_{0:N-1}\right).$
(13)

We assume that every detector identifies the mode of the system only in $N$
time steps. Note that $N$ defines the the trade-off between the detection
quality and detection performance. Large values of $N$ decreases the
probability of making an incorrect decision during the transient of the
posterior probabilities. However, when the attack duration is too small
compared to the length of the detection horizon, the attack might remain
undetected.

##  IV Control Unit

###  IV-A Control Objective Function

The control aim is to track the desired reference $r_{k}\in\mathbb{R}^{m}$
while penalizing the control effort. The control objective function can be
formulated as

$\displaystyle J_{c}(u_{0:N-1})=\text{E}\Bigg{[}$
$\displaystyle\sum\limits_{k=0}^{N}\left\|y_{k}-r_{k}\right\|^{2}_{Q}+\sum\limits_{k=0}^{N-1}\left\|u_{k}\right\|^{2}_{R}\Bigg{]}$
(14)

where $Q=Q^{T}\in\mathbb{R}^{m\times m}$ is a positive semi-definite matrix
and $R=R^{T}\in\mathbb{R}^{p\times p}$ is a positive definite matrix.

$\displaystyle F_{1}=$
$\displaystyle\left[\begin{matrix}B_{\mu_{i}}^{T}(A^{k-1})^{T}C^{T}QCA^{k-1}B_{\mu_{i}}&\cdots&B_{\mu_{i}}^{T}(A^{k-1})^{T}C^{T}QCB_{\mu_{i}}\\\
\vdots&\ddots&\vdots\\\
(B_{\mu_{i}})^{T}C^{T}QCA^{k-1}B_{\mu_{i}}&\cdots&(B_{\mu_{i}})^{T}C^{T}QCB_{\mu_{i}}\end{matrix}\right],$
(18) $\displaystyle F_{2}=$ $\displaystyle
2(\bar{x}_{0|\mu_{i}})^{T}(A^{k})^{T}Q\left[\begin{matrix}A^{k-1}B_{\mu_{i}}&A^{k-2}B_{\mu_{i}}&\cdots&B_{\mu_{i}}\end{matrix}\right]-2r_{k}^{T}QC\left[\begin{matrix}A^{k-1}B_{\mu_{i}}&A^{k-2}B_{\mu_{i}}&\cdots&B_{\mu_{i}}\end{matrix}\right],$
(19) $\displaystyle F_{3}=$
$\displaystyle\sum\limits_{k=0}^{N}\sum\limits_{i=0}^{2^{p}}P(\mu_{i})(\bar{x}_{0|\mu_{i}})^{T}(A^{k})^{T}C^{T}QCA^{k}\bar{x}_{0|\mu_{i}}+\text{Tr}\left(Q\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})H_{y,(k,k)|\mu_{i}}\right)+\sum\limits_{k=0}^{N}r_{k}^{T}Qr_{k}$
$\displaystyle+\sum\limits_{k=0}^{N}\sum\limits_{i=0}^{2^{p}}P(\mu_{i})r^{T}QCA^{k}\bar{x}_{0|\mu_{i}}.$
(20)

The objective function (14) is a stochastic function, where the uncertainties
are due to noises and the initial condition. In this paper, instead of using
deterministic approaches (i.e., assuming uncertainties as upper-bounded
signals), we will focus on probabilistic approaches, where available a priori
information can be used in obtaining the optimal control sequence. In
particular, we will use the open loop approach. This approach is based on the
information available at the beginning of each detection horizon (i.e.,
$\bar{x}_{0}$ and $H_{x,(0,0)}$), while the measurements during the detection
horizon are not used.

###### Theorem 1

Consider system (3)-(4), and control objective function (14). Suppose that the
open loop approach is used to determine the control sequence, i.e., the entire
control sequence $u_{0:N-1}$ is determined at the beginning of the prediction
horizon. Then, control objective function (14) can be expressed as an explicit
function of the control sequence.

###### Proof:

When the open loop approach is used, the objective function (14) can be
expressed as

$\displaystyle J_{c}(\cdot)=$
$\displaystyle\text{Tr}\Bigg{(}Q\sum\limits_{k=0}^{N}\text{E}\Big{[}y_{k}y_{k}^{T}\Big{|}u_{0:N-1}\Big{]}\Bigg{)}+\sum\limits_{k=0}^{N}r_{k}^{T}Qr_{k}$
$\displaystyle-2\sum\limits_{k=0}^{N}r_{k}^{T}Q\text{E}\Big{[}y_{k}\Big{|}u_{0:N-1}\Big{]}+\sum\limits_{k=0}^{N-1}u_{k}^{T}Ru_{k},$
(13)

where $\text{Tr}(\cdot)$ is the trace function. We know
that222$\text{Cov}(Y,Y)=\text{E}\\{YY^{T}\\}-\text{E}\\{Y\\}\left(\text{E}\\{Y\\}\right)^{T}$
for the random vector $Y$.

$\displaystyle\text{E}\Big{[}y_{k}y_{k}^{T}\big{|}u_{0:N-1}\Big{]}$
$\displaystyle=\sum\limits_{i=1}^{2^{p}}P(\mu_{i})\left(\bar{y}_{k|\mu_{i}}\bar{y}_{k|\mu_{i}}^{T}+H_{y,(k,k)|\mu_{i}}\right),$
(14)

where $\text{Cov}(\cdot)$ is the covariance function. Thus, the control
objective function can be rewritten as

$\displaystyle J_{c}(\cdot)=$
$\displaystyle\text{Tr}\left(Q\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})\left(\bar{y}_{k|\mu_{i}}\bar{y}_{k|\mu_{i}}^{T}+H_{y,(k,k)|\mu_{i}}\right)\right)$
$\displaystyle+\sum\limits_{k=0}^{N-1}u_{k}^{T}Ru_{k}+\sum\limits_{k=0}^{N}r_{k}^{T}Qr_{k}$
$\displaystyle-2\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})r_{k}^{T}Q\bar{y}_{k|\mu_{i}},$
(15)

which due to the fact that $\bar{y}_{k|\mu}=C\bar{x}_{k|\mu}$, it implies
that:

$\displaystyle J_{c}(\cdot)=$
$\displaystyle\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})\bar{x}_{k|\mu_{i}}^{T}C^{T}QC\bar{x}_{k|\mu_{i}}+\sum\limits_{k=0}^{N-1}u_{k}^{T}Ru_{k}$
$\displaystyle-2\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})r_{k}^{T}QC\bar{x}_{k|\mu_{i}}+\sum\limits_{k=0}^{N}r_{k}^{T}Qr_{k}$
$\displaystyle+\text{Tr}\left(Q\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})H_{y,(k,k)|\mu_{i}}\right).$
(16)

Finally, according to (6), (IV-A) can be rewritten as

$\displaystyle J_{c}(\cdot)=$
$\displaystyle\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})u_{0:k-1}^{T}F_{1}u_{0:k-1}+\sum\limits_{k=0}^{N-1}u_{k}^{T}Ru_{k}$
$\displaystyle+\sum\limits_{k=0}^{N}\sum\limits_{i=1}^{2^{p}}P(\mu_{i})F_{2}u_{0:k-1}+F_{3},$
(17)

where $F_{1}$, $F_{2}$, and $F_{3}$ are given in (18)-(20), respectively. This
completes the proof. ∎

###  IV-B Detection Objective Function

Suppose that the detector given in (13) is used to identify the mode of the
system. In order to determine the control sequence during the detection
horizon such that the probability of an incorrect identification is minimized,
we can use the following detection objective function

$J_{d}(u_{0:N})\triangleq\text{E}\big{[}\sigma(\hat{\mu}^{1})\big{]},$ (21)

where $\sigma(\hat{\mu}^{1})$ is zero when the identified mode is the actual
mode of the system, and is non-zero (we set it to 1 for simplicity) otherwise.
Note that since the control sequence $u_{0:N-1}$ is assumed to be applied at
time $k=0$, only the first detection horizon is taken into account in the
formulation of the detection objective function. In other words, the
determined control signal is optimal only for Detector#1.

###### Theorem 2

Consider system (3)-(4), and detection objective function (21). Suppose that
the open loop approach is used to determine the control sequence. Then, the
detection objective function can be upper bounded with an explicit function of
the control sequence.

###### Proof:

By using the open loop approach, the detection objective function (21) can be
expressed as

$\displaystyle J_{d}(\cdot)=$
$\displaystyle\text{E}\big{[}\sigma(\hat{\mu}^{1})\big{|}u_{0:N-1}\big{]}$
$\displaystyle=$
$\displaystyle\int\limits_{\mathbb{R}^{m(N+1)}}\sum\limits_{i=1}^{2^{p}}\sigma(\hat{\mu}^{1})P\left(\mu_{i}\big{|}y_{0:N},u_{0:N-1}\right)\cdot$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}P(y_{0:N}\big{|}u_{0:N-1})dy_{0:N},$
(22)

which according to Bayes’ theorem, it implies that

$\displaystyle J_{d}(\cdot)=$
$\displaystyle\int\limits_{\mathbb{R}^{m(N+1)}}\sum\limits_{i=1}^{2^{p}}\sigma(\hat{\mu}^{1})P(y_{0:N}\big{|}\mu_{i},u_{0:N-1})P\left(\mu_{i}\right)dy_{0:N},$
(23)

which is concluded due to the
fact333$P(\mu_{i}|u_{0:N-1})=\frac{P(u_{0:N-1}|\mu_{i})}{P(u_{0:N-1})}P(\mu_{i})=P(\mu_{i})$,
since $u_{0:N-1}$ is deterministic, and consequently
$P(u_{0:N-1}|\mu_{i})=P(u_{0:N-1})=1$. that the probability of the mode
$\mu_{i}$ conditioned by only input data is equal to the probability of the
mode $\mu_{i}$.

The right side of (23) cannot be computed analytically and its numerical
evaluation is computationally expensive. Due to this reason, in the following
we will find an upper bound for the detection objective function
$J_{d}(u_{0},\cdots,u_{N-1})$.

Since $0\leq\sigma(\hat{\mu}^{1})\leq 1$, it implies that:

$\displaystyle J_{d}(\cdot)\leq$
$\displaystyle\int\limits_{\mathbb{R}^{m(N+1)}}\sum\limits_{i=1}^{2^{p}}P(y_{0:N}\big{|}\mu_{i},u_{0:N-1})P\left(\mu_{i}\right)dy_{0:N}.$
(24)

Following the same arguments presented in [30], the right side of (24) can be
upper bounded as

$\displaystyle\int\limits_{\mathbb{R}^{m(N+1)}}\sum\limits_{\mu=1}^{2^{p}}P(y_{0:N}\big{|}\mu_{i},u_{0:N-1})P\left(\mu_{i}\right)dy_{0:N}\leq\hat{J}_{d}(u_{0:N}),$
(25)

where

$\displaystyle\hat{J}_{d}(u_{0:N})=\sum\limits_{i=1}^{2^{p}}\sum\limits_{j=i+1}^{2^{p}}\sqrt{P(\mu_{i})P(\mu_{j})}e^{-\phi_{ij}},$
(26)

with

$\displaystyle\phi_{ij}$
$\displaystyle=\frac{1}{4}\left(\bar{y}_{0:N|\mu_{j}}-\bar{y}_{0:N|\mu_{i}}\right)^{T}\left(H_{y|\mu_{i}}+H_{y|\mu_{j}}\right)^{-1}\cdot$
$\displaystyle\left(\bar{y}_{0:N|\mu_{j}}-\bar{y}_{0:N|\mu_{i}}\right)+\frac{1}{2}\ln\left(\frac{\text{det}\left(\frac{H_{y|\mu_{i}}+H_{y|\mu_{j}}}{2}\right)}{\sqrt{\text{det}(H_{y|\mu_{i}})\text{det}(H_{y|\mu_{j}})}}\right)$
(27)

where $\text{det}(\cdot)$ is the determinant function. It is noteworthy that
according to (6)-(7), (12), and since
$\bar{y}_{k|\mu_{i}}=C\bar{x}_{k|\mu_{i}}$, the upper bound
$\hat{J}_{d}(\cdot)$ given in (26) is an explicit function of the control
sequence. This completes the proof. ∎

###  IV-C Constraints

By using the open loop approach, the expectational constraints given in (5)
take the following form:

$\displaystyle\text{E}\left[G_{x}x_{k|\mu_{i}}+G_{u}u_{k}\big{|}u_{0:k-1}\right]\leq
g\Rightarrow$ $\displaystyle
G_{x}(A^{k}\bar{x}_{0|\mu_{i}}+A^{k-1}B_{\mu_{i}}u_{0}+\cdots+B_{\mu_{i}}u_{k-1})$
$\displaystyle+G_{u}u_{k}\leq g,$ (28)

which is an explicit function of the control sequence.

###  IV-D Proposed Solution

One possible way to pursue both control and detection aims is to let one of
the objective functions to take arbitrary value up to a known upper limit
value, and then to enforce this as a constraint and minimize the other
objective function. Therefore, the following two optimization problems can be
considered

$u_{0:N-1}^{\ast}=\left\\{\begin{array}[]{ll}&\arg\,\min\limits_{u_{0:N-1}}J_{c}(\cdot)\text{
given in \eqref{eq:ControlOFFinal}}\\\
\text{s.t.}&\eqref{eq:constraintsFinal}\text{ is satisfied }\forall
i,~{}\forall k\geq 0\\\ &\hat{J}_{d}(\cdot)\text{ given in
\eqref{eq:Jdhat}}\leq\bar{J}_{d}\end{array}\right.,$ (29)

or

$u_{0:N-1}^{\ast}=\left\\{\begin{array}[]{ll}&\arg\,\min\limits_{u_{0:N-1}}\hat{J}_{d}(\cdot)\text{
given in \eqref{eq:Jdhat}}\\\ \text{s.t.}&\eqref{eq:constraintsFinal}\text{ is
satisfied }\forall i,~{}\forall k\geq 0\\\ &J_{c}(\cdot)\text{ given in
\eqref{eq:ControlOFFinal}}\leq\bar{J}_{c}\end{array}\right.,$ (30)

where $\bar{J}_{d}$ and $\bar{J}_{c}$ are maximum acceptable levels of the
detection and control objective functions, respectively.

The objective function (IV-A) and constraints given in (28) are convex in
$u_{0:N-1}$. The objective function (26) is concave, as $\phi_{ij}$ as in (27)
is a quadratic function of $u_{0:N-1}$ (with a positive definite matrix), and
consequently convex in $u_{0:N-1}$. Thus, problems (29)-(30) are in general
non-convex. We use bmibnb [31] to numerically compute their solutions.

##  V Simulation Study– Irrigation Channel

In this section we will use the developed method to control the level of water
in pools 9 and 10 of the Haughton main channel, as shown in Fig. 1. The water
levels in the channel are controlled by overshot gates located along the
channel. The stretch of a channel between two gates is referred to as a reach
or a pool. We assume that the communication between the controller and the
gates is through internet.

The water level in the $g$-th pool ($g\in\\{9,10\\}$) of the irrigation
channel can be modeled as [32]

$\displaystyle\dot{y}_{g}(t)=\alpha_{g-1,in}h_{g-1}^{3/2}(t-\tau_{g-1})-\alpha_{g,out}h_{g}^{3/2}(t)+d_{g-1}(t),$
(31)

where $y_{g}(t)$ is the water level in the pool, $h_{g}(t)$ is the head over
the gate (the height of water above the gate), $\tau_{g}$ is the time delay
which accounts for the time it takes for the water to travel from the upstream
gate to the downstream gate in the $g$-th pool, $d_{g}(t)$ represents offtakes
to farms and side channels, and $\alpha_{g,in}$ and $\alpha_{g,out}$ are
constants which incorporates the effect of the discharge coefficients. The
real value of the parameters is given in TABLE I. For the sake of simplicity
we assume there is no offtake, i.e., $d_{g}=0,~{}g=8,9$.

The sampling time is 10 [min], and the control signal is the head over the
gate. We assume that initial water level in pool 9 and 10 is 6.60 [m] and 5.60
[m], respectively. Also, we assume that
$w_{k}\sim\mathcal{N}\left(\textbf{0},0.3I_{8}\right)$ and
$v_{k}\sim\mathcal{N}\left(\textbf{0},0.3I_{2}\right)$, where $I_{2}$ is the
$2\times 2$ identity matrix and 0 is the zero vector with appropriate size.
The water level in pools should not exceed 15 [m]. The system is subject to
actuator saturation [33], i.e., the control signals cannot be negative.

Since there are three actuators, eight different modes can be defined. We
assume that a priori probability of the $i$-th mode is
$P(\mu_{i})=0.125,~{}\forall i$.

Suppose that $Q=I_{2}$, $R=I_{3}$, the detection horizon is 200 [min], and the
level of detection and control objective functions must not exceed $1$ and
$2000$, respectively.

Figure 1: Side view of the Haughton main channel; pools 9 and 10. TABLE I: Parameters of the Haughton Main Channel [34]. Parameter | $g=8$ | $g=9$ | $g=10$
---|---|---|---
$\alpha_{g,in}$ [1/m2] | 0.0208 | 0.0700 | 0.0142
$\alpha_{g,out}$ [1/m2] | 0.0278 | 0.0614 | 0.0156
$\tau_{g}$ [min] | 6 | 3 | 16

We assume that for $k\in[0,80]$, $[200,300]$, $[360,480]$, $[580,700]$ the
mode of the system is 1, for $k\in[80,200]$ the mode of the system is 8, for
$k\in[300,360]$ the mode of the system is 2, and for $k\in[480,580]$ the mode
of the system is 7. Note that for comparison purposes, we also simulate a pure
control formulation, i.e.,

$u_{0:N-1}^{\ast}=\left\\{\begin{array}[]{ll}&\arg\,\min\limits_{u_{0:N-1}}\,J_{c}(\cdot)\text{
given in \eqref{eq:ControlOFFinal}}\\\
\text{s.t.}&\eqref{eq:constraintsFinal}\text{ is satisfied }\forall
i,~{}\forall k\geq 0\end{array}\right..$ (33)

The achieved normalized values of the control and detection objection
functions are shown in Fig. 2 and 3, where control and detection costs
obtained by the formulation (33) are assumed as the base unit quantity for
control and detection costs, respectively. As expected, compromise between
control and detection aims increases the control cost $J_{c}$. However, it
decreases the detection cost $\hat{J}_{d}$ which means that probability of
misidentification is minimized.

##  VI Conclusion

This paper proposed an optimization approach for active attack detection and
control of constrained CPSs systems subject to expectational linear
constraints. This paper mainly focused on PA2 attack, where the attacker
prevents the exchange of information between the controller and the actuators.
A set of parallel detectors based on hypothesis testing approach was proposed.
Using a probabilistic approach to deal with uncertainties, the detection and
control aims were formulated as two separate stochastic objective functions.
The open loop approach was deployed to transfer the stochastic functions to
deterministic ones. Two alternative compromise between detection and control
aims were presented in the form of a constrained optimization problem. The
effectiveness of the proposed active approach was validated through simulation
studies.

Figure 2: Normalized control cost by formulations (33), (29), and (30). Figure
3: Normalized detection cost by formulations (33), (29), and (30).

## References

  * [1] Y. Mo, J. P. Hespanha, and B. Sinopoli, “Resilient detection in the presence of integrity attacks,” _IEEE Trans. Signal Process._ , vol. 62, no. 1, pp. 31–43, Jan. 2014.
  * [2] Y. Mo, E. Garone, A. Casavola, and B. Sinopoli, “False data injection attacks against state estimation in wireless sensor networks,” in _Proc. 49th IEEE Conf. Decision and Control_ , Hilton Atlanta Hotel, Atlanta, GA, USA, Dec. 15-17, 2010, pp. 5967–5972.
  * [3] C.-Z. Bai, F. Pasqualetti, and V. Gupta, “Data-injection attacks in stochastic control systems: Detectability and performance tradeoffs,” _Automatica_ , vol. 82, pp. 251–260, Aug. 2017.
  * [4] G. Loukas, _Cyber-Physical Attacks: A Growing Invisible Threat_. Butterworth-Heinemann, 2015.
  * [5] L. K. Carvalho, Y.-C. Wu, R. Kwong, and S. Lafortune, “Detection and prevention of actuator enablement attacks in supervisory control systems,” in _Proc. 13th Int. Workshop Discrete Event Systems_ , Xi’an, China, May 30-Jun. 1, 2016, pp. 298–305.
  * [6] F. Stajano and R. Anderson, “The resurrecting duckling: Security issues for ad-hoc wireless networks,” in _Proc. Cambridge Int. Workshop Security Protocols_ , Cambridge, United Kingdom, Apr. 19-21, 1999, pp. 172–182.
  * [7] Y. Ha, S.-H. Jang, K.-W. Kim, and J. W. Yoon, “Side channel attack on digital door lock with vibration signal analysis: Longer password does not guarantee higher security level,” in _Proc. IEEE Int. Conf. Multisensor Fusion and Integration for Intelligent Systems_ , Daegu, South Korea, Nov. 16-18, 2017, pp. 103–110.
  * [8] H. Wu, J. Liu, J. Liu, M. Cui, X. Liu, and H. Gao, “Power grid reliability evaluation considering wind farm cyber security and ramping events,” _Appl. Sci._ , vol. 9, no. 15, Jul. 2019.
  * [9] M. Hosseinzadeh and F. R. Salmasi, “Determination of maximum solar power under shading and converter faults- a prerequisite for failure-tolerant power management systems,” _Simul. Model Pract. Theory_ , vol. 62, pp. 14–30, Mar. 2016.
  * [10] Y. Mo, R. Chabukswar, and B. Sinopoli, “Detecting integrity attacks on SCADA systems,” _IEEE Trans. Control Syst. Technol._ , vol. 22, no. 4, pp. 1396–1407, Jul. 2014.
  * [11] H. Zhang, P. Cheng, L. Shi, and J. Chen, “Optimal DoS attack scheduling in wireless networked control system,” _IEEE Trans. Control Syst. Technol._ , vol. 24, no. 3, pp. 843–852, May 2016.
  * [12] Y. Li, D. E. Quevedo, S. Dey, and L. Shi, “SINR-based DoS attack on remote state estimation: A game-theoretic approach,” _IEEE Trans. Control Netw. Syst._ , no. 3, pp. 632–642, Sep. 2017.
  * [13] Y. Li, D. Shi, and T. Chen, “False data injection attacks on networked control systems: A stackelberg game analysis,” _IEEE Trans. Autom. Control_ , vol. 63, no. 10, pp. 3503–3509, Oct. 2018.
  * [14] R. Zhang and P. Venkitasubramaniam, “False data injection and detection in lqg systems: A game theoretic approach,” _IEEE Trans. Netw. Syst._ , 2019.
  * [15] Y. Mo, S. Weerakkody, and B. Sinopoli, “Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs,” _IEEE Control Syst. Mag._ , vol. 35, no. 1, pp. 93–109, Feb. 2015\.
  * [16] T. Irita and T. Namerikawa, “Detection of replay attack on smart grid with code signal and bargaining game,” in _Proc. 2017 American Control Conf._ , Seattle, WA, USA, May 24-26, 2017, pp. 2112–2117.
  * [17] H. S. Sánchez, D. Rotondo, T. Escobet, V. Puig, J. Saludes, and J. Quevedo, “Detection of replay attacks in cyber-physical systems using a frequency-based signature,” _J. Franklin Inst._ , vol. 356, no. 5, pp. 2798–2824, Mar. 2019.
  * [18] Y. Mo and B. Sinopoli, “Secure control against replay attacks,” in _Proc. 47th Annu. Allerton Conf. Communication, Control, and Computing_ , Allerton House, Illinois, USA, Sep. 30-Oct. 2, 2009, pp. 911–918.
  * [19] M. Hosseinzadeh, A. Cotorruelo, D. Limon, and E. Garone, “Constrained control of linear systems subject to combinations of intersections and unions of concave constraints,” _IEEE Control Syst. Lett_ , vol. 3, no. 3, pp. 571–576, Jul. 2019.
  * [20] M. Hosseinzadeh, B. Sinopoli, and A. F. Bobick, “An explicit reference governor for time-varying linear constraints,” in _Proc. 59th IEEE Conf. Decision and Control_ , Jeju Island, Korea (South), Dec. 14-18, 2020, pp. 3323–3328.
  * [21] M. Hosseinzadeh and E. Garone, “An explicit reference governor for the intersection of concave constraints,” _IEEE Trans. Autom. Control_ , vol. 65, no. 1, pp. 1–11, Jan. 2020.
  * [22] M. Hosseinzadeh, B. Sinopoli, and E. Garone, “Feasibility and detection of replay attack in networked constrained cyber-physical systems,” in _Proc. 2019 57th Annual Allerton Conference on Communication, Control, and Computing_ , Allerton Park and Retreat Center,Monticello, IL, USA, Sep. 24-27, 2019, pp. 712–717.
  * [23] I. Punčochář, J. Široký, and M. Šimandl, “Constrained active fault detection and control,” _IEEE Trans. Autom. Control_ , vol. 60, no. 1, pp. 253–258, Jan. 2015.
  * [24] M. Šimandl and I. Punčochář, “Active fault detection and control: Unified formulation and optimal design,” _Automatica_ , vol. 45, no. 9, pp. 2052–2059, Sep. 2009.
  * [25] S. Fekri, M. Athans, and A. Pascoal, “RMMAC: A novel robust adaptive control scheme– Part I: Architecture,” in _Proc. 43rd IEEE Conf. Decision and Control_ , Atlantis, Paradise Island, Bahamas, Dec. 14-17, 2004, pp. 1134–1139.
  * [26] V. Hassani, A. P. Aguiar, M. Athans, and A. M. Pascoal, “Multiple model adaptive estimation and model identification using a minimum energy criterion,” in _Proc. 2009 American Control Conf._ , St. Louis, MO, USA, Jun. 10-12, 2009, pp. 518–523.
  * [27] N. Sadati, M. Hosseinzadeh, and G. A. Dumont, “Multi-model robust control of depth of hypnosis,” _Biomed. Signal Process._ , vol. 40, pp. 443–453, Feb. 2018.
  * [28] S. Fekri, M. Athans, and A. Pascoal, “Issues, progress and new results in robust adaptive control,” _Int. J. Adapt. Control Signal Process._ , vol. 20, no. 10, pp. 519–579, Dec. 2006.
  * [29] D. Rotondo, V. Hassani, and A. Cristofaro, “A multiple model adaptive architecture for the state estimation in discrete-time uncertain LPV systems,” in _Proc. 2017 American Control Conf._ , Sheraton Seattle Hotel, May 24-26, 2017, pp. 2393–2398.
  * [30] L. Blackmore and B. Williams, “Finite horizon control design for optimal discrimination between several models,” in _Proc. 45th IEEE Conf. Decision and Control_ , San Diego, CA, USA, Dec. 13-15, 2006, pp. 1147–1152.
  * [31] J. Lofberg, “YALMIP: A toolbox for modeling and optimization in MATLAB,” in _Proc. IEEE Int. Symp. Computer Aided Control Syst. Design_ , Taipei, Taiwan, Sep. 2-4, 2004, pp. 284–289.
  * [32] P. Zhang and E. Weyer, “A reference model approach to performance monitoring of control loops with applications to irrigation channels,” _Int. J. Adapt. Control Signal Process._ , vol. 19, no. 10, pp. 797–818, Dec. 2005.
  * [33] M. Hosseinzadeh and M. J. Yazdanpanah, “Robust adaptive passivity-based control of open-loop unstable affine nonlinear systems subject to actuator saturation,” _IET Control Theory A._ , vol. 11, no. 16, pp. 2731–2742, Nov. 2017.
  * [34] E. Weyer, “System identification of an open water channel,” _Control Eng. Pract._ , vol. 9, no. 12, pp. 1289–1299, Dec. 2001.


# Accelerating differential evolution algorithm with Gaussian sampling based
on estimating the convergence points

1st Rui Zhong  Graduate School of Information Science and Technology
Hokkaido University
Sapporo, Japan
2nd Masaharu Munetomo  Information Initiative Center
Hokkaido University
Sapporo, Japan

###### Abstract

In this paper, we propose a simple strategy for estimating the convergence
point approximately by averaging the elite sub-population. Based on this idea,
we derive two methods, which are ordinary averaging strategy, and weighted
averaging strategy. We also design a Gaussian sampling operator with the mean
of the estimated convergence point with a certain standard deviation. This
operator is combined with the traditional differential evolution algorithm
(DE) to accelerate the convergence. Numerical experiments show that our
proposal can accelerate the DE on most functions of 28 low-dimensional test
functions on the CEC2013 Suite, and our proposal can easily be extended to
combine with other population-based evolutionary algorithms with a simple
modification.

###### Index Terms:

Estimation of convergence point, Gaussian sampling, Acceleration, Averaging
strategy

##  I Introduction

As an important part of artificial intelligence, evolutionary computation has
achieved great success in solving continuous[1], large-scale[2],
constraint[3], and multi-objective optimization problems[4] in the past few
decades. However, according to the No Free Lunch Theorem (NFLT)[5], there is
no optimization algorithm that can solve all optimization problems perfectly.
NFLT proves that the average performance of any pair of algorithms $A$ and $B$
is identical on all possible problems. Therefore, if an algorithm performs
well on a certain class of problems, it must pay for that with performance
degradation on the remaining problems, since this is the only way for all
algorithms to have the same performance on average across all functions. This
hypothesis has a great impact on the scientific community for the notion that
there is no universal optimizer. In addition, a single algorithm often cannot
balance the exploration and the exploitation well. For several reasons, hybrid
algorithms or memetic algorithms appears, the hybrid algorithm was first
proposed to solve the Traveling Salesman Problem (TSP) by modifying the
genetic algorithm with a local search operator[6]. However, NFLT also limits
the hybrid algorithms with identical average performance on all possible
problems, thus the features of the problem become the starting point for
designing a suitable optimization algorithm.

Since hybrid algorithms were not proposed as specific optimization algorithms,
but as a broad class of algorithms inspired by the diffusion of ideas and
composed of multiple existing operators, the community started showing
increasing attention toward these algorithmic structures as a general
guideline for addressing specific problems.

In addition, many studies[7, 8] show that it is a promising strategy to find
trusting regions or elites in the fitness landscape, which can be fully
applied to guide the direction of evolution well. Murata et al[9]. first
proposed that a mathematical method could be used to calculate the global
optimum using the information of two subsequent generations. Yu et al[10].
proposed that in the differential evolution algorithm, the differential vector
from the parent individual to its offspring (or the vector from an individual
with poor fitness to an individual with higher fitness) is defined as the
moving vector, and the convergence point is estimated according to the moving
vector.

In this paper, we employ a simple strategy to approximately estimate the
convergence point, which is roughly estimated by the averaging strategy and
weighted averaging strategy of the elite sub-population. Besides, after each
generation optimized by DE, we apply a Gaussian sampling operator with the
mean of the estimated convergence point. In experiment design, we evaluate our
proposal with 28 benchmark functions from the CEC 2013 Suite[11]. Finally, we
summarize our work and provide some open topics for future discussions.

The rest of the paper is organized as follows. Section II includes
preliminaries and related works. Section III provides a detailed description
of our proposal. Section IV provides the numerical experiments and experiments
analysis. Finally, Section V concludes the paper and shows the future
directions.

##  II Preliminaries and Related works

In this section, we first introduce the preliminaries of this work, including
the Random Search (RS), Genetic Algorithm (GA), Differential Evolution
Algorithm (DE), Evolution Strategy (ES), and Particle Swarm Optimization
(PSO). Finally, we introduce the original estimation of a convergence point.

###  II-A Preliminaries

####  II-A1 RS

RS was proposed by Rastrigin in 1963, and an early introduction to RS with
basic mathematical analysis was given in the paper [12]. The principle of RS
is to iteratively move to better positions in the search space, these
positions are randomly sampled around the current optimum. RS is an
optimization method that does not require the gradients of the fitness
landscape, so RS can be employed for discrete or differentiable problems. The
Pseudocode of RS is shown in Algorithm 1.

Input: ${\rm Dimension}:D;{\rm Search\ space}:S;{\rm Generation}:T$

Output: ${\rm Best\ solution}:E$

Function  _RS(_$D,S,T$_)_:

$t\leftarrow 0$

$\blacktriangleright$ (Solution initialization)

$P_{t}\leftarrow\textbf{Initial}(D,S)$

$FP\leftarrow\textbf{Evaluate}(P_{t})$

$E\leftarrow P_{t}$

while  _$t <T$ and not convergence_ do

$O_{t}\leftarrow\textbf{Sampling}(P_{t})$

$FO\leftarrow\textbf{Evaluate}(O_{t})$

if  _$FO >FP$_ then

$P_{t}\leftarrow O_{t}$

$E\leftarrow O_{t}$

end if

$t\leftarrow t+1$

end while

$\textbf{return}\ E$

Algorithm 1 RS

####  II-A2 GA

GA[13] simulates the biological evolution process of chromosomes with
selection, crossover, and mutation. Chromosomes present problem solutions and
are evaluated based on the fitness function to select parents. Selection is an
important process for choosing parents to produce the offspring and can affect
the convergence of GA, Mutation increases the diversity of the population by
randomly modifying the genes in the chromosome with a certain probability. The
Pseudocode of GA is shown in Algorithm 2

Input: ${\rm Dimension}:D;{\rm Search\ space}:S;{\rm Population\ size}:PS;{\rm
Generation}:T$

Output: ${\rm Best\ solution}:E$

Function  _GA(_$D,S,PS,T$_)_:

$t\leftarrow 0$

$\blacktriangleright$ (Population initialization)

$P_{t}\leftarrow\textbf{Initial}(D,S,PS)$

$F_{t}\leftarrow\textbf{Evaluate}(P_{t})$

$E\leftarrow\textbf{bestIndividual}(P_{t},F_{t})$

while  _$t <T$ and not convergence_ do

$\blacktriangleright$ (Selection, Crossover and Mutation)

$O_{t}\leftarrow\textbf{Selection}(P_{t},F_{t})$

$O_{t}\leftarrow\textbf{Crossover}(O_{t})$

$P_{t+1}\leftarrow\textbf{Mutation}(O_{t})$

$F_{t+1}\leftarrow\textbf{Evaluate}(P_{t+1})$

$E\leftarrow\textbf{bestIndividual}(P_{t+1},F_{t+1})$

$t\leftarrow t+1$

end while

$\textbf{return}\ E$

Algorithm 2 GA

####  II-A3 DE

DE was first proposed by Storn and Price in 1995[14]. DE has been applied to
solve many complex optimization problems, and the procedure is similar to the
GA, DE mainly includes mutation, crossover, and selection. The Pseudocode of
DE is shown in Algorithm 3.

Input: ${\rm Dimension}:D;{\rm Search\ space}:S;{\rm Population\ size}:PS;{\rm
Generation}:T$

Output: ${\rm Best\ solution}:E$

Function  _DE(_$D,S,PS,T$_)_:

$t\leftarrow 0$

$\blacktriangleright$ (Population initialization)

$P_{t}\leftarrow\textbf{Initial}(D,S,PS)$

$F_{t}\leftarrow\textbf{Evaluate}(P_{t})$

$E\leftarrow\textbf{bestIndividual}(P_{t},F_{t})$

while  _$t <T$ and not convergence_ do

$\blacktriangleright$ (Selection, Crossover and Mutation)

$O_{t}\leftarrow\textbf{Mutation}(P_{t},F_{t})$

$O_{t}\leftarrow\textbf{Crossover}(O_{t})$

$F_{t+1}\leftarrow\textbf{Evaluate}(O_{t})$

$P_{t+1}\leftarrow\textbf{Selection}(O_{t},F_{t+1})$

$E\leftarrow\textbf{bestIndividual}(P_{t+1},F_{t+1})$

$t\leftarrow t+1$

end while

$\textbf{return}\ E$

Algorithm 3 DE

Mutation is an essential process in DE. We briefly introduce 2 common mutation
strategies:

$\displaystyle\textbf{DE/rand/1}:V_{i}(g)=X_{p1}(g)+F\cdot(X_{p2}(g)-X_{p3}(g))$
$\displaystyle\textbf{DE/best/1}:V_{i}(g)=X_{best}(g)+F\cdot(X_{p1}(g)-X_{p2}(g))$

$X_{pi}(g)$ is the different individuals randomly selected from the current
population, $X_{best}(g)$ is the best individual in the current population.
$F$ is the scaling factor. The differential vectors between individuals are
calculated in the mutation is the origin of the name of DE.

####  II-A4 ES

ES was proposed in 1963[15], as an optimization algorithm, ES imitates the
mechanism of biological evolution. Under the hypothesis that no matter what
changes occur in genes, the results or traits always follow a Gaussian
distribution with zero mean and a certain variance. Although ES employs
mutation and crossover to generate offspring which is similar to GA, ES
emphasizes the phenotype rather than the genotype, and ES also applies the
real coding instead of binary coding in GA. The Pseudocode of ES is shown in
Algorithm 4.

Input: ${\rm Dimension}:D;{\rm Search\ space}:S;{\rm Population\ size}:PS;{\rm
Generation}:T$

Output: ${\rm Best\ solution}:E$

Function  _ES(_$D,S,PS,T$_)_:

$t\leftarrow 0$

$\blacktriangleright$ (Population initialization)

$P_{t}\leftarrow\textbf{Initial}(D,S,PS)$

$FP_{t}\leftarrow\textbf{Evaluate}(P_{t})$

$E\leftarrow\textbf{bestIndividual}(P_{t},F_{t})$

while  _$t <T$ and not convergence_ do

$\blacktriangleright$ (Selection, Crossover and Mutation)

$O_{t}\leftarrow\textbf{Selection}(P_{t})$

$O_{t}\leftarrow\textbf{Crossover}(O_{t})$

$O_{t}\leftarrow\textbf{Mutation}(O_{t})$

$FO_{t}\leftarrow\textbf{Evaluate}(O_{t})$

$P_{t+1},FP_{t+1}\leftarrow\textbf{Offspring}(O_{t},FO_{t},P_{t},FP_{t})$

$E\leftarrow\textbf{bestIndividual}(P_{t+1},FP_{t+1})$

$t\leftarrow t+1$

end while

$\textbf{return}\ E$

Algorithm 4 ES

####  II-A5 PSO

PSO simulates the behavior of bird flocking, fish schooling, and swarming
theory to realize the optimization process[16]. The original version of PSO as
developed by the authors comprises a very simple concept, and paradigms can be
implemented in a few lines of computer code. It requires only primitive
mathematical operators and is computationally inexpensive in terms of both
memory requirements and speed. The Pseudocode of PSO is shown in Algorithm 5.

Input: ${\rm Dimension}:D;{\rm Search\ space}:S;{\rm Population\ size}:PS;{\rm
Generation}:T$

Output: ${\rm Best\ solution}:gbest$

Function  _PSO(_$D,S,PS,T$_)_:

$t\leftarrow 0$

$\blacktriangleright$ (Initialize the position and velocity)

$X_{t},V_{t}\leftarrow\textbf{Initial}(D,S,PS)$

for  _$i=0\ to\ PS$_ do

$pbest_{i}\leftarrow X_{t,i}$

end for

$F_{t}\leftarrow\textbf{Evaluate}(X_{t})$

$gbest\leftarrow\textbf{bestIndividual}(X_{t},F_{t})$

while  _$t <T$ and not convergence_ do

$\blacktriangleright$ (Optimization)

for  _$i=0\ to\ PS$_ do

$X_{t+1,i}\leftarrow\textbf{Update}(X_{t,i},V_{t,i})$

$V_{t+1,i}\leftarrow\textbf{Update}(V_{t,i},pbest_{i},gbest)$

end for

$F_{t+1}\leftarrow\textbf{Evaluate}(X_{t+1})$

for  _$i=0\ to\ PS$_ do

if  _$F_{t+1,i} >F_{pbest,i}$_ then

$pbest_{i}\leftarrow X_{t+1,i}$

end if

if  _$F_{pbest,i} >F_{gbest,i}$_ then

$gbest_{i}\leftarrow pbest_{i}$

end if

end for

$t\leftarrow t+1$

end while

$\textbf{return}\ gbest$

Algorithm 5 PSO

We update the position and velocity by Eq (1) in lines 13 and 14 of the
Algorithm 5

$\displaystyle x_{id}^{k+1}=x_{id}^{k}+v_{id}^{k+1}$ (1) $\displaystyle
v_{id}^{k+1}=wv_{id}^{k}+c_{1}r_{1}(p_{id,pbest}^{k}-x_{id}^{k})+$
$\displaystyle c_{2}r_{2}(p_{d,gbest}^{k}-x_{id}^{k})$

$i$ denotes the $i^{th}$ individual, $d$ denotes in the $d^{th}$ dimension,
$k$ denotes in the $k^{th}$ generation. In the velocity update equation, $w$
is a inertia weight, $c_{1}$ and $c_{2}$ are learning factor, $r_{1}$ and
$r_{2}$ are random number.

###  II-B Estimation of a convergence point

The concept of estimation of a convergence point was first proposed by Murata.
The paper[9] hypothesize that: in a population-based optimization algorithm,
all individuals are moving towards the global optimal solution. According to
this hypothesis, the nearest point to the extensions of all their movement
directions should locate near the global optimal point. Fig.1 shows how the
estimation of a convergence point works. But in practice, the estimated point
is not exactly on the global optimum due to incorrect directions or inaccurate
directions of population movements. However, it is highly expected that the
estimation point is close to the global optimum at least for unimodal
functions.

Figure 1: Moving vector $\bm{d}_{i}(=\bm{o}_{i}-\bm{p}_{i})$ is calculated
from a parent (worse) individual $\bm{p}_{i}$ and its offspring (better)
$\bm{o}_{i}$. The $\star$ is the estimated convergence point.

Let us derive how to estimate the convergence point mathematically. First, we
define parent (worse) individual $\bm{p}_{i}$, offspring (better)
$\bm{o}_{i}$, and moving vector $\bm{d}_{i}$ as describe in Fig 1. The unit
direction vector of $\bm{d}_{i}$ is given as
$\bm{d}_{0i}=\frac{\bm{d}_{i}}{||\bm{d}_{i}||}$,i.e., $\bm{d}_{0i}^{\rm
T}\bm{d}_{0i}=1$. $\mathbf{X}$ denotes the estimated convergence point, and
$\bm{p}_{i}+t_{i}\bm{d}_{i}$ represents the expansion from parent individual
$\bm{p}_{i}$ with the direction $\bm{d}_{i}$. $L(\mathbf{X},{t_{i}})$ in Eq
(2) becomes the minimum.

$\displaystyle\min(L(\mathbf{X},{t_{i}}))=\min(\sum_{i=1}^{n}||\bm{p}_{i}+t_{i}\bm{d}_{i}-\mathbf{X}||^{2})$
(2)

As the minimum line segment from the convergence point $\mathbf{X}$ to the
expansion line segments is the orthogonal projection from $\mathbf{X}$, we can
apply the Eq (3) into Eq (2) to remote $t_{i}$.

$\displaystyle\bm{d}_{i}^{T}(\bm{p}_{i}+t_{i}\bm{d}_{i}-\mathbf{X})=0\ ({\rm
orthogonal\ condition})$ (3)

Finally, the convergence point $\mathbf{X}$ can be calculated by Eq (4). See
detail expansion of equations in paper [9].

$\displaystyle\widehat{\mathbf{X}}=\left\\{\sum_{i=1}^{n}(\bm{I}_{d}-\bm{d}_{0i}\bm{d}_{0i}^{\rm
T})\right\\}^{-1}\left\\{\sum_{i=1}^{n}(\bm{I}_{d}-\bm{d}_{0i}\bm{d}_{0i}^{\rm
T})\bm{p}_{i}\right\\}$ (4)

##  III Proposal

In this Section, we will introduce our proposal in detail. First, we propose
two approximate methods for estimating the convergence points with Gaussian
sampling, which are shown in Fig 2

Figure 2: (a).elite sub-population averaging strategy for approximating the
convergence point. (b).elite sub-population weighted averaging strategy for
approximating the convergence point. The orange circle is the Gaussian
sampling space with the mean of estimated convergence point.

First, we proposed a hypothesis that there is a higher probability for a
optimum exists around an elite sub-population (individuals with higher
fitness). Based on this hypothesis, we select the elite sub-population from
the individual population of each generation from DE, and apply the averaging
strategy and the weighted averaging strategy to approximately estimate the
convergence point. Algorithm 6 and Algorithm 7 show the procedure of averaging
strategy and weighted averaging strategy.

Input: ${\rm Population}:P;{\rm Elite\ rate}:r;{\rm Dimension}:D$

Output: ${\rm Estimated\ convergence\ point}:C$

Function  _AS(_$P,r$_)_:

$EP\leftarrow\textbf{bestIndividuals}(P,r)$

$s\leftarrow\textbf{size}(EP)$

for  _$i=0\ to\ D$_ do

$C_{i}\leftarrow 0$

for  _$j=0\ to\ s$_ do

$C_{i}\leftarrow C_{i}+EP_{i,j}$

end for

$C_{i}\leftarrow C_{i}/s$

end for

$\textbf{return}\ C$

Algorithm 6 Averaging strategy

Input: ${\rm Population}:P;{\rm Elite\ rate}:r;{\rm Dimension}:D$

Output: ${\rm Estimated\ convergence\ point}:C$

Function  _WAS(_$P,r$_)_:

$EP\leftarrow\textbf{bestIndividuals}(P,r)$

$EF\leftarrow\textbf{Evaluate}(EP)$

$s\leftarrow\textbf{size}(EP)$

$F\leftarrow 0$

for  _$i=0\ to\ s$_ do

$F\leftarrow F+EF_{i}$

end for

$\blacktriangleright$ (Calculate the weights)

for  _$i=0\ to\ s$_ do

$w_{i}\leftarrow EF_{i}/F$

end for

for  _$i=0\ to\ D$_ do

$C_{i}\leftarrow 0$

for  _$j=0\ to\ s$_ do

$C_{i}\leftarrow C_{i}+w_{j}*EP_{i,j}$

end for

end for

$\textbf{return}\ C$

Algorithm 7 Weighted averaging strategy

In Algorithm 7, we calculate the weights based on normalized fitness of
individuals. After the convergence point is estimated, we apply a Gaussian
sampling operator with the mean of the convergence point and a certain
standard deviation $\sigma$ to find $k$ individuals. Then we select the best
$k$ individuals among 1 estimated convergence point and $k$ individuals found
by Gaussian sampling operator and worst $k$ individuals generated by DE as a
part of offspring to participate the optimization.

##  IV Numerical Experiment

Here, we define shortened names of our proposal:

(a) P1: proposal with averaging strategy and Gaussian sampling operator.

(b) P2: proposal with weighted averaging strategy and Gaussian sampling
operator.

We apply 28 benchmark functions from the CEC2013 Suite[11] with independent 30
trial runs in this evaluation experiment. Table I shows the details of these
functions.

TABLE I: CEC2013 Suite: Uni=unimodal, Multi=multimodal, Comp=Composition Func. | Types | Characteristics | Optimum
---|---|---|---
$f_{1}$ | Uni | Sphere function | -1400
$f_{2}$ | Rotated high conditioned elliptic function | -1300
$f_{3}$ | Rotated Bent Cigar function | -1200
$f_{4}$ | Rotated discus function | -1100
$f_{5}$ | Different powers function | -1000
$f_{6}$ | Multi | Rotated Rosenbrock fs function | -900
$f_{7}$ | Rotated Schaffers function | -800
$f_{8}$ | Rotated Ackley fs function | -700
$f_{9}$ | Rotated Weierstrass function | -600
$f_{10}$ | Rotated Griewank fs function | -500
$f_{11}$ | Rastrigin fs function | -400
$f_{12}$ | Rotated Rastrigin fs function | -300
$f_{13}$ | Non-continuous rotated Rastrigin fs function | -200
$f_{14}$ | Schwefel fs function | -100
$f_{15}$ | Rotated Schwefel fs function | 100
$f_{16}$ | Rotated Katsuura function | 200
$f_{17}$ | Lunacek bi-Rastrigin function | 300
$f_{18}$ | Rotated Lunacek bi-Rastrigin function | 400
$f_{19}$ | Expanded Griewank fs plus Rosenbrock fs function | 500
$f_{20}$ |  Expanded Schaffer fs $f_{6}$ function | 600
$f_{21}$ | Comp |  Composition function 1 ($n$ = 5, rotated) | 700
$f_{22}$ |  Composition function 2 ($n$ = 3, unrotated) | 800
$f_{23}$ |  Composition function 3 ($n$ = 3, rotated) | 900
$f_{24}$ |  Composition function 4 ($n$ = 3, rotated) | 1000
$f_{25}$ |  Composition function 5 ($n$ = 3, rotated) | 1100
$f_{26}$ |  Composition function 6 ($n$ = 5, rotated) | 1200
$f_{27}$ |  Composition function 7 ($n$ = 5, rotated) | 1300
$f_{28}$ |  Composition function 8 ($n$ = 5, rotated) | 1400

And we compare our proposal with RS, GA, DE, ES and PSO. The parameter of
these algorithm is shown in Table II

TABLE II: The parameter of comparison methods Alg. | Parameter | Value
---|---|---
Common parameters | dimension | 2-D, 10-D, 30-D
population size | 50 * dimension
max evaluation times | 1000 * dimension
search space | [-100, 100] * dimension
RS | sampling strategy | random sampling
GA | strategy | Elitist GA
crossover rate | 0.5
mutation rate | 0.1
DE | strategy | DE/current-to-best/1
scale factor | 0.7
crossover rate | 0.9
ES | strategy | (1+1)-ES
mutation strength | 0.2
PSO | $w$ | 0.9
$c_{1},c_{2}$ | 2
Proposal | strategy | DE/current-to-best/1
scale factor | 0.7
crossover rate | 0.9
elite proportion | 0.05
$\sigma$ in Gaussian sampling | 5

At the end of optimization, we apply the Kruskal-Wallis test and the Holm
multiple comparison test on the fitness values. Table III shows their results
of the statistical tests. In Holm multiple comparison test, if our proposal is
significant better than the second-best algorithm or original DE, we apply
$\gg$ (p$<$0.01) and $>$ (p$<$0.05) to denote the significance, $\approx$
represents there is no significance between two comparing methods.

TABLE III: CEC2013 Suite: Uni=unimodal, Multi=multimodal, Comp=Composition Func. | 2-D | 10-D | 30-D
---|---|---|---
P1 | P2 | P1 with DE | P2 with DE | P1 with P2 | P1 | P2 | P1 with DE | P2 with DE | P1 with P2 | P1 | P2 | P1 with DE | P2 with DE | P1 with P2
$f_{1}$ |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $>$ P2 |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{2}$ |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\gg$ DE |  P2 $\gg$ DE |  P2 $\gg$ P1 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{3}$ |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{4}$ |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $>$ DE |  P2 $>$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{5}$ |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $>$ DE |  P2 $>$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{6}$ |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $>$ DE |  P2 $>$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{7}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{8}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2
$f_{9}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $>$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2
$f_{10}$ |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\approx$ P1 |  PSO $\approx$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{11}$ |  P1 $\approx$ PSO |  P2 $>$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ GA |  P2 $\approx$ GA |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{12}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{13}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{14}$ |  P1 $\approx$ GA |  P2 $\approx$ GA |  P1 $\approx$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  GA $\gg$ P1 |  GA $\gg$ P2 |  P1 $\approx$ DE |  P2 $>$ DE |  P1 $\approx$ P2
$f_{15}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $>$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{16}$ |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\approx$ GA |  P2 $\approx$ GA |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\approx$ ES |  P2 $\approx$ ES |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2
$f_{17}$ |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2
$f_{18}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $>$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{19}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\gg$ PSO |  P2 $>$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{20}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $>$ DE |  P2 $>$ DE |  P1 $\approx$ P2
$f_{21}$ |  P1 $\gg$ PSO |  P2 $\gg$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $>$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{22}$ |  P1 $\approx$ GA |  P2 $\approx$ GA |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $>$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  GA $\gg$ P1 |  GA $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{23}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2
$f_{24}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ ES |  P2 $\approx$ ES |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{25}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $>$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2
$f_{26}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  P1 $\approx$ GA |  P2 $>$ GA |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\approx$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{27}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\gg$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\approx$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2
$f_{28}$ |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $\approx$ DE |  P2 $\approx$ DE |  P1 $\approx$ P2 |  P1 $\approx$ PSO |  P2 $\approx$ PSO |  P1 $>$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2 |  PSO $\gg$ P1 |  PSO $\gg$ P2 |  P1 $\approx$ DE |  P2 $\gg$ DE |  P1 $\approx$ P2

From the experiment results, we can see that our proposal can significantly
accelerate the DE in most test functions on the CEC2013 Suite, and P1 has no
siginificant difference with P2 in most functions. In the 2-D space, our
proposal is competitive with the comparing 5 methods, while in the 10-D and
50-D space, although PSO outperforms our proposal combined with DE on most
test functions, this is mainly due to the limited performance of DE on these
functions, our proposal can still accelerate the DE in optimization. Different
from the method proposed by Yu et al[10], it is unnecessary for our proposal
to calculate the moving vector, so it has stronger scalability for our
proposal combined with other evolutionary algorithms.

##  V Conclusion

In this paper, we propose a simple strategy to estimate the convergence point
approximately. Based on the averaging strategy, ordinary averaging strategy
and weighted averaging strategy are derived. We also design a Gaussian
sampling operator based on the estimated convergence point and combine it with
traditional DE. Numerical experiments show our proposal can improve the
performance of traditional DE. Although our proposal combined with DE is worse
than PSO in some functions with 10-D and 50-D, the main reason is the
performance of DE is limited in these functions, and our proposal has strong
scalability that can be extended to all population-based heuristic algorithms
in theory.

In future research, a well-performed local search operator replacing the
Gaussian sampling is a promising topic, such as the introduction of CMA-ES or
the adaptive evolution strategy. And in this paper, we only apply our proposal
in low dimensional test functions. In future research, we will combine our
proposal with cooperative co-evolution to solve large-scale optimization
problems. Finally, the estimation of a convergence point combined with the
Gaussian sampling operator has a broad prospect to solve optimization problems
and can be easily extended with other population-based evolutionary
algorithms.

##  VI Acknowledgement

This work was supported by JSPS KAKENHI Grant Number JP20K11967.

## References

  * [1] M. A. Ahandani, N. P. Shirjoposh, and R. Banimahd, “Three modified versions of differential evolution algorithm for continuous optimization,” _Soft Computing_ , vol. 15, no. 4, pp. 803–830, 2010.
  * [2] Z. Yang, K. Tang, and X. Yao, “Differential evolution for high-dimensional function optimization,” 2007, pp. 3523–3530.
  * [3] F. Wang, W. Wang, T. Feng, and H. Li, “Parallel test sheets generation using differential evolution algorithm with constraint effective encoding and either-or mutation,” in _2015 2nd International Conference on Information Science and Control Engineering_. IEEE, 2015, pp. 376–379.
  * [4] K. Deb, “Multi-objective optimisation using evolutionary algorithms: an introduction,” in _Multi-objective evolutionary optimisation for product design and manufacturing_. Springer, 2011, pp. 3–34.
  * [5] D. Wolpert and W. Macready, “No free lunch theorems for optimization,” _IEEE Transactions on Evolutionary Computation_ , vol. 1, no. 1, pp. 67–82, 1997.
  * [6] M. Norman, P. Moscato, and L. Plata, “A competitive-cooperative approach to complex combinatorial search,” 1991.
  * [7] A. Abdolmaleki, B. Price, N. Lau, L. P. Reis, and G. Neumann, “Deriving and improving cma-es with information geometric trust regions,” in _Proceedings of the Genetic and Evolutionary Computation Conference_ , 2017, pp. 657–664.
  * [8] F. Fröhlich and P. K. Sorger, “Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models,” _PLOS Computational Biology_ , vol. 18, no. 7, p. e1010322, 2022.
  * [9] N. Murata, R. Nishii, H. Takagi, and Y. Pei, “Analytical estimation of the convergence point of populations,” in _2015 IEEE Congress on Evolutionary Computation, CEC 2015 - Proceedings_. Institute of Electrical and Electronics Engineers Inc., Sep. 2015, pp. 2619–2624.
  * [10] J. Yu, Y. Li, Y. Pei, and H. Takagi, “Accelerating evolutionary computation using a convergence point estimated by weighted moving vectors,” _Complex and Intelligent Systems_ , vol. 6, 05 2019.
  * [11] J. Liang, B. Qu, P. Suganthan, and A. HernÃ¡ndez-DÃaz, “Problem definitions and evaluation criteria for the cec 2013 special session on real-parameter optimization,” _Technical Report 201212, Computational Intelligence Laboratory, Zhengzhou University, Zhengzhou China_ , 01 2013.
  * [12] L. A. Rastrigin, “Convergence of random search method in extremal control of multi-parameter systems,” _Avtomat.i Telemekh_ , p. 1467â1473, 1963.
  * [13] J. H. Holland, _Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence_. MIT press, 1992.
  * [14] R. Storn, “On the usage of differential evolution for function optimization,” in _Proceedings of north american fuzzy information processing_. Ieee, 1996, pp. 519–523.
  * [15] L. B. Slobodkin, “The strategy of evolution,” _American scientist_ , vol. 52, no. 3, pp. 342–357, 1964.
  * [16] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in _Proceedings of ICNN’95 - International Conference on Neural Networks_ , vol. 4, 1995, pp. 1942–1948 vol.4.