# CCO-VOXEL: Chance Constrained Optimization over Uncertain Voxel-Grid
Representation for Safe Trajectory Planning

Sudarshan S Harithas1, Rishabh Dev Yadav1, Deepak Singh1, Arun Kumar Singh2, K
Madhava Krishna1 1are with RRC, IIIT Hyderabad, India<EMAIL_ADDRESS>{sudarshan.s<EMAIL_ADDRESS><EMAIL_ADDRESS>with University of Tartu, Estonia
<EMAIL_ADDRESS>authors thank the anonymous reviewers for their helpful
comments and Ministry of Electronics and Information Technology (MeitY),
Government of India, for their generous financial support. This work was
supported in part by the European Social Fund via IT Academy program in
Estonia, and grant PSG753 from Estonian Research Council. †Project page:
https://github.com/sudarshan-s-harithas/CCO-VOXEL

###### Abstract

We present CCO-VOXEL: the very first chance-constrained optimization (CCO)
algorithm that can compute trajectory plans with probabilistic safety
guarantees in real-time directly on the voxel-grid representation of the
world. CCO-VOXEL maps the distribution over the distance to the closest
obstacle to a distribution over collision-constraint violation and computes an
optimal trajectory that minimizes the violation probability. Importantly,
unlike existing works, we never assume the nature of the sensor uncertainty or
the probability distribution of the resulting collision-constraint violations.
We leverage the notion of Hilbert Space embedding of distributions and Maximum
Mean Discrepancy (MMD) to compute a tractable surrogate for the original
chance-constrained optimization problem and employ a combination of A* based
graph-search and Cross-Entropy Method for obtaining its minimum. We show
tangible performance gain in terms of collision avoidance and trajectory
smoothness as a consequence of our probabilistic formulation vis a vis state-
of-the-art planning methods that do not account for such non-parametric noise.
Finally, we also show how a combination of low-dimensional feature embedding
and pre-caching of Kernel Matrices of MMD allow us to achieve real-time
performance in simulations as well as in implementations on on-board commodity
hardware that controls the quadrotor flight.

©2022 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works

##  I INTRODUCTION

Many aerial navigation systems compute trajectories over Octomap [1], or
Euclidean Signed-distance field (ESDF) [2] derived from a three-dimensional
voxel grid representation of the world. These approaches assume accurate
sensor readings during the process of map creation. However, sensor
measurements (e.g., point clouds) are noisy and often non-parametric in
nature. These noisy measurements get translated to incorrect estimates of the
distance to the closest obstacle and gradient field. Consequently,
deterministic planners [3, 4, 5, 6] who rely on these two pieces of
information being accurate are expected to lose their performance and safety
guarantees under perception noise. For example, noisy gradients coupled with
an incorrect estimate of the distance to the nearest obstacle can result in
the UAV moving towards the obstacle and colliding. This paper proposes a
chance-constrained optimization (CCO) over VOXEL grids as an antidote to this
duress. CCO-VOXEL acknowledges that the distance to the closest occupied voxel
is noisy and consequently takes this uncertainty into account while computing
optimal trajectory plans. As a result, CCO-VOXEL avoids collisions where
gradient-based deterministic planner fails (cf. LABEL:teaser).

Previous methods leveraging CCO for handling sensor noise had to rely on the
assumption that the underlying random variables that constitute the chance
constraint function are parametric or Gaussian [7, 8, 9, 10]. However, as seen
in Fig. (2) the distribution over the nearest voxel occupancy is non-
parametric. Moreover, existing works on CCO such as [10] assume that the
obstacles have a pre-defined shape which is crucial for obtaining tractable
reformulations of the chance constraints. In this paper, we relax both the
assumptions of parametric uncertainty and the analytical collision model. Yet,
we maintain real-time performance in both simulation and hardware experiments
with on-board computation. The core innovations behind our approach are
summarized below.

Algorithmic Contribution: For the first time, we present a CCO-based
probabilistically safe trajectory planning that is agnostic to both the nature
of the underlying uncertainty and the obstacle geometry. Inspired by our prior
works [11, 12], we use the notion of distribution embedding in Reproducing
Kernel Hilbert Space (RKHS) and Maximum Mean Discrepancy (MMD) to reformulate
CCO as a distribution matching problem. However, we overcome one of the key
limitations of [11, 12]: estimating the so-called desired (or target)
distribution for computing the distribution matching cost. We show that by
embedding the distribution of constraint violation in RKHS instead of
collision avoidance constraint themselves, as done in [11, 12], we can by-pass
the need to estimate the desired distribution. Furthermore, distribution of
constraint violation induces efficient structure in the algebraic form for
MMD, allowing for pre-computation and caching of most computationally
expensive parts. We also improve the computational performance of our approach
by projecting the samples drawn from the distribution of constraint violations
into some latent dimension before computing the MMD. We learn this projection
through auto-encoder-based supervised learning. Finally, we minimize a
combination of MMD and kinodynamic costs through the Cross-Entropy Method
(CEM) to compute a smooth and probabilistically safe trajectory.

State-of-the-art Performance: We show tangible performance gain over
deterministic gradient planners [3] and methods that handle uncertainty by
growing the obstacle map by the co-variance of the distribution [13, 14]. The
benefits are seen in the form of almost negligible collisions under the duress
of varying levels of sensor noise even as the previous methods reports high
collision rates.

##  II Related Work

Trajectory planning for quadrotor systems has a long history in the areas of
optimal trajectory generation [15, 7, 16] and integrated perception and
navigation [17, 18]. Gradient based planning over voxel representations have
become popular in recent years [3, 6] though they have not handled uncertainty
in their formulations. Meanwhile, Chance Constrained Optimization with
suitable surrogates too have been well studied, however their inadequacy in
handling non-parametric uncertainty was detailed in [11, 12]. Specifically
there appears no known methods that handle non parametric uncertainty of
discretized volumetric voxel grids other than the proposed. While [19]
provides for a chance constrained framework for nearest distance to obstacle
it formulates over well defined obstacle shapes and not over voxel grids.
Moreover, its computation of safe trajectories runs into minutes while the
proposed method is in the order of tens of milliseconds.

Figure 2: The error (difference between the true collision distance and the
measured value) distribution histogram that is obtained from a noisy voxel
grid. The OptiTrack Motion Capture system was used to determine the ground
truth distance. Octomap with Euclidean Distance Transform (dynamicEDT3D) is
used for collision distance measurements.

##  III Problem Formulation

Assumptions: We assume that the robot has precise motion capabilities but
noisy perception. As a result, the built ESDF has imprecise information about
the gradient field and distance to the closest obstacle. The distribution over
the distance has the form $d=d_{measured}+\epsilon_{i}$, where $\epsilon_{i}$
are samples drawn from some black-box distribution whose parametric form is
not known.

Trajectory Parametrization: We parametrize quadrotor trajectory in the
following manner

$\small\begin{bmatrix}x(t_{1})\\\ x(t_{2})\\\ \dots\\\
x(t_{n})\end{bmatrix}=\textbf{P}\textbf{c}_{x},\begin{bmatrix}\dot{x}(t_{1})\\\
\dot{x}(t_{2})\\\ \dots\\\
\dot{x}(t_{n})\end{bmatrix}=\dot{\textbf{P}}\textbf{c}_{x},\begin{bmatrix}\ddot{x}(t_{1})\\\
\ddot{x}(t_{2})\\\ \dots\\\
\ddot{x}(t_{n})\end{bmatrix}=\ddot{\textbf{P}}\textbf{c}_{x}.\vspace{-1mm}$
(1)

where, $\textbf{P},\dot{\textbf{P}},\ddot{\textbf{P}}$ are matrices formed
with time-dependent basis functions (e.g polynomials) and $\textbf{c}_{x}$ are
the coefficients associated with the basis functions. Similar expressions can
be written for $y(t),z(t)$ as well in terms of coefficients
$\textbf{c}_{y},\textbf{c}_{z}$, respectively.

Safe Motion Planning: Using the parametrization discussed above, we formulate
our probabilistically safe motion planning as the following optimization
problem, wherein $P(\cdot)$ represents probability.

$\displaystyle\min l(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})$ (2a)
$\displaystyle P(f(d)\leq 0)\geq\eta,\forall t$ (2b) $\displaystyle
f(d)=r_{safe}-d(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z},t)\vspace{-4mm}$
(3)

Our cost function (2a) encourages smoothness in the trajectory by mapping the
trajectory coefficients to penalties on higher-order motion derivatives such
as jerk and penalising the violation of acceleration and velocity limit. The
algebraic form is similar to that used in [3]. The function $d(\cdot)$ in (3)
maps the coefficients to distance to the closest obstacle at a given time $t$.
Consequently, the inequality (2b) ensures that the probability of the closest
obstacle distance being less than safe value $r_{safe}$ is below some
specified threshold $\eta$. Constraints of the form (2b) are known as chance
constraints [20]. Generally, chance-constrained optimization (CCO) are
considered computationally intractable. Thus, existing works focuses on
deriving tractable approximations of (2b) under two predominant assumption.
Firstly, the obstacles are assumed to have defined geometric shapes (e.g
ellipsoid) and as a result it becomes possible to derive an analytical
expression for $d(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z},t)$. Second,
the perception uncertainty are often assumed to be Gaussian. When both the
assumptions are made simultaneously it even becomes possible to formulate a
convex approximation of (2b) [10], [21].

Our work relaxes both the above mentioned assumptions and thus substantially
expands the applicability of chance-constrained optimization based motion
planning in real-world environment. We handle chance-constraints on collision
avoidance by treating $d(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z},t)$ as a
black-box model that can only be queried from a voxel-grid representation of
the world. Secondly, we make no assumption on the underlying perception
uncertainty.

Reformulation from [11], [12] and Limitations: Let $p_{f}$ represent the
probability distribution of $f(d)$, i.e $p_{f}=P(f(d))$. Our prior work [11],
[12] imagined CCO as the problem of choosing an appropriate
$(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})$ such that the $p_{f}$ at any
given time $t$ takes on an appropriate shape. This naturally gives rise to the
notion of desired distribution $p_{f_{des}}$, i.e a distribution that $p_{f}$
should resemble as closely as possible (cf. Fig. 2(a)). The entire mass of the
$p_{f_{des}}$ lies to the left of $f(d)=0$. As $p_{f}$ becomes similar to
$p_{f_{des}}$, its mass starts getting shifted to the left and the probability
of $P(f(d)\leq 0)$ goes up. This intuition led to the following reformulation
of optimization (2a)-(2b)

$\displaystyle\min
l(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})+wl_{dist}(p_{f},p_{f_{des}}),\vspace{-3.0mm}$
(4)

where, $l_{dist}$ is a cost that quantifies the similarity between two
distributions. The constant $w$ controls the trade-off between minimizing the
distribution similarity (or probability of collision avoidance) with the
primary cost. We discuss suitable form for $l_{dist}$ in the next section.

The main limitation of [11], [12] is that it requires one to solve an
optimization problem to estimate $p_{f_{des}}$. While tractable, in a reactive
one-step setting, this estimation process becomes highly challenging for
multi-step trajectory planning such as the one considered in this paper. In
the next section, we introduce our main algorithmic results: a novel RKHS
reformulation of the CCO (3)-(2b) that does not require us to estimate
$p_{f_{des}}$.

(a) The setting from [11] where the goal is to bring the distribution $p_{f}$
to the left of $f(d)=0$. This is achieved by constructing a desired
distribution $p_{f_{des}}$ and minimizing the distance of $p_{f}$ from it in
RKHS. Importantly, estimation of $p_{f_{des}}$ required solving an
optimization problem

(b) We work with the distribution of the collision constraint violation rather
than the constraint themselves. As a result, the desired distribution in our
case is simply the Dirac Delta distribution $p_{\delta}$. The distribution
$p_{\overline{f}}$ approaches the desired Dirac function $p_{\delta}$ for an
optimal $\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z}$

Figure 3: An intuitive understanding of the RKHS embedding of [11] and our
proposed method

##  IV Main Algorithmic Results

###  IV-A Distribution over Constraint Violations

Our main idea is to work with the distribution of violations of the collision
avoidance constraints rather than the constraints themselves. More precisely,
let $\overline{f}$ be the constraint violation function defined as

$\displaystyle\overline{f}=\max(0,f(d(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z},t)))\vspace{-3mm}$
(5)

A distribution over $d$ can be mapped to $\overline{f}$ and we represent it as
$p_{\overline{f}}$. Although, computing the analytical form for
$p_{\overline{f}}$ is intractable, we can make the following remark regarding
the best possible shape that it can take in the context of collision
avoidance.

######  Remark 1.

The best possible shape of $p_{\overline{f}}$ is given by a Dirac Delta
distribution $p_{\delta}$ . In other words, $p_{\delta}$ acts as the desired
distribution for $p_{\overline{f}}$

######  Remark 2.

As $p_{\overline{f}}$ becomes more and more similar to $p_{\delta}$, the
probability of collision avoidance goes up.

Remark 1 and 2 elucidates our motivation behind shifting from the distribution
of constraints, as done in [11], [12], to distribution of constraint
violations. The desired distribution in our case is always fixed and most
importantly, exactly known. Remark 2 forms the basis of the following proposed
reformulation of CCO.

$\min_{\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z}}l(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})+wl_{dist}(p_{\overline{f}},p_{\delta})\vspace{-2mm}$
(6)

Note that $p_{\overline{f}}$ is in fact a function of $d$ which in turn is a
function of the trajectory coefficients
$(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})$

###  IV-B Maximum Mean Discrepancy and Sample Approximation

A common method to quantify similarity between two distributions is through
Kullback Liebler Divergence (KLD). However, it requires us to know the
analytical form for the distributions under consideration and thus is not
suitable to quantify the similarity between $p_{\overline{f}}$ and
$p_{\delta}$. One possible workaround is provided by the notion of Maximum
Mean Discrepancy (MMD) that quantifies the similarity of two distributions in
Reproducing Kernel Hilbert Space (RKHS). Importantly, both the embedding into
RKHS and MMD can be obtained by just sample level information of
distributions. To this end, let $\mu_{p_{\overline{f}}}$ and
$\mu_{p_{\delta}}$ represent the RKHS embedding of $p_{\overline{f}}$ and
$p_{\delta}$ computed in the following manner.
$\small\mu_{p_{\overline{f}}}=\sum_{i=0}^{N}\alpha_{i}k(\overline{f}(d_{i}),\cdot)$
(7) $\mu_{p_{\delta}}=\sum_{i=0}^{N}\beta_{i}k(0,\cdot)$ (8)

where, $d_{i}$ refers to the $i^{th}$ sample drawn from the distribution of
$d$ for a given $(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z},t)$. Please
note that (8) follows from the fact the samples from a Dirac Delta
distribution are all zeroes. Using (7) and (8), we can formulate $l_{dist}$ in
the following manner.

$\small
l_{dist}=\left\|\mu_{p_{\overline{f}}}-\mu_{p_{\delta}}\right\|_{2}^{2}\vspace{-1mm}$
(9)

It is worth reiterating that MMD (9) is a function of
$(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z},t)$. Each possible choice leads
to a different distribution of $d$ and consequently MMD. Thus the goal is to
come up with the right trajectory coefficients which minimizes both MMD and
the primary cost function. We achieved this through Cross-Entropy Method [22]
[23] and is discussed in Sec. IV-D.

###  IV-C Pre-Computation and Improving Sample Complexity

Computing MMD through Kernel Matrices: Let ${d_{0},d_{1},d_{2},\dots,d_{n}}$
be the samples drawn from the distribution of the closest obstacle at any
given time $t$ for a given $\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z}$. The
distance samples can then be mapped to samples from $p_{\overline{f}}$ through
(5). In this part, we derive the matrix representation of the MMD (8) to show
how various parts of it can be pre-computed and cached. To this end, we expand
(9) in the following manner.

$\displaystyle\small\lVert\mu_{p_{\overline{f}}}-\mu_{p_{\delta}}\lVert^{2}=$
$\displaystyle\langle{\mu_{p_{\overline{f}}}(d)},{\mu_{p_{\overline{f}}}(d)}\rangle-2\langle{\mu_{p_{\overline{f}}}},{\mu_{p_{\delta}}}\rangle+$
$\displaystyle\langle{\mu_{p_{\delta}}},{\mu_{p_{\delta}}}\rangle\vspace{-2mm}$
(10)

Substituting kernel mean functions (7) and (8) in (10)

$\displaystyle\langle\mu_{p_{\overline{f}}}(d),\mu_{p_{\overline{f}}}\rangle=\langle\sum_{i=0}^{N}{\alpha_{i}k(\overline{f}(d_{i}),\cdot},\sum_{j=0}^{N}{\alpha_{j}k(\overline{f}(d_{j}),\cdot}\rangle$
(11a)
$\displaystyle\langle\mu_{p_{\overline{f}}}(d),\mu_{p_{\delta}}\rangle=\langle\sum_{i=0}^{N}{\alpha_{i}k(\overline{f}(d_{i}),\cdot},\sum_{j=0}^{N}\beta_{j}{k(0,\cdot)}\rangle$
(11b)
$\displaystyle\langle\mu_{p_{\delta}},\mu_{p_{\delta}}\rangle=\langle\sum_{i=0}^{N}{\beta_{i}k(0,\cdot)},\sum_{j=0}^{N}{\beta_{j}k(0,\cdot)}\rangle$
(11c)

Using the kernel trick Eq. (11) can be reduced to a matrix form and can be
expressed as

$\small\lVert\mu_{p_{\overline{f}}}-\mu_{p_{\delta}}\lVert^{2}=\mathbf{C_{\alpha}K_{\overline{f}\overline{f}}C_{\alpha}^{T}-2C_{\alpha}K_{\overline{f}0}C_{\beta}^{T}+C_{\beta}K_{00}C_{\beta}^{T}}\vspace{-1mm}$
(12)

where $C_{\alpha}$ and $C_{\beta}$ are the weight vectors given by

$\displaystyle
C_{\alpha}=\begin{bmatrix}\alpha_{0},\alpha_{1},...,\alpha_{n}\end{bmatrix};C_{\beta}=\begin{bmatrix}\beta_{0},\beta_{1},...,\beta_{n}\end{bmatrix}\vspace{-2mm}$

and $K_{\overline{f}\overline{f}}$, $K_{\overline{f}\delta}$ and
$K_{\delta\delta}$ are matrices defined as

$\footnotesize\mathbf{K_{\overline{f}\overline{f}}}=\begin{bmatrix}k(\overline{f}(d_{0}),\overline{f}(d_{0}))&\ldots&k(\overline{f}(d_{0}),\overline{f}(d_{n}))\\\
k(\overline{f}(d_{1}),\overline{f}(d_{0}))&\ldots&k(\overline{f}(d_{1}),\overline{f}(d_{n}))\\\
\vdots&\ddots&\vdots&\\\
k(\overline{f}(d_{n}),\overline{f}(d_{0}))&\ldots&k(\overline{f}(d_{n}),\overline{f}(d_{n}))\\\
\end{bmatrix}$
$\footnotesize\mathbf{K_{f\delta}}=\begin{bmatrix}k(\overline{f}(d_{0}),0)&\ldots&k(\overline{f}(d_{0}),0)\\\
k(\overline{f}(d_{1}),0)&\ldots&k(\overline{f}(d_{1}),0)\\\
\vdots&\ddots&\vdots&\\\
k(\overline{f}(d_{n}),0)&\ldots&k(\overline{f}(d_{n}),0)\vspace{-2mm}\end{bmatrix}$

$\mathbf{K_{\delta\delta}}=\mathbf{1}_{n\times n}\vspace{-2mm}$

We use Radial Basis Functions (RBF) during implementation. The use of the
Dirac-Delta function as the desired function leverages us with two key
properties which would further decrease the computation time.

  1. 1.

Pre computing matrices: The matrix $K_{\delta\delta}$ is essentially a set of
ones and it does not need explicit computation at run time. Furthermore, if a
polynomial kernel is used the matrix $K_{\bar{f}\delta}$ would also reduce to
set of ones and only one matrix needs computation at run time.

  2. 2.

Symmetric matrix: It can be observed the kernel matrices are square and
symmetric, this implies that we can compute values only for one triangle of
the matrix and copy the values into the other triangle.

The use of these techniques in addition to feature space embedding, and
vectorized matrix-algebra enables the computation of MMD in real-time.

Feature space embedding of distance measurements: We need a total of $n$
samples of $\overline{f}$ to compute its embedding in RKHS and consequently
MMD. A large $n$ $(\approx 100)$ gives a good estimate of MMD but at the same
time increases the computational burden (the size of Kernel matrix in (12)
increases). In this part, we use an auto-encoder to improve the sample and
computation complexity. More concretely, we stack the $n$ samples of
$\overline{f}$ at a given point on the trajectory and project it into a low
dimensional latent space. Our auto-encoder based embedding can be seen as the
more sophisticated version of reduced-set method used in [11].

We learn the low-dimensional mapping through supervised learning. Formally, we
seek to determine an encoder $W_{1}$ and decoder $W_{2}$ pair such that they
perform an identity operation. The training is performed by optimizing the
loss function given in (13) using Stochastic Gradient Descent.The training
data consists of $\mathnormal{n}$ sets of $\overline{f}$, each with
$\mathnormal{m}$ samples, stacked to form the data matrix D of dimension
$\boldsymbol{R^{n\times m}}$.

$\displaystyle\min_{\textbf{W}_{1},\textbf{W}_{2}}\left\|\textbf{D}\textbf{W}_{1}\textbf{W}_{2}-\textbf{D}\right\|^{2}\vspace{-5mm}$
(13)

The final latent space dimension where MMD would be computed is given by
$\textbf{D}\textbf{W}_{1}$

###  IV-D Trajectory Optimization

We minimize (6) in two steps. At the first step, we use graph-search technique
to find a trajectory that approximately minimizes (6). This resulting
trajectory is used to initialize the Cross Entropy Method (CEM) optimizer for
computing an optimal solution of (6).

Initial Trajectory Search: Our initial trajectory search builds on Kinodynamic
A* proposed in [3] but with a crucial difference that we use MMD as a part of
the Edge Cost that connects two nodes of the graph. The typical process is as
follows. Several motion primitives are generated by sampling different control
inputs $\textbf{u}_{t}$ (acceleration). These primitives are then assigned an
edge-cost (14), wherein $\Delta M$ (15) is the difference between MMD at the
current, and the next node and $\tau$ is the time for which the control input
is applied. The constant $\rho$ trades-off control costs with time. Apart from
the edge cost, there is also a heuristic cost that measures the goal-reaching
aspect of each primitive. The primitive with the best edge and heuristic cost
combination is chosen for expansion in the next iteration.

In contrast to our approach, [3] uses the distance to the closest obstacle (or
simply occupancy in an inflated map) to filter out unsafe primitives while the
edge-cost only has the first two terms from (14). As mentioned earlier, if
this distance information is noisy, the A* process itself can lead to unsafe
trajectories (see Table III).

$\small e_{c}=(\lVert u_{t}\lVert^{2}+\rho)\tau+\Delta M$ (14)

$\Delta M=M_{i+1}-M_{i}$ (15)

CEM Based Refinement: The trajectory produced by the graph search can be sub-
optimal and might have a higher arc length. We use Cross Entropy Minimization
(CEM) to further refine the trajectory and improve its smoothness while
maintaining clearance from the obstacle. The various steps of CEM are
summarized in Algorithm 1. The input to the CEM is the sub-optimal trajectory
obtained from graph-search. The initial trajectory is converted into an
initial estimate for trajectory coefficients
$(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})$ through simple curve-fitting
(line 1). The first step in Algorithm 1 is drawing $N$ samples of trajectory
coefficients, leading to $N$ possible trajectories (line 1). At each point of
every sampled trajectory, we computed the cost function (6) (line 1). Note
that as shown section IV-C, the MMD evaluation itself requires drawing $n$
samples from the collision distance distribution at each point of all the
sampled $N$ trajectories. We select the top $q$ trajectory coefficient samples
(line 1 ) that led to lowest cost and use them to compute the mean and
variance for the next iteration of CEM (line 1). The output of the CEM is the
coefficients of the optimal trajectory.

In a standard CEM implementation after every iteration of the inner loop (line
1 to 1) all the elite samples are discarded. This decreases the efficiency and
increases the convergence time. Authors in [23] propose CEM with memory where
a small subset of the elite set is stored and added to the pool of samples for
the next iteration. Our CEM uses this insight in line 1 of Algorithm 1.

1 Input: Sub optimal trajectory ($x(t),y(t),z(t)$) from initial path search
and trajectory execution time $t_{exe}$

2 $N_{CEM}$ = Number of CEM interations

3 $\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z}$ = Trajectory coefficients

4
$\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z}\leftarrow\textbf{FitPolynomial}($x(t),
y(t), z(t)$)$

5 $\boldsymbol{\mu}_{0}=(\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z})$

6 for  _$i=0,i\leq N_{CEM},i++$_ do

7 Initialize $CostList$ = []

8 $Pool=\text{N samples from }\mathcal{N}(\mu_{i},diag(\sigma_{i}^{2}))$

9 for  _$j=0,j\leq\text{size}(Pool),j++$_ do

10 $cost\leftarrow\textbf{DetermineCost}(S_{j})$

11 append ${cost}$ to $CostList$

12 end for

13 $EliteSet\leftarrow\textbf{SelectElites}(CostList,Pool)$

14 $Pool\leftarrow\textbf{UpdateSamples}(EliteSet)$

15 $(\mu_{i+1},\sigma_{i})\leftarrow\text{fit Gaussian distribution }EliteSet$

16 $\textbf{c}_{x},\textbf{c}_{y},\textbf{c}_{z}=\boldsymbol{\mu}_{i}$

17 end for

18$\text{Optimized
Trajectory}\leftarrow\mathbf{P}\mathbf{c_{x}},\textbf{P}\textbf{c}_{y},\textbf{P}\textbf{c}_{z}$

Algorithm 1 CEM based Trajectory Refinement

##  V Experiments and Validation

In this section, we demonstrate the advantages of CCO-VOXEL over deterministic
methods and other conservative approaches to handle uncertainty through
qualitative and quantitative comparisons. We further show that with the
convergence of CEM to an optimal solution the collision constraint violation
distribution will also converge to the dirac delta function. Ablations confirm
the real-time competency of CCO-VOXEL, while the actual experimental runs
indicate its robustness and efficacy for on-board implementations.

Benchmark Comparisons: We benchmark our CCO-VOXEL against fast planner [3] and
bounding volume [13, 14] algorithms. The comparison is done on two highly
cluttered simulation environments named Box Cylinder world and  Wall Grid
World (cf. Fig. 6) of dimensions $30\times 30\times 7$ and $7\times 30\times
7$ respectively. Over $200$ trials were conducted and the results presented in
Table I are the cumulative statistics of these trials in both the testing
environments. Maximum velocity and acceleration were set to $2m/s$ and
$3m/s^{2}$ respectively. We sample $50$ random trajectories at every iteration
of CEM. We implemented CCO-VOXEL in C++11, and used the PX4 SITL (Software in
the Loop) for simulations. During simulation, a Gaussian noise of zero mean
and $0.2m$ variance is added into the point cloud for testing purposes. Note
that although the point cloud noise is Gaussian, its mapping to the
distribution of the closest obstacle will be non-parametric.

The Fast-Planner [3] is a state-of-the-art algorithm known to generate optimal
and safe trajectories. However, in a cluttered and noisy environment the
distance to obstacle measurements are no longer reliable. Furthermore, the
gradients are noisy and can push the trajectory towards the obstacle leading
to a collision (cf. Fig. LABEL:teaser(a) and LABEL:teaser(b)) and a drop in
the success rate and increase in the smoothness cost of the trajectory. On
other hand, CCO-VOXEL is gradient free, and the MMD formulation takes into
account the noise uncertainty in the collision distance (cf. Fig.
LABEL:teaser(c)). As a result, it achieves a higher success rate and a lower
smoothness cost.

A possible approach to account for uncertainty is to construct a bounding
volume around the obstacle, which is obtained by increasing the robot’s
footprint by the size of the covariance of the distance distribution. This
technique has been used in several approaches such as [13]. Due to the
increase in a safety margin, the bounding volume approaches have a higher
accuracy in comparison to the fast planner [3], but they tend to generate
conservative trajectories. Furthermore, this method prevents the planner from
navigating in narrow spaces where it is difficult to maintain an appropriate
safe margin after the covariance-based inflation.

As shown in Table I the proposed method out performs both the benchmarks by a
fair margin. The results of our experiments demonstrate a 33.17% and 44.94%
decrease in the overall trajectory smoothness cost in comparison to Fast-
Planner and bounding volume based approaches respectively. Furthermore, the
success rate of CCO-VOXEL is approximately 3 times higher in comparison to the
benchmarked approaches. But our method has a slightly higher computation time,
due to the relatively high computation requirement of the evaluating the MMD
cost in CEM.

Real World Autonomous Flight We have conducted a fully autonomous flight test
in unknown and cluttered environment. The demonstration results can be found
here111https://github.com/sudarshan-s-harithas/CCO-VOXEL/tree/main/Demo. We
used a customized self developed quadrotor platform, which is equipped with an
Intel Realsense D415 depth camera [24], a Pixhawk flight controller is used
along with the PX4 firmware, all the computation modules which include
mapping, state estimation and motion planning is conducted online on a Intel
NUC 10 board [25] with Core i7-10710U Processor,1.1 GHz – 4.7 GHz Turbo, 6
core, the system is equipped with a 16GB RAM and 1TB SSD.

TABLE I: Benchmark Comparison Method | Computation | Smoothness | Success
---|---|---|---
| Time(s) | ($m^{2}/s^{5}$) | %
| mean: $0.0853$ | $\mathbf{12.45}$ | $\mathbf{100}$
Proposed | std: $0.014$ |  |
| mean: $0.012$ | $18.63$ | $32.37$
FastPlanner | std: $0.026$ |  |
| mean: $0.007$ | $22.54$ | $36.76$
Bounding Volumes | std: $0.001$ |  |

(a) Iteration 0 : Initilization

(b) Iteration 5

Figure 4: With every iteration of our CEM based trajectory optimization we
observe that the distribution of the collision-constraint violation approaches
the Dirac delta function and therefore the MMD cost is minimized. Note that
since plotting is based on approximate kernel density estimation from finite
samples, a tiny part of distribution appears to the left of 0 as well.

###  V-A Convergence analysis

Cross Entropy Minimization: In this section we demonstrate the ability of the
CEM to compute the minimum of our unconstrained reformulation of CCO (6). The
CEM takes the output of the initial trajectory search and further refines it.
The convergence of CEM is validated in Fig. 5. The figure presents the cost
profile of the mean trajectory after every iteration; it can be seen that the
optimal solution can be obtained in around four iterations with an average
computation time less than $0.015s$. The convergence can be verified by
observing the variance updates after every iteration. A monotonic decrease in
variance implies that after every iteration it converges to low-cost region
and thus it is more useful to sample trajectories close to mean trajectory.

Convergence of the Constraint Violation: We demonstrate that the convergence
of CEM (and minimization of MMD) implies a greater match between the
distribution of collision constraint violation and the Dirac Delta function.
This in turn leads to an increased probability of collision avoidance (recall
Remark 2). We observe from Fig. 4 that five iterations were enough for the
constraint violation distribution to converge to the Dirac distribution. This
convergence results have been empirically verified over a large number of
trials.

Figure 5: CEM Optimization: We observe that the normalized cost and the trace
of the covariance has monotonically decreased and the trajectory has converged
to an optimal solution

###  V-B Ablation Study

To study the impact of each component of our approach, we conducted an
ablation study between the possible methods.

MMD Computation: In Sec. IV-C we described the use of an auto-encoder
architecture to decrease the sample complexity and consequently time required
for an MMD computation. As presented in Table II we compare this approach
against the baseline method where RKHS embedding (Eq. 12) are directly
computed for a large number of samples without any low dimensional projection.
It can be observed that the use of an auto-encoder decreases the computation
time by approximately 12 times. Furthermore, we observe from Fig.4 that the
increase in computation time did not compromise the performance, we were still
able to minimize the probability of constraint violations within a very few
iterations of CEM.

TABLE II: Comparison between MMD computation methods  Method | Computation Time(s)
---|---
Autoencoder Embedding | $\mathbf{0.00012}$
Traditional Method | 0.00146

Initial Trajectory Searching Algorithm: Here, we validate the effectiveness of
using MMD as the edge-cost in A* based initial trajectory search. To this end,
we compare it against the standard kinodynamic A* used in [3]. As shown in
Table III, while our inclusion of MMD as edge cost increases the computation
time, it comes up with a massive 2.82 times improvement in success rate. It is
worth recalling that our computation time is higher because MMD is computed
over the several samples of closest distance to the obstacle drawn from a
black-box distribution.

TABLE III: Initial Trajectory Searching Comparison Method | Computation Time(s) | Success %
---|---|---
Kinodynamic A* with MMD | 0.0798 | $\mathbf{100}$
Kinodynamic A* [3] | 0.017 | 35.4
Figure 6: Gazebo Environment (a)  Wall Grid World (b) Box Cylinder World

##  VI Conclusions and Future Work

In this paper, we made a fundamental contribution towards bringing the
benefits of CCO-based safe trajectory planning to real-world settings commonly
encountered in quadrotor navigation. Our CCO formulation works directly on the
voxel grid representation of the world, wherein the collision constraints are
not known in analytical form but rather in the form of a black-box query
function. Existing approaches that assume known parametric nature in the
underlying uncertainty and geometric collision avoidance models are not
equipped to work under such minimalist assumptions. We outperformed state-of-
the-art deterministic planners in success rate and smoothness metric. In the
future, we aim to extend our work to handle state uncertainty in both static
and dynamic environments.

## References

  * [1] A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “Octomap: An efficient probabilistic 3d mapping framework based on octrees,” _Autonomous robots_ , vol. 34, no. 3, pp. 189–206, 2013.
  * [2] H. Oleynikova, Z. Taylor, M. Fehr, R. Siegwart, and J. Nieto, “Voxblox: Incremental 3d euclidean signed distance fields for on-board mav planning,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 1366–1373.
  * [3] B. Zhou, F. Gao, L. Wang, C. Liu, and S. Shen, “Robust and efficient quadrotor trajectory generation for fast autonomous flight,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3529–3536, 2019.
  * [4] B. Zhou, J. Pan, F. Gao, and S. Shen, “Raptor: Robust and perception-aware trajectory replanning for quadrotor fast flight,” _IEEE Transactions on Robotics_ , 2021.
  * [5] V. Usenko, L. Von Stumberg, A. Pangercic, and D. Cremers, “Real-time trajectory replanning for mavs using uniform b-splines and a 3d circular buffer,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 215–222.
  * [6] B. Zhou, F. Gao, J. Pan, and S. Shen, “Robust real-time uav replanning using guided gradient-based optimization and topological paths,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 1208–1214.
  * [7] L. Blackmore, H. Li, and B. Williams, “A probabilistic approach to optimal robust path planning with obstacles,” in _2006 American Control Conference_. IEEE, 2006, pp. 7–pp.
  * [8] B. Gopalakrishnan, A. K. Singh, and K. M. Krishna, “Closed form characterization of collision free velocities and confidence bounds for non-holonomic robots in uncertain dynamic environments,” in _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 4961–4968.
  * [9] B. Gopalakrishnan, A. K. Singh, M. Kaushik, K. M. Krishna, and D. Manocha, “Prvo: Probabilistic reciprocal velocity obstacle for multi robot navigation under uncertainty,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 1089–1096.
  * [10] H. Zhu and J. Alonso-Mora, “Chance-constrained collision avoidance for mavs in dynamic environments,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 776–783, 2019.
  * [11] S. N. J. Poonganam, B. Gopalakrishnan, V. S. S. B. K. Avula, A. K. Singh, K. M. Krishna, and D. Manocha, “Reactive navigation under non-parametric uncertainty through hilbert space embedding of probabilistic velocity obstacles,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2690–2697, 2020.
  * [12] B. Gopalakrishnan, A. K. Singh, K. M. Krishna, and D. Manocha, “Solving chance-constrained optimization under nonparametric uncertainty through hilbert space embedding,” _IEEE Transactions on Control Systems Technology_ , 2021.
  * [13] D. Claes, D. Hennes, K. Tuyls, and W. Meeussen, “Collision avoidance under bounded localization uncertainty,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 1192–1198.
  * [14] J. Snape, J. Van Den Berg, S. J. Guy, and D. Manocha, “The hybrid reciprocal velocity obstacle,” _IEEE Transactions on Robotics_ , vol. 27, no. 4, pp. 696–706, 2011.
  * [15] F. Gao and S. Shen, “Quadrotor trajectory generation in dynamic environments using semi-definite relaxation on nonconvex qcqp,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 6354–6361.
  * [16] D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in _2011 IEEE international conference on robotics and automation_. IEEE, 2011, pp. 2520–2525.
  * [17] C. Richter, A. Bry, and N. Roy, “Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments,” in _Robotics research_. Springer, 2016, pp. 649–666.
  * [18] L. von Stumberg, V. Usenko, J. Engel, J. Stückler, and D. Cremers, “Autonomous exploration with a low-cost quadrocopter using semi-dense monocular slam,” _arXiv preprint arXiv:1609.07835_ , 2016.
  * [19] J. J. Johnson and M. C. Yip, “Chance-constrained motion planning using modeled distance-to-collision functions,” _arXiv preprint arXiv:2107.10953_ , 2021\.
  * [20] M. Castillo-Lopez, P. Ludivig, S. A. Sajadi-Alamdari, J. L. Sanchez-Lopez, M. A. Olivares-Mendez, and H. Voos, “A real-time approach for chance-constrained motion planning with dynamic obstacles,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3620–3625, 2020.
  * [21] S. Boyd, “Stochastic programming,” _Lecture Notes, Stanford University_ , 2008\.
  * [22] R. J. Moss, “Cross-entropy method variants for optimization,” _arXiv preprint arXiv:2009.09043_ , 2020.
  * [23] C. Pinneri, S. Sawant, S. Blaes, J. Achterhold, J. Stueckler, M. Rolinek, and G. Martius, “Sample-efficient cross-entropy method for real-time planning,” _arXiv preprint arXiv:2008.06389_ , 2020.
  * [24] I. RealSense, “Intel realsense.” [Online]. Available: https://www.intelrealsense.com/depth-camera-d415/
  * [25] I. NUC, “Intel nuc.” [Online]. Available: https://www.intel.in/content/www/in/en/products/details/nuc.html


# On-line Search History-assisted Restart Strategy for Covariance Matrix
Adaptation Evolution Strategy ††thanks: This research was supported in part by
a grant from the Hong Kong Research Grants Council under GRF Grant CityU
125313 and GRF Grant 11200317, the National Natural Science Foundation of
China (Project No. 61603275, 61601329), and the Tianjin Higher Education
Creative Team Funds Program.

Yang Lou, Shiu Yin Yuen, Guanrong Chen  Department of Electronic Engineering
City University of Hong Kong
Hong Kong, China
<EMAIL_ADDRESS>
{kelviny.ee<EMAIL_ADDRESS>Xin Zhang  Tianjin Key Laboratory of
Wireless
Mobile Communications and Power Transmission
Tianjin Normal University
Tianjin, China
<EMAIL_ADDRESS>
###### Abstract

Restart strategy helps the covariance matrix adaptation evolution strategy
(CMA-ES) to increase the probability of finding the global optimum in
optimization, while a single run CMA-ES is easy to be trapped in local optima.
In this paper, the continuous non-revisiting genetic algorithm (cNrGA) is used
to help CMA-ES to achieve multiple restarts from different sub-regions of the
search space. The CMA-ES with on-line search history-assisted restart strategy
(HR-CMA-ES) is proposed. The entire on-line search history of cNrGA is stored
in a binary space partitioning (BSP) tree, which is effective for performing
local search. The frequently sampled sub-region is reflected by a deep
position in the BSP tree. When leaf nodes are located deeper than a threshold,
the corresponding sub-region is considered a region of interest (ROI). In HR-
CMA-ES, cNrGA is responsible for global exploration and suggesting ROI for
CMA-ES to perform an exploitation within or around the ROI. CMA-ES restarts
independently in each suggested ROI. The non-revisiting mechanism of cNrGA
avoids to suggest the same ROI for a second time. Experimental results on the
CEC 2013 and 2017 benchmark suites show that HR-CMA-ES performs better than
both CMA-ES and cNrGA. A positive synergy is observed by the memetic
cooperation of the two algorithms.

###### Index Terms:

non-revisiting genetic algorithm, covariance matrix adaptation evolution
strategy, memetic algorithm, on-line history

##  I Introduction

In the field of evolutionary computation, evaluation of solutions is
considered as the most time consuming and computational resource consuming
process. A well-designed evolutionary algorithm (EA) should cost negligible
computational resource itself, compared to the cost of solution evaluations.
Re-evaluation of any solution in EAs is wasteful and distorts the true
performance of the EA[1].

In the real-valued problems, since the search space is defined in the
continuous domain, the probability of generating two exactly the same
solutions is rare. However, EAs do not search purely randomly in the space,
and certain strategies in EAs may incur a risk of revisits. For example, it is
observed that the absorbing scheme of boundary handling method may generate
revisiting solutions, especially when the global/local optima is on or near a
boundary[2]. The absorbing scheme sets the variable to the boundary value if
it is generated outside the boundary. It is commonly adopted by artificial bee
colony (ABC)[3], differential evolution (DE)[4], and particle swarm
optimization (PSO)[5, 6], etc, and thus, a non-negligible number of revisits
on boundaries may exist in these EAs.

The non-revisiting genetic algorithm (NrGA)[7] records the entire on-line
search history by a binary space partitioning (BSP) tree, thus, re-evaluation
of any solution is completely avoided. The continuous version cNrGA (where ‘c’
stands for continuous) inherits the BSP tree to store the entire search
history[8, 9]. Crossover is used in cNrGA for the first-stage searching, which
introduces no new gene into the population, and thus identical solutions
(revisits) are prone to be generated by recombinations. Selection further
increases the probability of revisiting by giving the crossover opportunity to
partially fit genes. Whenever a revisit happens, the chromosome (solution)
would skip the re-evaluation, but instead, an adaptive mutation within a
certain sub-region is performed. The adaptive mutation introduces new gene(s)
into the population. More than just to avoid revisiting, the BSP tree archive
naturally offers parameter-less adaptive mutation operations for local search.

The cNrGA with a proper memory management makes efficiently use of the entire
on-line search history; cNrGA with least recently used (LRU) pruning (cNrGA-
LRU) prunes the LRU leaves in the BSP tree, if the memory threshold is reached
but the optimization is not terminated, and thus keeping the overall memory
use constant. It is shown that cNrGA-LRU can maintain the performance of
cNrGA, up to the empirical limit when 90% of the search history is missing[2].

The covariance matrix adaptation evolution strategy (CMA-ES)[10, 11] is a
stochastic derivative-free optimizer which specializes in ill-conditioned
local search. Though there is negligible or no revisit occurring in the search
process of CMA-ES by nature[2], it is an optimizer that would frequently be
trapped in local optima (premature convergence). An efficient remedy strategy
for CMA-ES being trapped in local optima is to restart it[12]. The niching
strategy[13, 14] for CMA-ES can be considered as a restart strategy that
focuses on performing parallel (re)start, while a sequential restart strategy
could consider when, where, and how. CMA-ES with restart strategies achieved
great success in both CEC and BBOB competitions[12].

In this paper, a memetic algorithm of cNrGA and CMA-ES is proposed, namely
CMA-ES with history-assisted restart strategy (HR-CMA-ES). HR-CMA-ES uses
cNrGA as the global exploration method, and CMA-ES as the local exploitation
method. The cooperation of the two algorithms bring benefit to each algorithm:
It provides a powerful local optimizer for cNrGA, and it also provides CMA-ES
good suggestions on region of interest (ROI) to restart, based on the search
history.

On the other hand, instead of pruning a single memory unit as in cNrGA-LRU,
HR-CMA-ES prunes a whole sub-region (i.e., ROI) if it has been exploited by
CMA-ES. Thus, HR-CMA-ES will not suggest the same ROI for CMA-ES to exploit,
but it does not forbid CMA-ES to exploit within the pruned sub-region again.
The ROI suggested by cNrGA has two features: 1) the sub-region is actively
visited by cNrGA; and 2) the sub-region has never been suggested before.

The rest of the paper is organized as follows: Sec. II reviews cNrGA and CMA-
ES. Sec. III introduces the proposed HR-CMA-ES. Experimental study with
comparisons and discussions, is given in Sec. IV. Conclusions are drawn in
Sec. V.

##  II Preliminaries

###  II-A cNrGA

Search history is accumulatively stored at each generation in cNrGA: The BSP
tree expands and the whole search space $S$ is partitioned. Each evaluated
solution is represented by a leaf node in the BSP tree, which occupies a sub-
region in the partitioned space. Suppose the current number of evaluated
solutions is $M$, then $S$ is divided into $M$ sub-regions. For example, two
solutions partition the search space into two parts, and the partitioning line
is drawn on the dimension along which the two solutions have the maximum
difference. Figs. 1, 2 and 3 show an example of BSP tree expansion and search
space partitioning in a 2-dimensional space.

Figure 1: Two solutions $a$ and $b$ have been evaluated and inserted into the
BSP tree. $R$ denotes the root node. Integers on the right denote the depth of
tree nodes. Figure 2: Four solutions have been evaluated and inserted in the
BSP tree, where $a$ becomes a virtual node of $c$; and $b$ becomes a virtual
node of $e$. Figure 3: Search space partitioning: (a) two solutions $a$ and
$b$ partition the search space into 2 sub-regions; (b) four solutions
partition the space into 4 parts.

Figure 1 shows when there are two evaluated solutions $a$ and $b$, the BSP
tree has two leaf nodes at depth 1 (depth 0 represents the root node).
Accordingly, in Fig. 3 (a) the 2-dimensional space is divided into two
regions. Since the two solutions $a$ and $b$ have a greater difference along
the dimension $x_{1}$ than $x_{2}$ ($\Delta x_{1}=|x_{1}^{a}-x_{1}^{b}|>\Delta
x_{2}=|x_{2}^{a}-x_{2}^{b}|$), the partitioning line is drawn at the center of
$x_{1}^{a}$ and $x_{1}^{b}$, on the dimension $x_{1}$.

As the number of evaluated solutions increases, more nodes are inserted into
the BSP tree as shown in Fig. 2; and the search space $S$ is further divided
into smaller regions as shown in Fig. 3 (b). There are now four evaluated
solutions as shown in Fig. 2, which are all leaf nodes located in depth 2.
Nodes $a$ and $c$ store the same evaluated solution, but node $c$ has a
smaller sub-region size than node $a$ in the search space (i.e., nodes $c$ and
$d$ divide the sub-region of $a$). Node $a$ becomes a virtual node of node
$c$, which is an intermediate node to reach other leaf nodes.

###  II-B CMA-ES

At each iteration of CMA-ES, new solutions are sampled according to a
multivariate normal distribution, as shown in Eq. 1, where $x_{i,g+1}$ denotes
the $i$th solution at generation $g+1$. Solutions with good fitness values
survive and are used to update the center of the multivariate normal
distribution ($m_{g}$) for the next generation sampling. The movement of
center (towards region with good fitness) is used to update the evolution
paths, including the isotropic and anisotropic path. The isotropic evolution
path is then used to update the step size ($\sigma_{g}$); the anisotropic path
is used to update the covariance matrix ($C_{g}$).

$x_{i,g+1}\sim\mathcal{N}(m_{g},\sigma_{g}^{2}C_{g})\sim
m_{g}+\sigma_{g}\times\mathcal{N}(0,C_{g})$ (1)

CMA-ES is particularly powerful if the problem is ill-conditioned, but is weak
if the problem landscape is rugged, when the covariance matrix may be trapped
in a local optima. As the searching process converges, the region near the
center of the multivariate normal distribution ($m_{g}$) has a much higher
probability of being sampled, while the region far away from the center may
not be sufficiently explored. Meanwhile, CMA-ES itself does not have a
mechanism to escape from local optima. However, the probability of CMA-ES
being stuck in local optima can be lowered when it is restarted, and thus the
restart strategy helps CMA-ES to become a good global optimizer, rather than a
local optimizer.

CMA-ES has developed a set of stopping criteria: Before the maximum number of
evaluations is used up, CMA-ES would stop, for example, if the condition
number of covariance matrix is large (e.g., greater than $1\times 10^{14}$),
or if the best objective values is not changed in the last several, say
$10+\lceil 30D/\lambda\rceil$ ($D$ and $\lambda$ represent problem dimension
and population size, respectively), generations. In addition, successively
increased population sizes in CMA-ES would benefit exploration[12]. The
stopping criteria could answer the when to restart question; the population
size modification suggests how to restart, but neither considers where to
restart.

##  III Search History-Based Restart Strategy

The proposed HR-CMA-ES is an attempt to answer the question of where to
restart CMA-ES. A preferable strategy is to restart from an unexplored
position, which requires a full memory of search history. The non-revisiting
scheme records the entire search history to avoid revisiting in the discrete
domain, and to perform adaptive mutation in both the discrete and continuous
domains. Therefore, the non-revisiting scheme is suitable to be employed to
suggest the region to (re)start CMA-ES.

As illustrated in Sec. II-A, a BSP tree of depth $\tau$ is capable of
recording $2^{\tau}$ solutions, only if the BSP tree is balanced. Practically,
the formed BSP tree is not balanced. Suppose there are $2^{\tau}$ solutions in
the BSP tree; many leaf nodes are located deeper than $\tau$, while many nodes
are shallower than $\tau$. The solutions that are deeper than $\tau$ suggest a
region of interest (ROI) for exploitation, while the shallow ones suggest
relatively less promising sub-regions. For example, the least recently used
solutions in cNrGA-LRU[2] are found amongst the shallow solutions in the BSP
tree.

Figure 4: Once any leaf node reaches a depth of $lv+k$, then CMA-ES is
triggered: The traced-back node on depth $lv$ (i.e., node $j$) is regarded as
the root of the sub-tree, and all the leaf nodes under $j$ form the initial
population of CMA-ES.

Suppose the solutions are allocated in a balanced manner in the BSP tree.
Given the maximum number of fitness evaluations is $M_{e}$, the depth of a BSP
tree to record the entire search history is $lv=\lceil\log_{2}{M_{e}}\rceil$,
where $\lceil x\rceil$ represents the minimum integer that is greater than or
equal to the real number $x$. Thus, $lv$ is regarded as the average depth of
the BSP tree. Any solution that exploits deeper than $lv$ may suggest a more
promising sub-region than average for further exploitation. Fig. 4 shows an
example, where nodes $s$ and $t$ have reached the depth $lv+k$. Then, cNrGA
would be suspended and CMA-ES would be triggered next. Node $j$ is the root of
the sub-tree that contains nodes $s$ and $t$. All the leaf nodes under node
$j$ are considered solutions within an ROI, and employed as the initial
population for CMA-ES restart. The recommended population size ($\lambda$) of
CMA-ES is shown in Eq. 2, where $D$ represents the problem dimension.

$\lambda=4+\lfloor 3\times\ln{D}\rfloor$ (2)

In HR-CMA-ES, to collect approximately $\lambda$ solutions for the CMA-ES
initialization, a sub-tree with depth $k=\lceil\log_{2}{\lambda}\rceil$ under
the sub-root (at depth $lv$) is needed. Note that a sub-tree with depth $k$
may include $\lambda_{0}$ ($k+1\leq\lambda_{0}\leq\lambda$) solutions. If any
leaf node reaches the depth of $lv+k$, all the leaf nodes in the sub-tree
(with the sub-root at depth $lv$) are collected as the initial population of
CMA-ES, and meanwhile the running of cNrGA is suspended. As shown in Fig. 4,
node $j$ at depth $lv$ is the sub-root; all the leaf nodes under node $j$,
namely, nodes $q$, $r$, $o$, $s$, $t$, etc. are employed as the initial
population of CMA-ES restart.

Figure 5: The corresponding sub-region in the search space. The sub-region of
size: $x_{\sigma 1}\times x_{\sigma 2}$ is the sub-region for node $j$ (the
entire sub-region of the sub-trees); as it is further partitioned, each
divided piece represents a leaf node. Figure 6: After CMA-ES exploitation, the
entire sub-region is regarded a blocked region for cNrGA, thus, any solution
generated by cNrGA within the sub-region would be considered a revisit, and
discarded without evaluation. Instead, another random solution would be
generated outside the blocked region. Figure 7: When the tree node is blocked,
no further insertion is allowed under the node. In contrast, other regions
would be further exploited, if any of them becomes an ROI.

The mean of all these collected solutions is used as the initial center
position of CMA-ES. The initial step size is set to
$\sigma=0.3\times\max\\{x_{\sigma 1},x_{\sigma 2}\\}$. The CMA-ES takes over
the optimization until any one of its stopping criteria is met. After CMA-ES
is run, the sub-region is denoted a blocked region for cNrGA (see Figs. 6 and
7). Note that this sub-region is not blocked to CMA-ES, though CMA-ES may have
very low probability visiting it again. As illustrated in Fig. 8), CMA-ES is
allowed to sample within the sub-region that is blocked to cNrGA.

Figure 8: (a) The gray-shaded sub-region has been blocked to cNrGA. (b) The
cNrGA explores a new ROI, and then CMA-ES is restarted where. The gray-shaded
sub-region is not closed to CMA-ES. Figure 9: Pseudo code for HR-CMA-ES.

##  IV Experimental Study

###  IV-A Benchmark Comparison

The CEC 2013[15] and CEC 2017[16] bound constrained real-valued benchmark
suites are employed for testing. There are 28 problems in CEC 2013 suite
(denoted by f1301 to f1328), and 29 problems in CEC 2017 suite (denoted by
f1701, f1703 to f1730, where f1702 is excluded because it shows unstable
behavior[16]). The problem dimension is set to 10D and 30D. The maximum number
of problem evaluations is $1\times 10^{5}$ for 10D problems and $3\times
10^{5}$ for 30D problems. The number of independent runs is set to 30. The
significance test is performed by the Kruskal-Wallis test[17] with confidence
level $\alpha=0.05$.

Two algorithms, CMA-ES and cNrGA-LRU, are compared with HR-CMA-ES. All three
algorithms are implemented in Matlab. The source code of CMA-ES is downloaded
from [18] and cNrGA from[19]. The employed CMA-ES is a version 3.33 copy with
default parameter settings. Before the maximum number of evaluations is used
up, CMA-ES may stop and restart, if any of the termination criteria is met.
The recommended termination criteria are presented in Sec. B.3 of [11]. The
population size of cNrGA-LRU is 100 and the crossover rate is 0.5 as
recommended in[2]. The memory-pruning rate of cNrGA-LRU is 0.5, meaning that
half of the search history that is recognized as LRU information is discarded
during the optimization.

TABLE I: Comparison of the three EAs on CEC 2013 benchmarks. (+) represents the algorithm is superior to HR-CMA-ES, and (–) represents inferior to HR-CMA-ES. CEC 2013 | 10D | 30D
---|---|---
HR | CMA | LRU | HR | CMA | LRU
f1301 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1302 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1303 | 1.5 | 1.5 | 3 (–) | 2 | 1 | 3 (–)
f1304 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1305 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1306 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1307 | 2 | 1 | 3 (–) | 1 | 2 | 3 (–)
f1308 | 2 | 3 (–) | 1 (+) | 2 | 3 (–) | 1 (+)
f1309 | 1 | 2 | 3 (–) | 2 | 1 | 3 (–)
f1310 | 1.5 | 1.5 | 3 (–) | 1 | 2 | 3 (–)
f1311 | 2 | 3 (–) | 1 (+) | 2 | 3 (–) | 1 (+)
f1312 | 2 | 1 | 3 (–) | 2 | 1 | 3 (–)
f1313 | 2 | 1 | 3 (–) | 2 | 1 | 3 (–)
f1314 | 2 | 3 (–) | 1 (+) | 2 | 3 (–) | 1 (+)
f1315 | 1 | 2 | 3 (–) | 1 | 2 | 3 (–)
f1316 | 2 | 3 (–) | 1 | 2 | 3 (–) | 1 (+)
f1317 | 2 | 3 | 1 (+) | 2 | 3 | 1 (+)
f1318 | 1 | 2 | 3 (–) | 1 | 2 | 3 (–)
f1319 | 2 | 3 | 1 | 2 | 3 | 1 (+)
f1320 | 2 | 3 | 1 (+) | 2 | 3 (–) | 1 (+)
f1321 | 2 | 1 | 3 (–) | 2 | 1 | 3 (–)
f1322 | 2 | 3 (–) | 1 (+) | 2 | 3 (–) | 1 (+)
f1323 | 1 | 2 | 3 | 2 | 1 | 3 (–)
f1324 | 2 | 1 (+) | 3 (–) | 2 | 1 (+) | 3 (–)
f1325 | 2 | 1 | 3 (–) | 2 | 1 | 3 (–)
f1326 | 2 | 1 | 3 (–) | 3 | 1 (+) | 2
f1327 | 2 | 1 (+) | 3 (–) | 2 | 1 (+) | 3 (–)
f1328 | 2 | 1 | 3 | 2 | 1 (+) | 3 (–)
avg | 1.73 | 1.84 | 2.43 | 1.80 | 1.80 | 2.39
(+) |  | 2 | 6 |  | 4 | 8
(–) |  | 4 | 18 |  | 6 | 19

Table I shows the rank comparison of the three algorithms on solving the CEC
2013 problems. If two or more algorithms obtain the same result, then they
share the rank. For example, HR-CMA-ES and CMA-ES get the same result for
f1301, and thus they share the first and second ranks, with each algorithm
obtaining rank 1.5. Note that although the significance test is multi-sample
wise, the significance mark (+/–) is marked pair-wise for the algorithm
compared to HR-CMA-ES. For the 10D problems, both HR-CMA-ES and CMA-ES
outperform cNrGA-LRU significantly on solving the first 5 unimodal problems
(f1301 to f1305). For the 15 basic multimodal problems (f1306 to f1320), it is
observed that cNrGA-LRU performs bipolarly. Specifically, cNrGA-LRU
outperforms HR-CMA-ES significantly on solving 5 problems (f1308, f1311,
f1314, f1317, and f1320), while it is significantly worse than HR-CMA-ES on
solving 8 problems (f1306, f1307, f1309, f1310, f1312, f1313, f1315, and
f1318). CMA-ES performs significantly worse than HR-CMA-ES on solving 4
problems (f1308, f1311, f1314, and f1316), and without outperforming HR-CMA-ES
on any basic multimodal problem.

As for the last 8 composition problems (f1321 to f1328), CMA-ES is superior to
HR-CMA-ES on solving problems f1324 and f1327, but is inferior on f1322. In
contrast, cNrGA-LRU performs well on f1322, but badly on problems f1324 to
f1327. The overall performance of HR-CMA-ES on 10D problems is that 1) the
average rank obtained by HR-CMA-ES is the best over the three algorithms; and
2) HR-CMA-ES is inferior to CMA-ES and cNrGA-LRU on solving 2 and 6 problems
respectively, but is superior to them on 4 and 18 problems, respectively.
Similar ranks and significance test results are observed on solving the 30D
problems.

The detailed results of HR-CMA-ES and CMA-ES on solving CEC 2013 problems are
given in Table III (for 10D problems) and Table IV (for 30D problems).

TABLE II: Comparison of the three EAs on CEC 2017 benchmarks. (+) represents the algorithm is superior to HR-CMA-ES, and (–) represents inferior to HR-CMA-ES. CEC 2017 | 10D | 30D
---|---|---
HR | CMA | LRU | HR | CMA | LRU
f1701 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1703 | 1.5 | 1.5 | 3 (–) | 1.5 | 1.5 | 3 (–)
f1704 | 1.5 | 1.5 | 3 (–) | 2 | 1 | 3 (–)
f1705 | 2 | 1 | 3 (–) | 2 | 1 | 3
f1706 | 1 | 2 | 3 (–) | 2 | 3 | 1
f1707 | 2 | 1 | 3 (–) | 2 | 1 | 3 (–)
f1708 | 1 | 2 | 3 (–) | 1 | 2 | 3 (–)
f1709 | 1.5 | 1.5 | 3 (–) | 2 | 1 | 3 (–)
f1710 | 2 | 3 | 1 | 1 | 3 | 2
f1711 | 2 | 3 (–) | 1 | 1 | 2 (–) | 3 (–)
f1712 | 1 | 2 | 3 (–) | 2 | 1 | 3 (–)
f1713 | 1 | 2 | 3 (–) | 2 | 1 | 3 (–)
f1714 | 1 | 2 | 3 (–) | 1 | 2 | 3 (–)
f1715 | 1 | 2 | 3 (–) | 2 | 1 | 3 (–)
f1716 | 1 | 2 | 3 (–) | 1 | 2 | 3 (–)
f1717 | 1 | 2 | 3 | 1 | 2 (–) | 3 (–)
f1718 | 1 | 2 | 3 (–) | 1 | 2 | 3 (–)
f1719 | 2 | 1 | 3 (–) | 1 | 2 | 3 (–)
f1720 | 2 | 3 (–) | 1 (+) | 1 | 2 | 3
f1721 | 3 | 1 | 2 | 2 | 1 | 3 (–)
f1722 | 2 | 1 | 3 (–) | 1 | 3 (–) | 2 (–)
f1723 | 2 | 1 (+) | 3 (–) | 2 | 1 (+) | 3
f1724 | 2 | 1 (+) | 3 (–) | 3 | 1 (+) | 2
f1725 | 2 | 1 | 3 (–) | 1 | 2 | 3 (–)
f1726 | 2 | 1 | 3 (–) | 2 | 1 (+) | 3 (–)
f1727 | 1 | 2 (–) | 3 (–) | 1 | 2 | 3 (–)
f1728 | 2 | 1 | 3 (–) | 1 | 2 | 3 (–)
f1729 | 1 | 3 | 2 | 1 | 3 | 2
f1730 | 2 | 1 | 3 (–) | 1 | 2 | 3 (–)
avg | 1.59 | 1.69 | 2.72 | 1.48 | 1.72 | 2.79
(+) |  | 2 | 1 |  | 3 | 0
(–) |  | 3 | 23 |  | 3 | 22

Table II shows the rank comparison of the three algorithms on solving the CEC
2017 benchmark suite, 10D and 30D. As in the CEC 2013 comparison, both HR-CMA-
ES and CMA-ES outperform cNrGA-LRU on solving the 2 unimodal problems (f1701
and f1703). For the following 7 multimodal problems (f1704 to f1710), a
similar phenomenon as in CEC 2013 is observed, i.e., cNrGA-LRU either performs
the best (f1710-10D and f1706-30D), or performs the worst (the rest of
multimodal problems except f1710-30D). In contrast, HR-CMA-ES and CMA-ES
performs close to each other on these problems.

Problems f1711 to f1720 are hybrid problems, which are difficult for cNrGA-LRU
to handle. HR-CMA-ES performs quite well on solving these problems, not only
because it gains many rank one results, but also because it performs
significantly better than CMA-ES in 4 cases (f1711-{10D,30D}, f1717-30D, and
f1720-10D).

The last 10 problems in CEC 2017 suite (f1721 to f1730) are composition
problems. Though the ranks of HR-CMA-ES and CMA-ES are close to each other,
CMA-ES gains more significantly better results on solving these problems. CMA-
ES performs better on the 10D problems, while HR-CMA-ES performs better on the
30D problems.

In summary, HR-CMA-ES gains the best overall ranks on the CEC 2017 test.
Compared to CMA-ES, HR-CMA-ES gains one more superior results to CMA-ES than
inferior to CMA-ES. Compared to cNrGA-LRU, HR-CMA-ES is inferior to cNrGA-LRU
only once, but superior to cNrGA-LRU on 23 10D problems and on 22 30D
problems, respectively.

The detailed results of HR-CMA-ES and CMA-ES on solving CEC 2017 problems are
given in Table V (for 10D problems) and Table VI (for 30D problems).

###  IV-B Discussions

History-assisted restart strategy improves the performance of CMA-ES on
solving multimodal problems. For the 32 multimodal and hybrid problems (i.e.,
f1306 to f1320, f1704 to f1720), HR-CMA-ES performs significantly better than
CMA-ES on 6 10D problems (f1308, f1311, f1314, f1316, f1711, and f1720) and 7
30D problems (f1308, f1311, f1314, f1316, f1320, f1711, and f1717), without
performing significantly worse than CMA-ES.

In this paper, HR-CMA-ES is considered as a CMA-ES with search history-
assisted restart, i.e., cNrGA explores a suitable promising ROI for CMA-ES to
(re)start. Once CMA-ES has finished an exploitation in the ROI, it shuts down
the opportunity of cNrGA suggesting it again. On the other hand, HR-CMA-ES
could also be considered a novel cNrGA with CMA-ES as a powerful local search.
As soon as an ROI is explored by cNrGA, a more powerful local optimizer CMA-ES
is applied, instead of the adaptive mutation. Once exploited, the sub-region
is blocked to cNrGA, which is equivalent to a novel pruning strategy that
prunes a whole sub-region, rather than a single solution. It is a more
efficient pruning strategy in the continuous space. A positive algorithm
synergy[20] is observed in the cooperation of cNrGA and CMA-ES in this work.

Compared to other CMA-ES variants with restart, HR-CMA-ES suggests where to
(re)start, but is unable to suggest how to (re)start. For example, strategy
like population size modification is difficult to be applied in HR-CMA-ES,
since CMA-ES is (re)initiated in different sub-regions of the search space,
but the suitable population sizes for different sub-regions are difficult to
predict.

##  V Conclusions

The covariance matrix adaptation evolution strategy (CMA-ES) is prone to be
trapped in local optima if the problem landscape being solved is rugged.
Restart strategy turns out to be an efficient remedy to help CMA-ES escape
from being trapped. The continuous non-revisiting genetic algorithm (cNrGA)
utilizes the entire on-line search history to perform revisiting avoidance and
adaptive mutation. Since the question about when and how to restart CMA-ES has
been well studied, this paper attempts to answer where to restart. The
history-based restart CMA-ES (HR-CMA-ES) is proposed. cNrGA suggests region of
interest (ROI) for CMA-ES to perform exploitation. An ROI is a sub-region
where cNrGA generates solutions frequently. The non-revisiting mechanism
avoids suggesting the same ROI for CMA-ES repeatedly. Experimental results on
the CEC 2013 and 2017 benchmark suites show that HR-CMA-ES significantly
outperforms cNrGA with constant memory, and performs better than CMA-ES in
terms of average rank. Particularly, HR-CMA-ES shows powerful optimization
ability on solving multimodal problems. A positive synergy is observed by the
cooperation of the two algorithms, which makes HR-CMA-ES a better algorithm
than the two component algorithms.

## References

  * [1] D. H. Wolpert and W. G. Macready, “No free lunch theorems for optimization,” _IEEE Transactions on Evolutionary Computation_ , vol. 1, no. 1, pp. 67–82, 1997.
  * [2] Y. Lou and S. Y. Yuen, “Non-revisiting genetic algorithm with adaptive mutation using constant memory,” _Memetic Computing_ , vol. 8, no. 3, pp. 189–210, 2016, doi:10.1007/s12293-015-0178-6.
  * [3] B. Akay and D. Karaboga, “A modified artificial bee colony algorithm for real-parameter optimization,” _Information Sciences_ , vol. 192, pp. 120–142, 2012.
  * [4] A. H. Gandomi and X.-S. Yang, “Evolutionary boundary constraint handling scheme,” _Neural Computing and Applications_ , vol. 21, no. 6, pp. 1449–1462, 2012.
  * [5] W. Chu, X. Gao, and S. Sorooshian, “Handling boundary constraints for particle swarm optimization in high-dimensional search space,” _Information Sciences_ , vol. 181, no. 20, pp. 4569–4581, 2011.
  * [6] M. Zambrano-Bigiarini, M. Clerc, and R. Rojas, “Standard particle swarm optimisation 2011 at CEC-2013: A baseline for future PSO improvements,” in _IEEE Congress on Evolutionary Computation (CEC)_. IEEE, 2013, pp. 2337–2344.
  * [7] S. Y. Yuen and C. K. Chow, “A genetic algorithm that adaptively mutates and never revisits,” _IEEE Transactions on Evolutionary Computation_ , vol. 13, no. 2, pp. 454–472, 2009.
  * [8] C. K. Chow and S. Y. Yuen, “Continuous non-revisiting genetic algorithm with random search space re-partitioning and one-gene-flip mutation,” in _IEEE Congress on Evolutionary Computation (CEC)_. IEEE, 2010, pp. 1–8.
  * [9] ——, “Continuous non-revisiting genetic algorithm with overlapped search sub-region,” in _IEEE Congress on Evolutionary Computation (CEC)_. IEEE, 2012, pp. 1–8.
  * [10] N. Hansen, “The CMA evolution strategy: A comparing review,” in _Towards a New Evolutionary Computation_. Springer, 2006, pp. 75–102.
  * [11] ——, “The CMA evolution strategy: A tutorial,” _arXiv preprint arXiv:1604.00772_ , 2016.
  * [12] I. Loshchilov, “CMA-ES with restarts for solving CEC 2013 benchmark problems,” in _IEEE Congress on Evolutionary Computation (CEC)_. IEEE, 2013, pp. 369–376.
  * [13] M. Preuss, “Niching the CMA-ES via nearest-better clustering,” in _Genetic and Evolutionary Computation Conference_. ACM, 2010, pp. 1711–1718.
  * [14] O. M. Shir, M. Emmerich, and T. Bäck, “Adaptive niche radii and niche shapes approaches for niching with the CMA-ES,” _Evolutionary Computation_ , vol. 18, no. 1, pp. 97–126, 2010.
  * [15] J. Liang, B. Qu, P. Suganthan, and A. G. Hernández-Díaz, “Problem definitions and evaluation criteria for the CEC 2013 special session on real-parameter optimization,” Zhengzhou University, China and Nanyang Technological University, Singapore, Tech. Rep., 2013.
  * [16] N. Awad, M. Ali, J. Liang, B. Qu, and P. N. Suganthan, “Problem definitions and evaluation criteria for the CEC 2017 special session and competition on single objective real-parameter numerical optimization,” Nanyang Technological University, Singapore, Jordan University of Science and Technology, Jordan, Zhengzhou University, China, Tech. Rep., 2017.
  * [17] Y. Lou, S. Y. Yuen, and G. Chen, “Evolving benchmark functions using Kruskal-Wallis test,” in _Genetic and Evolutionary Computation Conference (GECCO)_. ACM, 2018, pp. 1337–1341.
  * [18] CMA-ES source code. [Online]. Available: http://cma.gforge.inria.fr/cmaes.m
  * [19] cNrGA source code. [Online]. Available: http://www.ee.cityu.edu.hk/ syyuen/Public/Code.html
  * [20] S. Y. Yuen, C. K. Chow, X. Zhang, and Y. Lou, “Which algorithm should I choose: An evolutionary algorithm portfolio approach,” _Applied Soft Computing_ , vol. 40, pp. 654–673, 2016.

TABLE III: The detailed results of HR-CMA-ES and CMA-ES on solving CEC 2013 10D problems. CEC 2013 | HR-CMA-ES | CMA-ES
---|---|---
best | worst | median | mean | std | best | worst | median | mean | std
f1301 | 0.00E+00 | 2.27E-13 | 2.27E-13 | 1.74E-13 | 9.78E-14 | 0.00E+00 | 2.27E-13 | 2.27E-13 | 1.52E-13 | 1.09E-13
f1302 | 0.00E+00 | 2.27E-13 | 2.27E-13 | 1.89E-13 | 8.62E-14 | 0.00E+00 | 2.27E-13 | 2.27E-13 | 1.97E-13 | 7.86E-14
f1303 | 2.27E-13 | 6.82E-13 | 2.27E-13 | 3.41E-13 | 1.55E-13 | 2.27E-13 | 9.09E-13 | 2.27E-13 | 3.71E-13 | 1.84E-13
f1304 | 0.00E+00 | 2.27E-13 | 2.27E-13 | 2.20E-13 | 4.15E-14 | 0.00E+00 | 2.27E-13 | 2.27E-13 | 1.97E-13 | 7.86E-14
f1305 | 4.55E-13 | 2.73E-12 | 9.09E-13 | 1.15E-12 | 5.38E-13 | 3.41E-13 | 2.73E-12 | 1.02E-12 | 1.13E-12 | 5.73E-13
f1306 | 0.00E+00 | 2.27E-13 | 1.14E-13 | 1.44E-13 | 6.63E-14 | 0.00E+00 | 2.27E-13 | 1.14E-13 | 1.14E-13 | 2.99E-14
f1307 | 4.40E-02 | 1.09E+01 | 4.15E-01 | 9.17E-01 | 2.02E+00 | 5.22E-03 | 4.55E+00 | 2.65E-01 | 7.73E-01 | 1.18E+00
f1308 | 2.03E+01 | 2.07E+01 | 2.05E+01 | 2.05E+01 | 1.18E-01 | 2.06E+01 | 2.11E+01 | 2.09E+01 | 2.08E+01 | 1.42E-01
f1309 | 2.88E-01 | 3.33E+00 | 1.67E+00 | 1.65E+00 | 8.26E-01 | 5.56E-01 | 3.34E+00 | 1.98E+00 | 1.94E+00 | 7.93E-01
f1310 | 0.00E+00 | 5.68E-14 | 5.68E-14 | 5.31E-14 | 1.44E-14 | 0.00E+00 | 1.14E-13 | 5.68E-14 | 5.68E-14 | 1.49E-14
f1311 | 1.14E-13 | 4.97E+00 | 1.99E+00 | 2.03E+00 | 1.02E+00 | 1.99E+00 | 7.96E+00 | 3.98E+00 | 4.48E+00 | 1.58E+00
f1312 | 9.95E-01 | 8.95E+00 | 4.97E+00 | 4.88E+00 | 1.82E+00 | 2.98E+00 | 6.96E+00 | 4.97E+00 | 4.58E+00 | 1.21E+00
f1313 | 1.17E+00 | 1.86E+01 | 8.64E+00 | 8.70E+00 | 4.19E+00 | 9.95E-01 | 1.72E+01 | 8.21E+00 | 7.95E+00 | 4.18E+00
f1314 | 3.66E+00 | 2.72E+02 | 9.50E+01 | 1.06E+02 | 7.23E+01 | 1.68E+01 | 8.18E+02 | 5.11E+02 | 4.82E+02 | 1.68E+02
f1315 | 1.02E+02 | 6.75E+02 | 4.27E+02 | 4.27E+02 | 1.60E+02 | 1.69E+02 | 7.72E+02 | 4.82E+02 | 4.51E+02 | 1.58E+02
f1316 | 3.69E-02 | 1.89E+00 | 9.80E-01 | 1.01E+00 | 6.20E-01 | 6.28E-02 | 3.69E+00 | 2.37E+00 | 2.04E+00 | 1.16E+00
f1317 | 6.20E+00 | 1.63E+01 | 1.25E+01 | 1.18E+01 | 2.61E+00 | 4.93E+00 | 1.88E+01 | 1.42E+01 | 1.35E+01 | 3.27E+00
f1318 | 2.73E+00 | 1.70E+01 | 1.39E+01 | 1.18E+01 | 4.21E+00 | 6.52E+00 | 1.84E+01 | 1.47E+01 | 1.45E+01 | 2.79E+00
f1319 | 1.88E-01 | 7.13E-01 | 5.11E-01 | 4.87E-01 | 1.26E-01 | 3.42E-01 | 7.37E-01 | 5.40E-01 | 5.15E-01 | 9.13E-02
f1320 | 2.63E+00 | 4.02E+00 | 3.62E+00 | 3.57E+00 | 3.56E-01 | 2.12E+00 | 4.34E+00 | 3.79E+00 | 3.69E+00 | 4.20E-01
f1321 | 2.00E+02 | 4.00E+02 | 2.00E+02 | 2.07E+02 | 3.66E+01 | 1.00E+02 | 2.00E+02 | 2.00E+02 | 1.57E+02 | 5.04E+01
f1322 | 4.10E+01 | 5.96E+02 | 2.39E+02 | 2.41E+02 | 1.28E+02 | 4.14E+02 | 8.78E+02 | 5.63E+02 | 6.00E+02 | 1.18E+02
f1323 | 1.36E+02 | 9.78E+02 | 5.72E+02 | 5.51E+02 | 1.77E+02 | 2.10E+02 | 9.03E+02 | 6.04E+02 | 5.62E+02 | 1.66E+02
f1324 | 2.00E+02 | 2.12E+02 | 2.07E+02 | 2.06E+02 | 3.77E+00 | 1.05E+02 | 1.37E+02 | 1.08E+02 | 1.11E+02 | 7.07E+00
f1325 | 2.00E+02 | 2.12E+02 | 2.03E+02 | 2.05E+02 | 4.09E+00 | 1.15E+02 | 2.10E+02 | 2.01E+02 | 1.99E+02 | 1.72E+01
f1326 | 1.01E+02 | 2.00E+02 | 1.06E+02 | 1.40E+02 | 4.67E+01 | 1.03E+02 | 1.14E+02 | 1.07E+02 | 1.08E+02 | 2.65E+00
f1327 | 3.00E+02 | 5.17E+02 | 3.02E+02 | 3.25E+02 | 5.62E+01 | 3.00E+02 | 3.09E+02 | 3.00E+02 | 3.01E+02 | 1.69E+00
f1328 | 3.00E+02 | 3.00E+02 | 3.00E+02 | 3.00E+02 | 3.04E-09 | 1.00E+02 | 3.00E+02 | 3.00E+02 | 2.67E+02 | 7.58E+01
TABLE IV: The detailed results of HR-CMA-ES and CMA-ES on solving CEC 2013 30D problems. CEC 2013 | HR-CMA-ES | CMA-ES
---|---|---
best | worst | median | mean | std | best | worst | median | mean | std
f1301 | 2.27E-13 | 4.55E-13 | 4.55E-13 | 3.79E-13 | 1.09E-13 | 2.27E-13 | 4.55E-13 | 2.27E-13 | 3.26E-13 | 1.15E-13
f1302 | 4.55E-13 | 6.82E-13 | 4.55E-13 | 5.08E-13 | 9.78E-14 | 4.55E-13 | 6.82E-13 | 4.55E-13 | 5.00E-13 | 9.25E-14
f1303 | 6.83E-10 | 2.44E+02 | 3.85E-05 | 8.31E+00 | 4.45E+01 | 1.43E-11 | 2.59E+00 | 3.14E-06 | 1.54E-01 | 4.93E-01
f1304 | 4.55E-13 | 6.82E-13 | 4.55E-13 | 5.46E-13 | 1.13E-13 | 2.27E-13 | 6.82E-13 | 4.55E-13 | 5.00E-13 | 1.10E-13
f1305 | 2.61E-12 | 3.11E-10 | 1.07E-11 | 3.14E-11 | 6.28E-11 | 2.73E-12 | 2.75E-10 | 1.34E-11 | 2.69E-11 | 5.12E-11
f1306 | 2.27E-13 | 3.41E-13 | 2.27E-13 | 2.46E-13 | 4.31E-14 | 2.27E-13 | 3.41E-13 | 2.27E-13 | 2.50E-13 | 4.63E-14
f1307 | 1.01E+00 | 2.42E+01 | 1.08E+01 | 1.10E+01 | 6.59E+00 | 3.51E+00 | 3.36E+01 | 1.23E+01 | 1.38E+01 | 7.11E+00
f1308 | 2.09E+01 | 2.11E+01 | 2.10E+01 | 2.10E+01 | 6.98E-02 | 2.11E+01 | 2.13E+01 | 2.13E+01 | 2.13E+01 | 6.23E-02
f1309 | 1.05E+01 | 1.70E+01 | 1.40E+01 | 1.38E+01 | 1.65E+00 | 9.72E+00 | 1.72E+01 | 1.39E+01 | 1.38E+01 | 1.78E+00
f1310 | 1.14E-13 | 7.40E-03 | 1.71E-13 | 2.47E-04 | 1.35E-03 | 5.68E-14 | 7.40E-03 | 1.71E-13 | 2.47E-04 | 1.35E-03
f1311 | 8.95E+00 | 2.19E+01 | 1.39E+01 | 1.48E+01 | 3.94E+00 | 1.89E+01 | 3.98E+01 | 3.08E+01 | 3.00E+01 | 4.74E+00
f1312 | 1.69E+01 | 4.12E+01 | 2.96E+01 | 3.00E+01 | 5.58E+00 | 1.89E+01 | 3.28E+01 | 2.84E+01 | 2.78E+01 | 3.59E+00
f1313 | 4.32E+01 | 8.67E+01 | 6.39E+01 | 6.47E+01 | 1.11E+01 | 1.48E+01 | 9.36E+01 | 6.06E+01 | 6.00E+01 | 1.52E+01
f1314 | 7.68E+02 | 1.66E+03 | 1.18E+03 | 1.22E+03 | 2.67E+02 | 2.21E+03 | 3.33E+03 | 2.72E+03 | 2.73E+03 | 2.71E+02
f1315 | 1.45E+03 | 3.06E+03 | 2.49E+03 | 2.41E+03 | 3.84E+02 | 1.55E+03 | 2.95E+03 | 2.44E+03 | 2.43E+03 | 3.38E+02
f1316 | 2.22E-01 | 3.75E+00 | 3.13E+00 | 3.01E+00 | 6.63E-01 | 2.80E+00 | 6.01E+00 | 4.82E+00 | 4.77E+00 | 8.83E-01
f1317 | 2.94E+01 | 5.11E+01 | 4.14E+01 | 4.09E+01 | 5.26E+00 | 2.42E+01 | 6.86E+01 | 5.52E+01 | 5.06E+01 | 1.19E+01
f1318 | 3.32E+01 | 7.84E+01 | 6.12E+01 | 5.81E+01 | 1.25E+01 | 3.06E+01 | 9.01E+01 | 6.19E+01 | 5.94E+01 | 1.37E+01
f1319 | 1.34E+00 | 2.98E+00 | 2.33E+00 | 2.23E+00 | 3.85E-01 | 1.43E+00 | 2.96E+00 | 2.32E+00 | 2.27E+00 | 3.33E-01
f1320 | 1.19E+01 | 1.37E+01 | 1.32E+01 | 1.31E+01 | 5.06E-01 | 1.32E+01 | 1.50E+01 | 1.44E+01 | 1.44E+01 | 4.37E-01
f1321 | 1.00E+02 | 2.00E+02 | 2.00E+02 | 1.90E+02 | 3.05E+01 | 1.00E+02 | 2.00E+02 | 2.00E+02 | 1.80E+02 | 4.07E+01
f1322 | 5.82E+02 | 2.83E+03 | 1.15E+03 | 1.19E+03 | 4.91E+02 | 1.87E+03 | 3.62E+03 | 2.98E+03 | 2.96E+03 | 4.17E+02
f1323 | 2.06E+03 | 3.38E+03 | 2.76E+03 | 2.82E+03 | 3.57E+02 | 1.78E+03 | 3.25E+03 | 2.82E+03 | 2.73E+03 | 4.05E+02
f1324 | 2.25E+02 | 2.46E+02 | 2.39E+02 | 2.39E+02 | 6.15E+00 | 1.33E+02 | 2.30E+02 | 2.11E+02 | 1.95E+02 | 3.43E+01
f1325 | 2.44E+02 | 2.60E+02 | 2.54E+02 | 2.54E+02 | 4.27E+00 | 2.43E+02 | 2.59E+02 | 2.52E+02 | 2.52E+02 | 4.26E+00
f1326 | 2.00E+02 | 3.37E+02 | 2.84E+02 | 2.67E+02 | 5.82E+01 | 1.32E+02 | 2.00E+02 | 2.00E+02 | 1.88E+02 | 2.25E+01
f1327 | 5.60E+02 | 7.28E+02 | 6.80E+02 | 6.70E+02 | 4.20E+01 | 4.00E+02 | 6.57E+02 | 5.72E+02 | 5.64E+02 | 6.44E+01
f1328 | 1.00E+02 | 3.00E+02 | 3.00E+02 | 2.87E+02 | 5.07E+01 | 1.00E+02 | 3.00E+02 | 1.00E+02 | 1.41E+02 | 8.10E+01
TABLE V: The detailed results of HR-CMA-ES and CMA-ES on solving CEC 2017 10D problems. CEC 2017 | HR-CMA-ES | CMA-ES
---|---|---
best | worst | median | mean | std | best | worst | median | mean | std
f1701 | 0.00E+00 | 1.42E-14 | 1.42E-14 | 1.23E-14 | 4.91E-15 | 0.00E+00 | 1.42E-14 | 1.42E-14 | 1.18E-14 | 5.39E-15
f1703 | 0.00E+00 | 5.68E-14 | 5.68E-14 | 4.36E-14 | 2.45E-14 | 0.00E+00 | 5.68E-14 | 5.68E-14 | 4.55E-14 | 2.31E-14
f1704 | 0.00E+00 | 5.68E-14 | 5.68E-14 | 5.49E-14 | 1.04E-14 | 0.00E+00 | 5.68E-14 | 5.68E-14 | 5.31E-14 | 1.44E-14
f1705 | 2.98E+00 | 8.95E+00 | 5.97E+00 | 5.87E+00 | 1.68E+00 | 1.99E+00 | 8.95E+00 | 4.97E+00 | 5.50E+00 | 1.58E+00
f1706 | 1.39E-05 | 1.38E-03 | 9.72E-05 | 2.11E-04 | 3.02E-04 | 2.68E-05 | 3.17E-03 | 2.15E-04 | 4.50E-04 | 6.88E-04
f1707 | 4.63E+00 | 1.70E+01 | 1.29E+01 | 1.20E+01 | 3.26E+00 | 4.61E+00 | 1.60E+01 | 1.33E+01 | 1.20E+01 | 3.60E+00
f1708 | 2.98E+00 | 6.96E+00 | 4.97E+00 | 4.88E+00 | 1.24E+00 | 2.98E+00 | 8.95E+00 | 4.97E+00 | 5.14E+00 | 1.48E+00
f1709 | 0.00E+00 | 1.14E-13 | 1.14E-13 | 1.06E-13 | 2.88E-14 | 0.00E+00 | 1.14E-13 | 1.14E-13 | 1.10E-13 | 2.08E-14
f1710 | 1.02E+01 | 5.25E+02 | 2.66E+02 | 2.74E+02 | 1.56E+02 | 5.22E+01 | 7.45E+02 | 4.35E+02 | 3.81E+02 | 1.59E+02
f1711 | 1.99E+00 | 2.19E+01 | 6.96E+00 | 9.06E+00 | 5.76E+00 | 3.98E+00 | 4.38E+01 | 1.19E+01 | 1.41E+01 | 8.43E+00
f1712 | 6.82E-13 | 4.40E+02 | 1.58E+02 | 1.79E+02 | 1.16E+02 | 1.14E-12 | 5.06E+02 | 1.61E+02 | 1.87E+02 | 1.52E+02
f1713 | 1.31E+01 | 2.22E+02 | 7.67E+01 | 8.27E+01 | 4.15E+01 | 2.80E+01 | 4.20E+02 | 1.06E+02 | 1.28E+02 | 8.90E+01
f1714 | 2.10E+01 | 5.18E+01 | 2.99E+01 | 3.15E+01 | 7.54E+00 | 2.20E+01 | 6.48E+01 | 3.64E+01 | 3.73E+01 | 1.11E+01
f1715 | 2.31E+00 | 8.47E+01 | 2.39E+01 | 2.75E+01 | 1.96E+01 | 2.00E+00 | 7.57E+01 | 2.75E+01 | 3.34E+01 | 2.14E+01
f1716 | 9.33E-01 | 1.22E+02 | 1.99E+00 | 1.05E+01 | 2.34E+01 | 1.08E+00 | 2.20E+02 | 2.19E+00 | 3.62E+01 | 6.50E+01
f1717 | 1.69E+01 | 5.52E+01 | 2.88E+01 | 3.05E+01 | 8.55E+00 | 2.10E+01 | 6.12E+01 | 3.05E+01 | 3.36E+01 | 1.03E+01
f1718 | 2.00E+01 | 1.02E+02 | 2.67E+01 | 3.36E+01 | 1.78E+01 | 2.03E+01 | 9.96E+01 | 3.28E+01 | 4.45E+01 | 2.46E+01
f1719 | 4.80E+00 | 3.19E+01 | 9.97E+00 | 1.31E+01 | 7.82E+00 | 4.19E+00 | 3.27E+01 | 1.14E+01 | 1.25E+01 | 6.51E+00
f1720 | 2.10E+01 | 6.01E+01 | 3.16E+01 | 3.38E+01 | 1.03E+01 | 2.20E+01 | 1.65E+02 | 4.94E+01 | 5.62E+01 | 3.04E+01
f1721 | 1.91E+02 | 2.11E+02 | 2.07E+02 | 2.06E+02 | 4.43E+00 | 1.00E+02 | 2.10E+02 | 2.07E+02 | 2.03E+02 | 1.95E+01
f1722 | 1.00E+02 | 1.00E+02 | 1.00E+02 | 1.00E+02 | 2.89E-11 | 1.42E-09 | 1.00E+02 | 1.00E+02 | 7.49E+01 | 3.76E+01
f1723 | 3.05E+02 | 3.14E+02 | 3.09E+02 | 3.09E+02 | 2.02E+00 | 3.05E+02 | 3.10E+02 | 3.07E+02 | 3.07E+02 | 1.31E+00
f1724 | 3.23E+02 | 3.40E+02 | 3.35E+02 | 3.35E+02 | 3.39E+00 | 1.00E+02 | 3.37E+02 | 3.31E+02 | 2.82E+02 | 7.61E+01
f1725 | 3.98E+02 | 3.98E+02 | 3.98E+02 | 3.98E+02 | 1.71E-01 | 1.00E+02 | 3.98E+02 | 3.98E+02 | 3.78E+02 | 7.56E+01
f1726 | 2.00E+02 | 1.19E+03 | 3.00E+02 | 3.26E+02 | 1.65E+02 | 2.00E+02 | 3.00E+02 | 3.00E+02 | 2.77E+02 | 4.30E+01
f1727 | 3.89E+02 | 3.96E+02 | 3.91E+02 | 3.92E+02 | 2.17E+00 | 3.90E+02 | 3.98E+02 | 3.95E+02 | 3.95E+02 | 1.59E+00
f1728 | 3.00E+02 | 5.84E+02 | 3.00E+02 | 4.27E+02 | 1.41E+02 | 3.00E+02 | 3.28E+02 | 3.00E+02 | 3.01E+02 | 5.05E+00
f1729 | 2.42E+02 | 3.08E+02 | 2.58E+02 | 2.64E+02 | 1.97E+01 | 2.37E+02 | 2.98E+02 | 2.67E+02 | 2.68E+02 | 1.68E+01
f1730 | 3.95E+02 | 8.18E+05 | 4.70E+02 | 5.49E+04 | 2.07E+05 | 3.95E+02 | 5.55E+02 | 4.50E+02 | 4.57E+02 | 5.16E+01
TABLE VI: The detailed results of HR-CMA-ES and CMA-ES on solving CEC 2017 30D problems. CEC 2017 | HR-CMA-ES | CMA-ES
---|---|---
best | worst | median | mean | std | best | worst | median | mean | std
f1701 | 1.42E-14 | 2.84E-14 | 2.84E-14 | 2.56E-14 | 5.78E-15 | 1.42E-14 | 2.84E-14 | 2.84E-14 | 2.61E-14 | 5.39E-15
f1703 | 5.68E-14 | 1.71E-13 | 1.14E-13 | 1.08E-13 | 2.73E-14 | 5.68E-14 | 1.71E-13 | 1.14E-13 | 1.17E-13 | 2.56E-14
f1704 | 1.14E-13 | 5.86E+01 | 2.27E-13 | 8.36E+00 | 2.01E+01 | 1.14E-13 | 5.86E+01 | 1.71E-13 | 4.31E+00 | 1.48E+01
f1705 | 2.65E+01 | 4.38E+01 | 3.20E+01 | 3.33E+01 | 4.01E+00 | 2.29E+01 | 3.88E+01 | 3.08E+01 | 3.09E+01 | 4.71E+00
f1706 | 9.26E-04 | 1.38E-01 | 1.55E-02 | 2.73E-02 | 3.30E-02 | 4.78E-04 | 5.49E-01 | 9.27E-03 | 4.26E-02 | 1.03E-01
f1707 | 3.00E+01 | 6.33E+01 | 5.25E+01 | 5.20E+01 | 8.59E+00 | 3.30E+01 | 6.51E+01 | 5.53E+01 | 5.19E+01 | 9.19E+00
f1708 | 2.19E+01 | 4.28E+01 | 3.26E+01 | 3.25E+01 | 5.64E+00 | 2.29E+01 | 4.08E+01 | 3.33E+01 | 3.30E+01 | 4.63E+00
f1709 | 1.14E-13 | 8.95E-02 | 3.41E-13 | 8.95E-03 | 2.73E-02 | 2.27E-13 | 8.95E-02 | 3.41E-13 | 5.97E-03 | 2.27E-02
f1710 | 1.65E+03 | 3.06E+03 | 2.48E+03 | 2.38E+03 | 3.75E+02 | 1.87E+03 | 3.12E+03 | 2.58E+03 | 2.53E+03 | 3.64E+02
f1711 | 4.08E+01 | 1.18E+02 | 7.31E+01 | 7.75E+01 | 1.78E+01 | 4.08E+01 | 1.47E+02 | 9.45E+01 | 9.64E+01 | 2.45E+01
f1712 | 1.64E+02 | 1.64E+03 | 1.02E+03 | 1.05E+03 | 4.12E+02 | 1.69E+02 | 1.71E+03 | 9.67E+02 | 9.90E+02 | 3.15E+02
f1713 | 7.13E+02 | 2.66E+03 | 1.51E+03 | 1.57E+03 | 4.98E+02 | 7.02E+02 | 2.42E+03 | 1.53E+03 | 1.54E+03 | 5.34E+02
f1714 | 7.06E+01 | 2.19E+02 | 1.33E+02 | 1.42E+02 | 3.92E+01 | 8.17E+01 | 2.25E+02 | 1.57E+02 | 1.54E+02 | 4.10E+01
f1715 | 1.12E+02 | 6.46E+02 | 2.90E+02 | 3.12E+02 | 1.43E+02 | 8.19E+01 | 5.78E+02 | 2.85E+02 | 2.91E+02 | 1.18E+02
f1716 | 1.24E+01 | 1.02E+03 | 3.81E+02 | 4.13E+02 | 2.24E+02 | 5.33E+01 | 8.53E+02 | 4.77E+02 | 4.77E+02 | 2.07E+02
f1717 | 5.54E+01 | 3.53E+02 | 1.13E+02 | 1.47E+02 | 8.09E+01 | 5.85E+01 | 4.04E+02 | 2.15E+02 | 2.16E+02 | 9.41E+01
f1718 | 6.00E+01 | 2.85E+02 | 1.59E+02 | 1.64E+02 | 6.20E+01 | 8.40E+01 | 3.30E+02 | 1.75E+02 | 1.72E+02 | 5.94E+01
f1719 | 5.92E+01 | 2.42E+02 | 1.38E+02 | 1.44E+02 | 4.56E+01 | 5.67E+01 | 3.01E+02 | 1.51E+02 | 1.55E+02 | 4.93E+01
f1720 | 1.12E+02 | 4.88E+02 | 2.42E+02 | 2.56E+02 | 1.05E+02 | 1.15E+02 | 5.14E+02 | 2.31E+02 | 2.67E+02 | 1.21E+02
f1721 | 2.19E+02 | 2.46E+02 | 2.35E+02 | 2.35E+02 | 6.95E+00 | 2.25E+02 | 2.44E+02 | 2.34E+02 | 2.34E+02 | 4.93E+00
f1722 | 1.00E+02 | 3.20E+03 | 1.00E+02 | 3.88E+02 | 8.84E+02 | 1.00E+02 | 3.20E+03 | 2.47E+03 | 1.93E+03 | 1.16E+03
f1723 | 3.75E+02 | 4.13E+02 | 3.93E+02 | 3.92E+02 | 9.03E+00 | 3.67E+02 | 3.99E+02 | 3.83E+02 | 3.83E+02 | 8.60E+00
f1724 | 4.46E+02 | 4.77E+02 | 4.61E+02 | 4.61E+02 | 7.23E+00 | 2.00E+02 | 4.58E+02 | 4.49E+02 | 4.40E+02 | 4.56E+01
f1725 | 3.83E+02 | 3.87E+02 | 3.83E+02 | 3.84E+02 | 1.55E+00 | 3.83E+02 | 3.87E+02 | 3.83E+02 | 3.84E+02 | 1.55E+00
f1726 | 1.10E+03 | 1.54E+03 | 1.44E+03 | 1.42E+03 | 9.25E+01 | 2.00E+02 | 1.41E+03 | 1.26E+03 | 1.09E+03 | 4.14E+02
f1727 | 5.01E+02 | 5.44E+02 | 5.15E+02 | 5.17E+02 | 1.04E+01 | 5.02E+02 | 5.36E+02 | 5.21E+02 | 5.19E+02 | 8.44E+00
f1728 | 3.00E+02 | 3.00E+02 | 3.00E+02 | 3.00E+02 | 1.49E-10 | 3.00E+02 | 3.00E+02 | 3.00E+02 | 3.00E+02 | 9.32E-11
f1729 | 5.16E+02 | 9.78E+02 | 7.35E+02 | 7.36E+02 | 1.19E+02 | 5.12E+02 | 9.67E+02 | 7.31E+02 | 7.38E+02 | 1.19E+02
f1730 | 2.12E+03 | 2.74E+03 | 2.41E+03 | 2.41E+03 | 1.76E+02 | 2.13E+03 | 2.87E+03 | 2.46E+03 | 2.49E+03 | 2.04E+02