Adaptive Ranking-Based Constraint Handling for Explicitly Constrained Black-Box Optimization

Sakamoto, Naoki; Akimoto, Youhei

doi:10.1162/evco_a_00310

Abstract

We propose a novel constraint-handling technique for the covariance matrix adaptation evolution strategy (CMA-ES). The proposed technique is aimed at solving explicitly constrained black-box continuous optimization problems, in which the explicit constraint is a constraint whereby the computational time for the constraint violation and its (numerical) gradient are negligible compared to that for the objective function. This method is designed to realize two invariance properties: invariance to the affine transformation of the search space, and invariance to the increasing transformation of the objective and constraint functions. The CMA-ES is designed to possess these properties for handling difficulties that appear in black-box optimization problems, such as non-separability, ill-conditioning, ruggedness, and the different orders of magnitude in the objective. The proposed constraint-handling technique (CHT), known as ARCH, modifies the underlying CMA-ES only in terms of the ranking of the candidate solutions. It employs a repair operator and an adaptive ranking aggregation strategy to compute the ranking. We developed test problems to evaluate the effects of the invariance properties, and performed experiments to empirically verify the invariance of the algorithm. We compared the proposed method with other CHTs on the CEC 2006 constrained optimization benchmark suite to demonstrate its efficacy. Empirical studies reveal that ARCH is able to exploit the explicitness of the constraint functions effectively, sometimes even more efficiently than an existing box-constraint handling technique on box-constrained problems, while exhibiting the invariance properties. Moreover, ARCH overwhelmingly outperforms CHTs by not exploiting the explicit constraints in terms of the number of objective function calls.

Explicit constraint, black-box optimization, invariance, CMA-ES

1 Introduction

We consider explicitly constrained black-box continuous minimization problems defined as

\underset{x \in R^{n}}{\arg \min} f (x) subject to g_{j} (x) ⩽ 0, \forall j = 1, \dots, m,

(1)

where

f : R^{n} \to R

is the objective function and

g_{j} : R^{n} \to R

(⁠

j = 1, \dots, m

⁠) are the constraint functions. The equality constraints

h_{k} (x) = 0, \forall k = 1, \dots, l

are assumed to be transformed into inequality constraints

g_{k} (x) = | h_{k} (x) | - ɛ_{eq}

⁠, with a numerical tolerance

ɛ_{eq} > 0

⁠. In the black-box optimization scenario, the evaluation of the objective function

f

often requires computationally expensive simulation, while the explicit constraints can be computed independently of

f

-calls, and their computational cost is significantly lower than that of

f

-calls. These can be presented in a relatively simple mathematical expression, and their gradients are available symbolically or can be estimated numerically with a relatively low computational cost. The most common example of an explicit constraint is the box constraint, in which each coordinate of the design variable

x

is constrained in a closed interval. Explicit constraints often appear in engineering optimization as prerequisites for executing the simulation to compute the objective function. Taking into account the abovementioned situations, in this study, we assume that

the computational time for $g_{j} (x)$ and its (numerical) gradient $\nabla g_{j} (x)$ are negligible compared to that for $f (x)$ ⁠; and
the objective function $f (x)$ is not necessarily defined for an infeasible solution $x$ ⁠, which violates some $g_{j} (x)$ ⁠.

The covariance matrix adaptation evolution strategy (CMA-ES) (Hansen and Ostermeier, 2001; Hansen et al., 2003; Hansen, 2016) is employed in this study as the baseline optimization method for solving the problem (1). The CMA-ES is known as one of the state-of-the-art zeroth-order optimization algorithms for unconstrained black-box continuous optimization. In particular, it has been demonstrated as efficient for difficult objective functions such as non-convex, ill-conditioned, and nonseparable functions. Since these difficulties of objective functions naturally appear in engineering optimization, whether or not constraints exist, it is desirable for search algorithms for constrained optimization to operate efficiently under these difficulties. Although the CMA-ES is designed for unconstrained optimization, it can be applied to constrained optimization with the aid of constraint-handling techniques (CHTs).

One key characteristic of the CMA-ES is its invariance to several transformations of the objective function and search space (Hansen et al., 2011). The invariance properties induce equivalent problem classes. All instances (problems) in an equivalent class are regarded as equivalent under the corresponding transformation of the initial search distribution. Owing to the invariance properties of the CMA-ES, if its initial state is transformed properly, it is empirically observed that the CMA-ES minimizes ill-conditioned and nonseparable quadratic functions as efficiently as well-conditioned and separable spherical functions. This is key to the success of the CMA-ES in real-world problems because it is designed to adapt the distribution to such a proper state. Although the invariance properties themselves do not imply algorithm efficacy, they are useful for generalizing observations. That is, they are essential for assessing the performance of algorithms empirically. The invariance further contributes to the quasi-parameter-free feature of the CMA-ES, which is an important characteristic of the CMA-ES that has attracted attention from practitioners. As opposed to many other evolutionary approaches to continuous black-box optimization, in which the hyper-parameters are required to be tuned depending on the problem characteristics to enable efficient performance (Karafotias et al., 2015), default values are prepared for all hyper-parameters of the CMA-ES, depending only on the search space dimension $n$ ⁠. As many advantages of the CMA-ES for unconstrained optimization originate from its invariance properties, we hypothesize that a CHT for explicitly constrained optimization is desirable to preserve the invariance as much as possible to use the CMA-ES for constrained optimization.

However, the CHTs for explicit constraints employed in variants of the ES are not designed to preserve the invariance of the baseline search algorithms, as we discuss briefly in Section 3. When the CMA-ES is applied to an explicitly constrained black-box optimization problem with these CHTs, it loses the invariance properties exhibited by the CMA-ES. Certain CHTs, such as death penalty or resampling techniques, are applicable to the CMA-ES and preserve its invariance properties. However, as these approaches do not exploit the fact that the constraint violations are cheap to evaluate and their gradient information is available, they are often inefficient, as observed in Section 6.

We propose a novel CHT for the CMA-ES in solving the explicitly constrained optimization problem (1), named adaptive ranking-based constraint handling (ARCH). ARCH is designed to include the following two types of invariance properties: invariance to the element-wise increasing transformation of the objective and constraint functions, and invariance to the affine transformation of the search space. ARCH replaces the evaluation step in the sampling–evaluation–update cycle of the CMA-ES, as follows: A candidate solution generated by the CMA-ES is first repaired on a boundary of the feasible domain. The objective function value is evaluated at the repaired point. A penalty for the repair operation is computed by the Mahalanobis distance between the original and repaired solutions under the covariance matrix of the current search distribution. The candidate solutions are ranked based on the adaptive weighted sum of the rankings of the objective function values and rankings of the penalty values. An adaptive penalty coefficient is introduced to control the balance between the rankings, and is adapted for the search distribution so as not to move away from the feasible domain in terms of the Mahalanobis distance, while allowing infeasible but near boundary solutions to exhibit high rankings.

The contributions of this study are summarized as follows. Firstly, we present ARCH, which can handle explicit and nonlinear constraints.^¹ We prove that ARCH is invariant to any element-wise increasing transformation of $f$ and $g_{j}$ ⁠, and to any affine transformation of the search space coordinates. To the best of the authors' knowledge, this is the first approach for explicit constraints that is invariant to these transformations. Secondly, we empirically evaluate the effectiveness of the proposed approach from an invariance perspective. We develop test problems to demonstrate the effects of the invariance to affine transformation. In test cases, we empirically observe that ARCH is invariant to affine transformations, and illustrate that ARCH is even more effective for box constrained optimization problems than an existing CHT specialized for the box constraint. Thirdly, we compare ARCH with existing CHTs for non-explicit constraints on problems in which the objective function is defined for the infeasible domain. We use the CEC 2006 constrained optimization benchmark suite to demonstrate that ARCH overwhelmingly outperforms other CHTs for non-explicit constraints. This indicates that the explicit constraints should be treated as explicit constraints even if the objective function is defined on an infeasible domain. Compared to the previous publication in Sakamoto and Akimoto (2019), we $i$ ⁠) improve the algorithm to prevent the search distribution from being unnecessarily biased toward the boundary if the population size is larger than the default, $i i$ ⁠) prove the invariance properties of ARCH, and $i i i$ ⁠) compare ARCH with existing approaches on the CEC 2006 testbed.

The remainder of this article is organized as follows. We summarize the mathematical notations applied throughout the article below. Section 2 introduces the baseline optimization algorithm, namely the CMA-ES, and related CHTs. We discuss the invariance properties for constrained optimization problems in Section 3. Our proposed CHT, ARCH, is described in Section 4. The invariant properties of ARCH are demonstrated in Section 5. In Section 6, we empirically demonstrate how the invariance properties operate in practice, by means of numerical experiments on linearly constrained quadratic problems, to observe the efficacy of ARCH. We present our comparison of ARCH with other CHTs using the CEC 2006 constrained optimization testbed in Section 7. The article is concluded in Section 8.

Notations In the following, $R$ is the set of real numbers and $R^{+}$ is the set of strictly positive real numbers. Let $x \in R^{n}$ be an $n$ -dimensional column vector, where $x^{T}$ is its transpose, $∥ x ∥$ denotes the Euclidean norm of $x$ ⁠, and ${[x]}_{i}$ denotes the $i$ th coordinate of $x$ ⁠. Note that the $i$ th coordinate of the $k$ th vector $x_{k}$ is denoted by ${[x_{k}]}_{i}$ ⁠. The (⁠ $i$ ⁠, $j$ ⁠)th element of a matrix $A$ is also denoted by ${[A]}_{i, j}$ ⁠. The identity matrix is denoted by $I_{n}$ ⁠. The indicator function $1_{{condition}}$ returns 1 if condition is true, and 0 otherwise. The sign function $sgn (a)$ returns 1 if $a > 0$ ⁠, $- 1$ if $a < 0$ ⁠, and 0 otherwise. The integer interval between and including $a$ and $b$ is denoted by $[[a, b]]$ ⁠.

2 CMA-ES and Related CHTs

In this section, we introduce the CMA-ES, which is our baseline algorithm for unconstrained continuous optimization, followed by an overview of the existing CHTs for the CMA-ES.

2.1 CMA-ES

The CMA-ES (Hansen et al., 2003; Hansen, 2016) is a stochastic multi-point search algorithm for black-box continuous minimization of $f : R^{n} \to R$ ⁠. The CMA-ES samples $λ$ candidate solutions $x_{k}$ for $k \in {1, \dots, λ}$ from the multivariate normal distribution $N (m, σ^{2} C)$ ⁠, where $m \in R^{n}$ is the mean vector, $σ \in R^{+}$ is the step size, and $C \in R^{n \times n}$ is the covariance matrix. These distribution parameters are updated using the candidate solutions and their ranking information.

Step0. Initialize $m^{(0)}, σ^{(0)}, C^{(0)}$ according to the initial problem search domain, and initialize two evolution paths $p_{c}^{(0)} = p_{σ}^{(0)} = 0$ and their correction factors $γ_{σ}^{(0)} = γ_{c}^{(0)} = 0$ ⁠. All parameters that appear in the following are set to the default values listed in Table 1. The meanings of these parameters are described in Hansen (2016). These are designed based on theoretical ES research (e.g., see Akimoto et al., 2020) and extensive experiments. The CMA-ES repeats the following steps at each iteration, $t = 0, 1, \dots$ ⁠, until a termination criterion is satisfied.

Table 1:

Default parameter setting for CMA-ES.

λ = 4 + ⌊ 3 ln (n) ⌋, μ = ⌊ λ / 2 ⌋

w_{i} = \frac{ln (\frac{λ + 1}{2}) - ln (i)}{\sum_{k = 1}^{μ} (ln (\frac{λ + 1}{2}) - ln (k))}, c_{σ} = \frac{μ_{w} + 2}{n + μ_{w} + 5}, c_{c} = \frac{4 + μ_{w} / n}{n + 4 + 2 μ_{w} / n}

c_{1} = \frac{2}{{(n + 1.3)}^{2} + μ_{w}}, c_{μ} = min (1 - c_{1}, \frac{2 (μ_{w} - 2 + 1 / μ_{w})}{{(n + 2)}^{2} + μ_{w}})

d_{σ} = 1 + c_{σ} + 2 \times max (0, \sqrt{\frac{μ_{w} - 1}{n + 1}} - 1)

Step1. Draw $λ$ samples $z_{k}$ for $k \in {1, \dots, λ}$ independently from $N (0, I_{n})$ ⁠. Compute $y_{k} = \sqrt{C^{(t)}} z_{k}$ and $x_{k} = m^{(t)} + σ^{(t)} y_{k}$ ⁠. Then, $x_{k}$ (⁠ $k = 1, \dots, λ$ ⁠) are the candidate solutions that are independently $N (m^{(t)}, {(σ^{(t)})}^{2} C^{(t)})$ distributed. Here, $\sqrt{C^{(t)}}$ is the symmetric matrix satisfying $C^{(t)} = {(\sqrt{C^{(t)}})}^{2}$ ⁠.

Step2. Evaluate the candidate solutions $x_{k}$ ⁠, for $k \in {1, \dots, λ}$ ⁠, on the loss function $L$ ⁠, and sort them in ascending order. In an unconstrained optimization scenario, usually, $L = f$ ⁠. Let the $i$ th best candidate solution be denoted by $x_{i : λ}$ ⁠. In the same manner, we denote the corresponding steps and normalized steps as $y_{i : λ}$ and $z_{i : λ}$ ⁠, respectively.

Step3. Compute the weighted sum of the

μ

best steps of the candidate solutions

{〈 y 〉}_{w} = \sum_{i = 1}^{μ} w_{i} y_{i : λ}

and update the mean vector

m^{(t)}

⁠, as follows:

m^{(t + 1)} = m^{(t)} + σ^{(t)} {〈 y 〉}_{w},

where

w_{i}

is the recombination weight for the

i

th best candidate, which satisfies

w_{1} ⩾ w_{2} ⩾ \dots ⩾ w_{μ} > 0

and

\sum_{i = 1}^{μ} w_{i} = 1

⁠.

Step4. Update the evolution paths according to

\begin{matrix} p_{σ}^{(t + 1)} & = & (1 - c_{σ}) p_{σ}^{(t)} + \sqrt{c_{σ} (2 - c_{σ}) μ_{w}} {(\sqrt{C^{(t)}})}^{- 1} {〈 y 〉}_{w}, \\ p_{c}^{(t + 1)} & = & (1 - c_{c}) p_{c}^{(t)} + h_{σ}^{(t + 1)} \sqrt{c_{c} (2 - c_{c}) μ_{w}} {〈 y 〉}_{w}, \end{matrix}

where

c_{σ}

and

c_{c}

are the cumulation factors for the evolution paths,

μ_{w} = 1 / \sum_{i = 1}^{μ} w_{i}^{2}

⁠,

h_{σ}^{(t + 1)} = \{\begin{matrix} 1 & if ∥ p_{σ}^{(t + 1)} ∥ < (1.4 + \frac{2}{n + 1}) {(γ_{σ}^{(t + 1)})}^{\frac{1}{2}} χ \\ 0 & otherwise, \end{matrix}

and

χ = E [∥ N (0, I_{n}) ∥] \approx n^{\frac{1}{2}} (1 - \frac{1}{4 n} + \frac{1}{21 n^{2}})

is the expectation of the norm of the

n

variate standard normal distribution. The Heaviside function

h_{σ}^{(t + 1)}

stalls the update of

p_{c}^{(t + 1)}

if

∥ p_{σ}^{(t + 1)} ∥

is large. The correction factors for the evolution paths^² are updated as follows:

\begin{matrix} γ_{σ}^{(t + 1)} & = & {(1 - c_{σ})}^{2} γ_{σ}^{(t)} + c_{σ} (2 - c_{σ}), \\ γ_{c}^{(t + 1)} & = & {(1 - c_{c})}^{2} γ_{c}^{(t)} + h_{σ}^{(t + 1)} c_{c} (2 - c_{c}) . \end{matrix}

Step5. Update the step size and covariance matrix as follows:

\begin{matrix} σ^{(t + 1)} & = & σ^{(t)} exp (\frac{c_{σ}}{d_{σ}} (\frac{∥ p_{σ}^{(t + 1)} ∥}{χ} - {(γ_{σ}^{(t + 1)})}^{\frac{1}{2}})), \\ C^{(t + 1)} & = & C^{(t)} + c_{1} (p_{c}^{(t + 1)} {(p_{c}^{(t + 1)})}^{T} - γ_{c}^{(t + 1)} C^{(t)}) + c_{μ} \sum_{i = 1}^{μ} w_{i} (y_{i : λ} {(y_{i : λ})}^{T} - C^{(t)}), \end{matrix}

where

d_{σ}

is the damping parameter for the step size adaptation, while

c_{1}

and

c_{μ}

are the learning rate for the rank-one and rank-

μ

updates of the covariance matrix, respectively.

2.2 CHTs for ESs

We briefly review the CHTs employed in variants of ESs.

2.2.1 Resampling and Death Penalty

The resampling technique is the simplest CHT. Candidate solutions are resampled repeatedly until $λ$ feasible candidate solutions are generated. To guarantee that the resampling stops within a finite time, the maximum number of sampling in one iteration is set to a finite number. If the number of feasible candidate solutions is less than $λ$ ⁠, the current population is filled with infeasible solutions and the worst loss value is assigned to them (for example, $+ \infty$ ⁠). When the maximum sampling number is set to $λ$ ⁠, the resampling technique is known as the death penalty.

The resampling and death penalty methods are easy to implement and are applicable to any constraint type. However, they are not appropriate if the optimum is located on the boundary of the feasible domain, as candidate solutions are biased in the feasible domain and the search distribution tends to approach the boundary slowly, as we observe in Section 6. Moreover, the search algorithm cannot conduct a meaningful ranking of the candidate solutions if the probability of sampling feasible solutions is rather low and the population is filled with infeasible solutions. If this is the case, the parameter update will result in random fluctuation.

2.2.2 Penalty Function Methods

CHTs based on the penalization of infeasible candidate solutions are the most extensively used CHTs for real-world engineering optimization problems. The main concept of penalty function methods is transforming a constrained optimization problem into an unconstrained optimization problem by defining the penalized loss function:

L (x) = f (x) + p (x),

(2)

where

p

is a penalty function. If the objective function value is not well defined for infeasible solutions, a repair operation needs to be applied, and the first term is replaced with the objective function value evaluated at the repaired solution. The penalty function is manually designed in many engineering optimization problems, and a typical choice is the weighted sum of the constraint violations

max (g_{j} (x), 0)

or its monotone transformation, such as square.

The adaptive penalty box constraint handling (AP-BCH) method (Hansen et al., 2009) is a penalty function-based box CHT. The loss function is defined as

L (x) = f (x^{feas}) + \frac{1}{n} \sum_{i = 1}^{n} γ_{i} {({[x]}_{i} - {[x^{feas}]}_{i})}^{2},

(3)

where

γ_{i}

is an adaptive penalty coefficient and

x^{feas}

is a feasible vector closest to the infeasible solution

x

⁠; that is,

x^{feas} = {\arg \min}_{y} ∥ x - y ∥

⁠. The feasible solution

x^{feas}

is only used for evaluating the objective function and the penalty; that is, the repair operator is used in the Darwinian manner. This approach does not assume that the objective function is well defined outside the feasible domain, but it is only applicable to box-constrained problems.

The adaptive augmented Lagrangian constraint handling (AL) method (Arnold and Porter, 2015; Atamna et al., 2016) adapts the augmented Lagrangian

L (x) = f (x) + \sum_{j = 1}^{m} \{\begin{matrix} γ_{j} g_{j} (x) + \frac{ω_{j}}{2} g_{j}^{2} (x) & if γ_{j} + ω_{j} g_{j} (x) ⩽ 0, \\ - \frac{γ_{j}^{2}}{2 ω_{j}} & otherwise, \end{matrix}

where

γ_{j} \in R

is a Lagrange factor,

ω_{j} \in R^{+}

is a penalty coefficient, and both are adapted during the optimization process. AL was initially proposed for (⁠

1 + 1

⁠)-ES (Arnold and Porter, 2015), and was extended to the CMA-ES in a single constraint case (Atamna et al., 2016), where the median success rule was applied for the step-size adaptation. The AL is designed for implicit constraints, in which the implicit constraint means that the evaluation is expensive or is computed at the same time as the objective function. It requires the objective function to be defined in the infeasible domain.

2.2.3 Ranking-Based Methods

Ranking-based CHTs aggregate the rankings of the objective function values and constraint function values to create the final rankings of the candidate solutions, instead of aggregating the function values. The stochastic ranking technique (Runarsson and Yao, 2000) attempts to balance the objective and constraint functions by sorting the candidate solutions according to the objective function with a probability $P_{f}$ ⁠, and the sum of the constraint function values with a probability $1 - P_{f}$ ⁠. The multiple constraint ranking (MCR) technique (de Paula Garcia et al., 2017) ranks the candidate solutions according to the sum of the rankings of the objective values and rankings of each constraint violation value. Other techniques included in ranking-based CHTs such as the lexicographic ordering and the $ɛ$ -lexicographic ordering have been proposed in GA and DE, and imported to ES (Oyman et al., 1999; Hellwig and Beyer, 2018). The advantage of these approaches over penalty function-based approaches is that they are invariant to strictly increasing transformation of the objective function and constraint functions, which we later refer to as element-wise increasing transformation. Therefore, practitioners do not need to tune the balance between the objective function values and constraint violation values manually. As with the AL, however, this method requires the objective function to be defined in the infeasible domain.

2.2.4 Active Constraint Handling

There are certain CHTs that modify the covariance matrix adaptation mechanism to shrink the variance actively in the direction of the constraint function gradient so as to decrease the likelihood of infeasible solutions. We refer to such methods as active constraint handling (ACH) in this article. Reference Arnold and Hansen (2012) proposed the use of the active covariance matrix update in the (1+1)-CMA-ES. Similar concepts have been employed in other variants of the CMA-ES, such as the (⁠ $μ, λ$ ⁠)-CMA-ES in Chocat et al. (2015), MA-ES in Spettel and Beyer (2019), and xCMA-ES in Krause and Glasmachers (2015). As ACH techniques use binary information, whether or not the constraint is violated, they are invariant to monotone transformations of the objective and constraint violations. Such methods can be applied to problems in which the constraints return the outcome that the given solution is either feasible or not. However, this approach tends to be inefficient compared to other CHTs on quantifiable constraints, as it does not utilize the amount of constraint violations.

2.2.5 Other Explicit CHTs

The linear constraint covariance matrix self-adaptation ES (lcCMSA-ES) (Spettel et al., 2019) handles explicit and linear constraints in a variant of CMA-ES, known as CMSA-ES. It basically samples only feasible solutions, and updates the distribution parameters by using the feasible solutions. Active-set ES (Arnold, 2016, 2017) is also designed for explicit, but not necessarily linear, constraints. This approach is the most relevant one from the perspective of the assumptions on the constraint problem. It is a (1+1)-ES based approach that applies a repair operator in the Lamarckian manner; that is, the repair solution is used as the candidate solution, and not only to compute the loss function value. Unfortunately, it is not possible to directly extend it to the state-of-the-art variant of the ES, namely (⁠ $μ$ ⁠, $λ$ ⁠)-CMA-ES.

2.3 Formal Classification of CHTs

We summarize the formal classification of the abovementioned CHTs in Table 2. The taxonomy proposed in Le Digabel and Wild (2015) classifies constraints as follows: Quantifiable/Nonquantifiable, Relaxable/Unrelaxable, Simulation-based/A priori, and Known/Hidden. A Quantifiable constraint is a constraint for which the amount of feasibility and/or violation can be quantified, while a Nonquantifiable constraint returns a binary output indicating whether or not the solution satisfies the constraint. A Relaxable constraint is a constraint that does not need to be satisfied to compute the objective function, while an Unrelaxable constraint is a prerequisite for executing the simulation for the objective function computation. An A priori constraint is a constraint for which the feasibility can be confirmed without running a simulation; that is, this constraint can be formulated using optimization variables such as $g (x) = \sum_{i = 1}^{n} {[x]}_{i} ⩽ 1$ ⁠, while a Simulation-based constraint is only computed through a computationally expensive simulation. A Known constraint is a constraint that is explicitly provided in the problem formulation; for example, $min f (x)$ s.t. $g (x) ⩽ 0$ ⁠. Constrained problems can be expressed by combining these initial letters as an acronym, such as QRSK. Refer to Le Digabel and Wild (2015) for further descriptions of each type of constraint and example situations. Our assumption on the constraints is QUAK in this terminology.

Table 2:

Classification of CHTs, where bnds/lc/nlc means that the CHT can handle bound constraints, linear constraints, or nonlinear constraints, respectively.

CHT (bnds/lc/nlc)	Taxonomy	Invariance
ARCH [proposed CHT, described in Section 4] (nlc)	QUAK	increasing / affine
AP-BCH (Hansen et al., 2009) (bnds)	QUAK	$\times$ / $\times$
lcCMSA-ES (Spettel et al., 2019) (lc)	QUAK	increasing / $\times$
Active-set ES (Arnold, 2017) (nlc)	QUAK	increasing / $\times$
AL (Atamna et al., 2016) (nlc)	QRSK	$\times$ / affine
Stochastic ranking (Runarsson and Yao, 2000) (nlc)	QRSK	$\times$ / affine
MCR (de Paula Garcia et al., 2017) (nlc)	QRSK	increasing / affine
(1+1)-CMA-ES with ACH (Arnold and Hansen, 2012) (nlc)	NUSK	increasing / affine
(⁠ $μ, λ$ ⁠)-CMA-ES with ACH (Chocat et al., 2015) (nlc)	NRSK	increasing / affine
xCMA-ES with ACH (Krause and Glasmachers, 2015) (nlc)	NUSK	increasing / affine
MA-ES with ACH (Spettel and Beyer, 2019) (nlc)	NUSK	increasing / affine
Resampling technique (nlc)	NUSH	increasing / affine

CHT (bnds/lc/nlc)	Taxonomy	Invariance
ARCH [proposed CHT, described in Section 4] (nlc)	QUAK	increasing / affine
AP-BCH (Hansen et al., 2009) (bnds)	QUAK	$\times$ / $\times$
lcCMSA-ES (Spettel et al., 2019) (lc)	QUAK	increasing / $\times$
Active-set ES (Arnold, 2017) (nlc)	QUAK	increasing / $\times$
AL (Atamna et al., 2016) (nlc)	QRSK	$\times$ / affine
Stochastic ranking (Runarsson and Yao, 2000) (nlc)	QRSK	$\times$ / affine
MCR (de Paula Garcia et al., 2017) (nlc)	QRSK	increasing / affine
(1+1)-CMA-ES with ACH (Arnold and Hansen, 2012) (nlc)	NUSK	increasing / affine
(⁠ $μ, λ$ ⁠)-CMA-ES with ACH (Chocat et al., 2015) (nlc)	NRSK	increasing / affine
xCMA-ES with ACH (Krause and Glasmachers, 2015) (nlc)	NUSK	increasing / affine
MA-ES with ACH (Spettel and Beyer, 2019) (nlc)	NUSK	increasing / affine
Resampling technique (nlc)	NUSH	increasing / affine

CHTs assuming Unrelaxable constraints can be applied to QUAK constraints. For example, the resampling technique can be applied to QUAK constraints. However, CHTs assuming weaker conditions on constraints utilize less information on the constraints than that available to the optimization approaches. Therefore, we expect that CHTs assuming QUAK constraints are more efficient for solving QUAK constrained optimization problems. Only two approaches in Table 2, including the proposed approach, are designed for QUAK nonlinear constraints. As it is not clear how active-set ES is extended to a variant of the CMA-ES, the proposed approach is the only approach applied to the CMA-ES. However, if the constraints are nonlinear but Relaxable, many of the CHTs listed in Table 2 can be applied.

3 Desired Invariance Properties for Constrained Optimization

We describe two invariance properties that are desirable for a CHT that is designed for variants of the CMA-ES.

3.1 Element-Wise Increasing Transformation of Functions

An strictly increasing transformation $h : R \to R$ is a function satisfying $h (t) < h (s)$ if $t < s$ ⁠. Invariance to an increasing transformation of the objective function in an unconstrained optimization scenario refers to the property whereby the algorithm does not change the behavior (which is possibly characterized by the sequence of the solutions generated by the algorithm) when solving $f$ and its composite $h \circ f$ ⁠. Algorithms that are invariant to any increasing transformation can solve, for example, a non-convex discontinuous function $h \circ f$ as easily as a convex continuous functions $f$ ⁠. The importance of this invariance property is extensively recognized in engineering optimization: if the search algorithm is not invariant, the objective function needs to be tuned for the search algorithm to perform effectively, which is time consuming. Numerous evolutionary algorithms, including the CMA-ES, are invariant to any increasing transformation of the objective function, because they use only the objective function value rankings of the candidate solutions.

In constrained optimization, we consider invariance to an element-wise increasing transformation $H = (h_{0}, \dots, h_{m}) : R^{m + 1} \to R^{m + 1}$ of the objective and constraint functions $F = (f, g_{1}, \dots, g_{m})$ ⁠, where $h_{j} : R \to R$ (⁠ $j = 0, \dots, m$ ⁠) is an increasing transformation and $h_{j}$ for $j = 1, \dots, m$ satisfies $h_{j} (0) = 0$ ⁠. Invariance to an element-wise increasing transformation of the objective and constraint functions refers to the property whereby the algorithm does not change the behavior when solving a constrained problem $F = (f, g_{1}, \dots, g_{m}) = (h_{0} \circ f, h_{1} \circ g_{1}, \dots, h_{m} \circ g_{m}) = H \circ F$ ⁠. The original constrained optimization problem $F$ and its transformation $H \circ F$ models the same optimization problem in the sense that they define the same feasible domain $X$ and the same total order on $X$ regarding the objective function value.^³

In real-world applications, the ranges of the objective function values and constraint violations are often quite different. Algorithms without this invariance will suffer from this difference and place implicit priority on an objective or certain constraints depending on their values. Therefore, practitioners may determine a reasonable transformation $H$ ⁠. As described above for unconstrained optimization cases, this can be time consuming, and requires domain knowledge of the problem and deep insight into the optimization algorithm.

Although this invariance property is a straightforward extension of invariance to the increasing transformation of the objective function and it is seemingly important, it is not exhibited by the frequently used penalty function-based techniques that take the sum of the objective and constraint function values as loss values. For example, it is clear that the AL does not exhibit this invariance: it is not invariant to increasing the transformation of $f$ and $g_{j}$ ⁠. Although the AP-BCH is invariant to any increasing transformation of $g_{j}$ ⁠, it is not invariant to an increasing transformation of $f$ in general, as the quadratic penalty term is added to the objective function value directly. It has been demonstrated in Sakamoto and Akimoto (2017) that the AP-BCH deteriorates when the objective function is, for example, an exponential function, where the objective function value is more sensitive than the quadratic penalty term.

3.2 Affine Transformation of Search Space Coordinates

To formulate the invariance properties for constrained optimization problems, we consider the fact that our minimization problem is defined on an

n

-dimensional inner product space

(V, 〈 \cdot, \cdot 〉)

on the real field

R

⁠:

\underset{p \in V}{\arg \min} f^{V} (p) s . t . g_{j}^{V} (p) ⩽ 0, \forall j = 1, \dots, m,

(4)

where

f^{V} : V \to R

and

g_{j}^{V} : V \to R

are the objective and constraint functions, respectively. The constrained problem () is considered a realization of () under an orthonormal basis

{e_{i} \in V}_{i = 1}^{n}

and a bias vector

e_{0} \in V

⁠, where

〈 e_{i}, e_{i} 〉 = 1

for

i = 1, \dots, n

and

〈 e_{i}, e_{j} 〉 = 0

for any

j \neq i

(⁠

i, j ⩾ 1

⁠), and

x \in R^{n}

corresponds to

p = e_{0} + \sum_{i = 1}^{n} {[x]}_{i} e_{i} \in V

⁠. Let

p \mapsto x

be denoted by

ψ_{e} : V \to R^{n}

⁠. If the optimization algorithm is defined on the inner product space, its behavior (which is possibly characterized by the sequence of the solutions generated by the algorithm) is identical on any basis

{v_{i} \in V}_{i = 1}^{n}

with a bias vector

v_{0} \in V

⁠, where the bases are not necessarily orthonormal to one another. Let

p = v_{0} + \sum_{i = 1}^{n} {[\tilde{x}]}_{i} v_{i} \mapsto \tilde{x}

be denoted by

ψ_{v} : V \to R^{n}

⁠. The map

A : ψ_{e} \circ ψ_{v}^{- 1}

from a coordinate system

ψ_{v}

to a coordinate system

ψ_{e}

is an affine transformation. The invariance to an affine transformation of the search space coordinates refers to the property whereby the algorithm behaves the same on the original coordinate system and its affine transformed coordinate system.

The objective function

f^{V} (p) = \frac{1}{2} 〈 p, p 〉

is expressed as

f^{e} (x) = f^{V} (ψ_{e}^{- 1} (x)) = \sum_{i = 1}^{n} {[x - ψ_{e} (e_{0})]}_{i}^{2}

in the system

ψ_{e}

⁠, while it is expressed as

f^{v} (\tilde{x}) = f^{e} (A (\tilde{x})) = \sum_{i = 1}^{n} {[A \tilde{x} - ψ_{e} (v_{0})]}_{i}^{2}

in the system

ψ_{v}

⁠. The former is the sphere function (well conditioned and separable), while the latter is a convex quadratic function that is ill conditioned if

A

also is, and is nonseparable if

A

is not diagonal. If the algorithm is invariant to any affine transformation of the search space coordinates, the algorithm solves the ill-conditioned and nonseparable function

f^{v}

as efficiently as it solves the well-conditioned and separable

f^{e}

under the corresponding transformation of the initial search distribution. Not all variants of the CMA-ES, including that presented in Section 2.1, are (proven to be) invariant to any affine transformation of the search space coordinate,^⁴ although we empirically observe statistically invariant behaviors under arbitrary affine transformations (see Figure 3). Numerous other evolutionary computation approaches are not invariant to these transformations, and we empirically observe rather different performances depending on the affine transformation properties.

Figure 1:

For example, J(x1)={1,2,3} and A(x1)=∅, then x1 is repaired by using Eq. (6).

View large Download slide

For example, $J (x_{1}) = {1, 2, 3}$ and $A (x_{1}) = \emptyset$ ⁠, then $x_{1}$ is repaired by using Eq. (6).

Figure 2:

View large Download slide

Same problem with different coordinate systems.

Figure 3:

View large Download slide

Median (line) and 25% to 75%-ile range (band) over 100 trials.

In real-world applications, the change in the coordinate system corresponds to the change in the features describing the object to be optimized, or the change in the unit in each feature. In a constrained optimization scenario, the constraint function $g_{j}$ and the objective function $f$ are transformed by the same transformation $A$ ⁠, resulting in $g_{j} \circ A$ and $f \circ A$ in the transformed coordinate system, respectively. The linear constraints are again linear in the transformed coordinate system. However, a box constraint will no longer be a box constraint, but rather a set of linear constraints (forming an $n$ -parallelotope shape).

Suppose that the underlying unconstrained optimization algorithm is invariant to any affine transformation of the search space coordinates. As summarized in Table 2, CHTs that only touch the loss function values, such as the resampling technique, ACH, and AL, do not disturb the invariance property of the underlying algorithm, resulting in invariance to the affine transformation of the search space under the constraints (see Atamna et al. (2020) for proof of the affine invariance of AL). However, CHTs for unrelaxable constraints that require a repair operator often lose the invariance property. For example, the AP-BCH does not exhibit this invariance property, as the repair operator used in this approach exploits the fact that the feasible domain is composed of an interval in each coordinate. Repair operators that use the inner product or the distance in Euclidean space are generally affected by the affine transformation of the coordinate system.

4 ARCH

We propose an explicit constraint-handling method based on adaptive ranking, known as ARCH. ARCH virtually transforms the constrained optimization problem into an unconstrained problem by defining an adaptive loss function $L$ ⁠. The loss $L$ is defined by the weighted sum of the rankings of the objective and constraint violations, denoted as $L = R_{T}$ below. The proposed algorithm is designed to exhibit the invariance properties listed in Section 3, and it does not require the objective function to be well defined in the infeasible domain.

4.1 Repair Operator

To make ARCH applicable to a problem in which the objective function values are not defined in the infeasible domain, we employ a repair operator in the Darwinian manner. Given a candidate solution $x$ ⁠, ARCH determines a repaired solution $\tilde{x}$ ⁠. The repaired solution $\tilde{x}$ is only used for the objective function value computation.

The repair operator

Repair (\cdot)

is defined as follows: Let

S = {x \in R^{n} ∣ g_{j} (x) ⩽ 0, \forall j \in [[1, m]]}

be the feasible domain. Given a solution

x

⁠, let

J (x) = {j \in [[1, m]] ∣ g_{j} (x) > 0}

be the set of indices of the violated constraints. Let

A (x) = {y \in R^{n} ∣ g_{j} (y) = 0, \forall j \in J (x)}

be the intersection of the violated constraint boundaries (see Figure 1). We introduce the repair operator defined as follows:

\begin{matrix} \tilde{x} = Repair (x) = {\arg \min}_{y} {∥ x - y ∥}_{Σ^{- 1}}^{2} s . t . y \in \{\begin{matrix} A (x) \cap S & if A (x) \cap S \neq \emptyset \\ S & otherwise \end{matrix}, \begin{matrix} (5) \\ (6) \end{matrix} \end{matrix}

where

{∥ x - y ∥}_{Σ^{- 1}}^{2} = {(x - y)}^{T} Σ^{- 1} (x - y)

is the Mahalanobis distance between

x

and

y

under the inverse matrix

Σ^{- 1} = {(σ^{2} C)}^{- 1}

⁠. This is a constrained minimization of a quadratic function, which is solved by standard numerical optimization routines. If the constraints are all linear, the problem is reduced to a quadratic programming problem. Equation (6) (see Figure 1) is the nearest feasible solution to

x

⁠, while Eq. () is the feasible solution nearest to

x

under the constraint that the active constraints at

x

remain active at

\tilde{x}

⁠. Eq. () is preferred over Eq. (6).

The condition that $A (x) \cap S \neq \emptyset$ is indeed simplified to $A (x) \neq \emptyset$ when the constraints are all linear and not redundant. Here, a redundant constraint $g_{j}$ is defined as a constraint such that the boundary of $g_{j}$ is never in contact with the feasible domain; that is, ${x \in R^{n} ∣ g_{j} (x) = 0} \notin S$ ⁠. Note that determining whether $A (x) \neq \emptyset$ is as easy as checking whether a system of linear equations has a solution. In explicit linear constraint situations, such a redundant constraint can be removed in advance.

In this approach, the Mahalanobis distance is employed rather than the Euclidean distance ${∥ x - y ∥}^{2} = {(x - y)}^{T} (x - y)$ ⁠, in order to achieve affine invariance of the search space coordinate system. This point is discussed formally in Section 5.

4.2 Total Ranking

The total ranking

R_{T} (x_{k})

for

k \in {1, \dots, λ}

at each iteration

t \in {0, 1, \dots}

is a weighted sum of the rankings of the candidate solutions based on the objective function values,

R_{f} (x_{k})

⁠, and based on the Mahalanobis distance to the feasible domain,

R_{g} (x_{k})

⁠, namely:

R_{T} (x_{k}) = R_{f} (x_{k}) + α^{(t + 1)} R_{g} (x_{k}),

(7)

where

α^{(t + 1)}

is known as the ranking coefficient that controls the balance between the objective and constraints. The rankings

R_{f} (x_{k})

and

R_{g} (x_{k})

are defined as follows:

The

f

-ranking,

R_{f} (x_{k})

⁠, is the number of better candidate solutions in terms of

f

⁠, plus the number of tie candidate solutions divided by 2, namely

R_{f} (x_{k}) = \sum_{l = 1}^{λ} 1_{{f ({\tilde{x}}_{l}) < f ({\tilde{x}}_{k})}} + \frac{1}{2} \sum_{l = 1}^{λ} 1_{{f ({\tilde{x}}_{l}) = f ({\tilde{x}}_{k})}},

(8)

where

{\tilde{x}}_{i} = Repair (x_{i})

⁠. Note that the sum of

R_{f} (x_{k})

for

k \in {1, \dots, λ}

is

λ^{2} / 2

⁠. Moreover, in an unconstrained optimization scenario, the probability of sampling tie solutions is often zero, while the second term above may be nonzero with a nonzero probability in our situation owing to the repair operation.

The

g

-ranking,

R_{g} (x_{k})

⁠, is analogously defined by simply replacing

f \circ Repair

with

g_{Σ}

⁠:

R_{g} (x_{k}) = \sum_{l = 1}^{λ} 1_{{g_{Σ} (x_{l}) < g_{Σ} (x_{k})}} + \frac{1}{2} \sum_{l = 1}^{λ} 1_{{g_{Σ} (x_{l}) = g_{Σ} (x_{k})}},

(9)

where

g_{Σ} (x) = {∥ x - \tilde{x} ∥}_{Σ^{- 1}}^{2}

and

Σ = {(σ^{(t)})}^{2} C^{(t)}

⁠.

4.3 Adaptation of Ranking Coefficient

The ranking coefficient $α$ controls the balance between $R_{f}$ and $R_{g}$ ⁠. If $α$ is too large, the search distribution is biased toward the feasible domain. If $α$ is too small, the search distribution is biased toward the infeasible domain. Therefore, the adaptation of $α$ significantly influences the search performance.

The adaptation of

α

is based on a theoretical study of the weighted recombination ES on a spherical function. It was reported in Akimoto et al. (2020) that, in the optimal situation in terms of quality gain (Beyer, 2001)

\frac{∥ m_{sph} - x_{sph}^{*} ∥}{n \cdot σ} = \frac{1}{{\bar{σ}}^{*}} \approx \frac{1}{\hat{σ}}, \hat{σ} = \frac{c \cdot n \cdot μ_{ef f}}{n - 1 + c^{2} \cdot μ_{ef f}},

(10)

where

m_{sph}

and

x_{sph}^{*}

are the mean vector and optimum on a spherical function, respectively;

{\bar{σ}}^{*}

is the optimal normalized step size, which is approximated by

\hat{σ}

⁠;

c = - \sum_{i = 1}^{μ} w_{i} E [N_{i : λ}]

is the weighted average of the expected value of the normal order statistics from

λ

samples and is usually in

O (1)

⁠; and

μ_{ef f} = {(\sum_{i = 1}^{μ} w_{i}^{2})}^{- 1}

is the variance effective selection mass.

Ideally, the CMA-ES with the proposed constraint handling should treat a constrained sphere function as though it is an unconstrained sphere function. That is, we wish (10) to hold even for a constrained sphere problem.

Assuming that the optimum of a constrained

n

-dimensional sphere problem having

n

constraints is located on the

n

boundaries, we estimate the left-hand side of () using the Mahalanobis distance between

m

violating

n

constraints and

\tilde{m} = Repair (m)

⁠:

\frac{{∥ m_{sph} - x_{sph}^{*} ∥}^{2}}{σ^{2}} \approx {∥ m - \tilde{m} ∥}_{Σ^{- 1}}^{2} .

(11)

We define

d_{m}^{(t + 1)}

using the parameter at iteration

t

(⁠

m = m^{(t)}, Σ^{- 1} = {({(σ^{(t)})}^{2} C^{(t)})}^{- 1}

⁠), as follows:

d_{m}^{(t + 1)} = \frac{{∥ m - \tilde{m} ∥}_{Σ^{- 1}}^{2} \cdot {\hat{σ}}^{2}}{n (n / 2 + c_{act})} exp (\frac{min (0, λ_{def} - λ)}{λ}),

(12)

where

c_{act} = | {j ∣ g_{j} (\tilde{m}) = 0, \forall j} |

is the number of active constraints at the repaired mean

\tilde{m}

⁠. From the perspective of () and (), we wish to maintain

d_{m} \approx 1

if the denominator is

n^{2}

⁠. However, in (), we replace

n^{2}

with

n (n / 2 + c_{act})

⁠, leading to a variation in the denominator in

[n^{2} / 2, 3 n^{2} / 2]

⁠. The motivation is to incorporate the fact that the numerator is expected to be smaller if

c_{act}

is smaller, as the projection is performed only on this

c_{act}

-dimensional subspace. Moreover, we introduce the exponential term to prevent the search distribution from being unnecessarily biased toward the boundary if the population size

λ

is larger than the default

λ_{def}

⁠. If

λ \in Ω (n)

⁠,

\hat{σ}

will be in

O (n)

⁠, while

\hat{σ} \in O (μ_{w})

if

λ ≪ n

⁠. Then, to maintain

d_{m} \approx 1

⁠,

m

needs to be closer to

\tilde{m}

⁠, as

λ

is greater. That is, a larger

λ

results in the search distribution being more biased toward the feasible domain, and less efficient search performance. To mitigate this problem, we allow

m

to be up to

e

times away from

\tilde{m}

by introducing the exponential term when

λ

is greater than

λ_{def}

⁠. We adapt

α

so that

d_{m}

will remain approximately 1, as follows:

Ranking Coefficient Adaptation Initialize

α^{(0)} = 1

⁠,

d_{m}^{(0)} = 0

⁠. At iteration

t

⁠,

α^{(t)}

is updated as

α^{(t + 1)} = α^{(t)} \cdot (\frac{sgn (d_{m}^{(t + 1)} - 1)}{n})

(13)

only if

sgn (d_{m}^{(t + 1)} - 1) = sgn (d_{m}^{(t + 1)} - d_{m}^{(t)})

or

d_{m}^{(t + 1)} = 0

⁠, the latter of which is necessary to decrease

α

when the mean vector remains in the feasible domain; that is,

d_{m}^{(t)} = d_{m}^{(t + 1)} = 0

⁠. Following the update,

α^{(t + 1)}

is clipped to

[1 / λ, λ]

⁠, because

α < 1 / λ

and

α > λ

result in

R_{g}

and

R_{f}

being ignored, respectively.^⁵

5 Invariance Properties of ARCH

We prove that ARCH is invariant to the problem transformations described in Section 3. ARCH receives the candidate solutions

{(x_{k})}_{k = 1}^{λ}

⁠, and distribution parameters

m

and

Σ = σ^{2} C

⁠, from the CMA-ES, and returns the total rankings

{(R_{T} (x_{k}))}_{k = 1}^{λ}

of the candidate solutions to the CMA-ES. Meanwhile, it maintains the penalty coefficient

α

and retains

d_{m}

for the next update. Therefore, the functionality of ARCH can be formulated as

({(R_{T} (x_{k}))}_{k = 1}^{λ}, α^{(t + 1)}, d_{m}^{(t + 1)}) = A R C H ({(x_{k})}_{k = 1}^{λ}, m, Σ^{- 1}, α^{(t)}, d_{m}^{(t)}, F),

(14)

where

F : x \mapsto (f (x), g_{1} (x), \dots, g_{m} (x))

is the vector-valued function consisting of the objective and constraints. In the following, we prove that the outputs of () remain unchanged by the abovementioned transformations.

5.1 ARCH on Inner Product Space

We begin by defining ARCH on an inner product space $(V, 〈 \cdot, \cdot 〉)$ ⁠.

A random vector $X_{V}$ on the inner product space $(V, 〈 \cdot, \cdot 〉)$ has a normal distribution $N (μ_{V}, Σ_{V})$ if $〈 v, X 〉$ has the normal distribution $N (〈 v, μ_{V} 〉, 〈 v, Σ_{V} (v 〉))$ for any $v \in V$ ⁠, where $μ_{V} \in V$ is the mean vector and $Σ_{V} : V \to V$ is the covariance (Eaton, 2007). Let $F^{V} = (f^{V}, g_{1}^{V}, \dots, g_{m}^{V})$ be the objective and constraint functions defined on $V$ as (4). ARCH on the inner product space receives these distribution parameters $μ_{V}$ and $Σ_{V}$ ⁠, candidate solutions $p_{1}, \dots, p_{λ}$ drawn from the normal distribution, and the objective and constraint functions $F^{V}$ ⁠.

Let

p

and

q

be the vectors in

V

⁠. The Mahalanobis distance between

p

and

q

under the covariance

Σ_{V}

is defined as

{∥ p - q ∥}_{Σ_{V}^{- 1}}^{2} : = 〈 p - q, Σ_{V}^{- 1} (p - q) 〉 .

Let

S^{V} = {p \in V ∣ g_{j}^{V} (p) ⩽ 0, \forall j \in [[1, m]]}

be the feasible domain. Given a solution

p \in V

⁠, let

J^{V} (p) = {j \in [[1, m]] ∣ g_{j}^{V} (p) > 0}

be the set of indices of the unsatisfied constraints, and let

A^{V} (p) = {q \in V ∣ g_{j}^{V} (q) = 0, \forall j \in J^{V} (p)}

be the intersection of the violated constraint boundaries. The repair operation is defined as follows:

{Repair}^{V} (p) = \underset{q \in V}{arg min} {∥ p - q ∥}_{Σ_{V}^{- 1}}^{2} s . t . q \in \{\begin{matrix} A^{V} (p) \cap S^{V} & if A^{V} (p) \cap S^{V} \neq \emptyset \\ S^{V} & otherwise . \end{matrix}

(15)

Let

d_{m}

and

α

be computed as in () and (), where

m

⁠,

\tilde{m} = Repair (m)

⁠,

Σ^{- 1}

⁠, and

c_{act}

are replaced with

μ_{V}

⁠,

{\tilde{μ}}_{V} = {Repair}^{V} (μ_{V})

⁠,

Σ_{V}^{- 1}

⁠, and

c_{act}^{V} = | {j ∣ g_{j}^{V} ({\tilde{μ}}_{V}) = 0, \forall j} |

⁠, respectively. Analogously to Eq. () to (9), we compute

R_{T}^{V} (p_{k}) = R_{f}^{V} (p_{k}) + α R_{g}^{V} (p_{k}),

(16)

where

R_{f}^{V} (p_{k}) = \sum_{l = 1}^{λ} 1_{{f^{V} ({\tilde{p}}_{l}) < f^{V} ({\tilde{p}}_{k})}} + \frac{1}{2} \sum_{l = 1}^{λ} 1_{{f^{V} ({\tilde{p}}_{l}) = f^{V} ({\tilde{p}}_{k})}}

(17)

R_{g}^{V} (p_{k}) = \sum_{l = 1}^{λ} 1_{{g_{Σ}^{V} (p_{l}) < g_{Σ}^{V} (p_{k})}} + \frac{1}{2} \sum_{l = 1}^{λ} 1_{{g_{Σ}^{V} (p_{l}) = g_{Σ}^{V} (p_{k})}},

(18)

in which

g_{Σ}^{V} (p) = {∥ p - \tilde{p ∥}}_{Σ_{V}^{- 1}}^{2}

and

\tilde{p} = {Repair}^{V} (p)

⁠.

The operation of ARCH on the inner product space

(V, 〈 \cdot, \cdot 〉)

is then expressed as

({(R_{T}^{V} (p_{k}))}_{k = 1}^{λ}, α^{(t + 1)}, d_{m}^{(t + 1)}) = A R C H^{V} ({(p_{k})}_{k = 1}^{λ}, μ_{V}, Σ_{V}^{- 1}, α^{(t)}, d_{m}^{(t)}, F^{V}) .

(19)

Note that

A R C H^{V}

is defined without relying on any coordinate system. This implies that

A R C H^{V}

is invariant to any coordinate system transformation.

5.2 Invariance to Affine Transformation of Search Space Coordinates

Firstly, we demonstrate that ARCH is invariant to an arbitrary affine transformation of the search space coordinate system. For this purpose, it is sufficient to prove that the operation of ARCH (14) on any given coordinate system $x = ψ (p)$ is equivalent to the operation of ARCH (19) on the inner product space $(V, 〈 \cdot, \cdot 〉)$ ⁠.

Note that the normal distribution $N (μ_{V}, Σ_{V})$ on $V$ corresponds to the normal distribution $N (μ_{ψ}, Σ_{ψ})$ on $R^{n}$ in the coordinate system $ψ : V \to R^{n}$ ⁠, where $μ_{ψ} = ψ (μ_{V})$ and ${[Σ_{ψ}]}_{i, j} = 〈 e_{i}, Σ_{V} (e_{j}) 〉$ ⁠. The objective and constraint function $F : R^{n} \to R^{m + 1}$ corresponding to $F^{V} : V \to R^{m + 1}$ is expressed as $F = F^{V} \circ ψ^{- 1}$ ⁠.

Theorem 1.

Let

ψ : V \to R^{n}

be an arbitrary coordinate system. Let

μ_{V} \in V

be an arbitrary vector and

Σ_{V} : V \to V

be an arbitrary positive definite symmetric linear transformation. Let

μ_{ψ} = ψ (μ_{V})

and

{[Σ_{ψ}]}_{i, j} = 〈 e_{i}, Σ_{V} (e_{j}) 〉

and

F = F^{V} \circ ψ^{- 1}

⁠, assume that the feasible domain

S

is convex. Then, for any

p_{k} \in V

for

k = 1, \dots, λ

⁠,

α > 0

and

d_{m} > 0

⁠,

A R C H ({(ψ (p_{k}))}_{k = 1}^{λ}, μ_{ψ}, Σ_{ψ}^{- 1}, α, d_{m}, F) = A R C H^{V} ({(p_{k})}_{k = 1}^{λ}, μ_{V}, Σ_{V}^{- 1}, α, d_{m}, F^{V}) .

5.3 Invariance to Element-Wise Increasing Transformation

Next, we demonstrate that ARCH is invariant to an arbitrary element-wise increasing transformation. That is, the ranking and internal parameter updates are not affected by the transformation.

Theorem 2.

Let

H = (h_{0}, h_{1}, \dots, h_{m})

be an arbitrary element-wise increasing transformation; assume that the feasible domain

S

is convex. For any

x_{k} \in R^{n}

for

k = 1, \dots, λ

⁠,

α > 0

⁠, and

d_{m} > 0

⁠,

A R C H ({(x_{k})}_{k = 1}^{λ}, m, Σ^{- 1}, α, d_{m}, F) = A R C H ({(x_{k})}_{k = 1}^{λ}, m, Σ^{- 1}, α, d_{m}, H \circ F) .

5.4 Invariance Properties of Entire Algorithm

Finally, we discuss the invariance properties of the search algorithm including ARCH.

Suppose that the underlying unconstrained optimization algorithm is defined on an inner product space, and the following steps are repeated:

$p_{1}, \dots, p_{λ} = s a m p l e (μ_{V}, Σ_{V}, θ_{V})$ ⁠;
$R_{f} (p_{1}), \dots, R_{f} (p_{λ}) = e v a l u a t e ({(p_{k})}_{k = 1}^{λ}, f^{V})$ ⁠; and
$μ_{V}, Σ_{V}, θ_{V} = u p d a t e ({(p_{k}, R_{f} (p_{k}))}_{k = 1}^{λ}, μ_{V}, Σ_{V}, θ_{V})$ ⁠,

where $θ_{V}$ contains all of the information used in the algorithm. By nature, it is invariant to any strictly increasing transformation of the objective function $f^{V}$ ⁠, as the outputs of $e v a l u a t e$ are the rankings of the objective function values, and it is invariant to any strictly increasing transformation $h$ ⁠. Moreover, as it is defined independently of a coordinate system, its implementation under a given coordinate system $ψ$ produces $ψ (p_{1}), \dots, ψ (p_{λ}), μ_{ψ}, Σ_{ψ}, θ_{ψ}$ at any iteration $t$ ⁠, where $p_{1}, \dots, p_{λ}$ are the candidate solutions mentioned above, while $μ_{ψ}$ ⁠, $Σ_{ψ}$ ⁠, and $θ_{ψ}$ are the expressions of $μ_{V}$ ⁠, $Σ_{V}$ ⁠, and $θ_{V}$ ⁠, respectively, in the coordinate system $ψ$ ⁠.

When solving an explicitly constrained minimization problem $F^{V}$ ⁠, we simply replace the evaluation step as follows:

$p_{1}, \dots, p_{λ} = s a m p l e (μ_{V}, Σ_{V}, θ_{V})$ ⁠;
$R_{T} (p_{1}), \dots, R_{T} (p_{λ}), α, d_{m} = A R C H^{V} ({(p_{k})}_{k = 1}^{λ}, μ_{V}, Σ_{V}^{- 1}, α, d_{m}, F^{V})$ ⁠; and
$μ_{V}, Σ_{V}, θ_{V} = u p d a t e ({(p_{k}, R_{T} (p_{k}))}_{k = 1}^{λ}, μ_{V}, Σ_{V}, θ_{V})$ ⁠.

As all of the operations are defined independently from a coordinate system, the entire algorithm is invariant to any affine coordinate system transformation. Moreover, because ${ARCH}^{V}$ is invariant to any element-wise increasing transformation $H$ ⁠, and the underlying algorithm only relies on the rankings of the candidate solutions, all of the operations are invariant to $H$ ⁠.

ARCH can be combined with a search algorithm that is not generalized to an inner product space. For example, our baseline CMA-ES, which is described in Section 2.1, does not generalize to an inner product space (as $p_{σ}$ is coordinate dependent), although we empirically observe quite uniform behaviors under different coordinate systems, as we observe in Section 6. ARCH does not disturb the uniform behavior, as ARCH itself is invariant to any affine coordinate transformation.

6 Experiments on Linearly Constrained Quadratic Problems

Our first numerical experiments are aimed at demonstrating the effects of the invariance to affine coordinate transformations of ARCH. We develop a set of linearly constrained test problems. We compare ARCH with other CHTs for the CMA-ES, demonstrating the manner in which affine coordinate transformations affect the performances of CHTs that are not invariant thereto, while ARCH performs equally effectively under different transformations. Moreover, we compare ARCH with CHTs that are specialized for box constraints. Note that box CHTs are expected to be more efficient for box-constrained optimization problems than other CHTs for general linear constraints. We observe that ARCH is competitive with and sometimes outperforms these, even on box-constrained problems.

6.1 Linearly Constrained Quadratic Minimization Problems

We consider a linearly constrained problem (P0), defined in the $n$ -dimensional inner product space $(V, 〈 \cdot, \cdot 〉)$ on the real field $R$ as (4), where the constraints are defined as $g_{i} (p) = 〈 - v_{i}, p 〉 - {[LB]}_{i}$ and $g_{n + i} (p) = 〈 v_{i}, p 〉 - {[UB]}_{i}$ ⁠, where $v_{i} \in V$ is the normal vector of the $i$ th constraint and these are orthogonal to one another; that is, $〈 v_{i}, v_{j} 〉 = δ_{i, j}$ ⁠, and ${[LB]}_{i}, {[UB]}_{i} \in R$ are the lower and upper bounds, respectively.

We consider the coordinate system

ψ : V \to R^{n}

with

{v_{i}}

as the basis vectors. That is,

ψ : p = \sum_{i = 1}^{n} {[x]}_{i} v_{i} \mapsto x

⁠. On this coordinate system, (P0) can be expressed as the box-constrained minimization problem (P1), as follows:

\underset{x \in R^{n}}{arg min} f_{(P 1)} (x) = f (\sum_{i = 1}^{n} {[x]}_{i} v_{i}) s . t . g_{(P 1)} (x) = A x - b ⪯ 0,

(20)

where

A = {[- I_{n}, I_{n}]}^{T}

and

= [{LB}^{T}, {UB}^{T}]

⁠, and

g_{(P 1)} = (g_{1}, \dots, g_{m})

is a vector form of

m = 2 n

constraint functions. We define an initial mean vector and initial covariance matrix as

m^{(0)}

and

C^{(0)}

⁠, respectively, on this coordinate system.

By taking another basis

{w_{i}}

⁠, the general linearly constrained optimization problem (P2) can be defined as

\underset{y \in R^{n}}{arg min} f_{(P 2)} (y) = f_{(P 1)} (P y) s . t . g_{(P 2)} (y) = A P y - b ⪯ 0 .

(21)

On this coordinate system, the coordinate vector

y \in R^{n}

is transformed as

x = P y

using a basis-transformation matrix

P

⁠, which transforms

{v_{i}}

into

{w_{i}}

⁠. The initial mean vector and covariance matrix are transformed as

P^{- 1} m^{(0)}

and

P^{- 1} C^{(0)} {(P^{- 1})}^{T}

⁠, respectively.

These optimization problems (P1) and (P2) are equivalent to (P0). However, (P1) has a box constraint, while (P2) has a set of linear constraints. Algorithms that are invariant to any affine transformation of the search space must perform equivalently on these problems. We use these problems to assess the invariance properties.

6.2 Settings

For the test problem (P1) in (20), we use three objective functions: sphere (⁠ $f_{sph} (x) = \sum_{i = 1}^{n} {[x]}_{i}^{2}$ ⁠), ellipsoid (⁠ $f_{ell} (x) = \sum_{i = 1}^{n} 10^{6 \frac{i - 1}{n - 1}} {[x]}_{i}^{2}$ ⁠), and rotated ellipsoid (⁠ $f_{rotell} (x) = f_{ell} (Q_{θ} x)$ ⁠), where $Q_{θ} \in R^{n \times n}$ is a block diagonal matrix such that each block is a $2 \times 2$ orthogonal (that is, rotation) matrix and all blocks share the same matrix. We let the counter-clockwise rotation angle of each block be $θ = π / 6$ in this experiment.

Problem (P1) is a box-constrained problem. The lower and upper bounds are set as $LB = {[- 1, 1, \dots, - 1, 1]}^{T}$ and $UB = LB + {[5, \dots, 5]}^{T}$ ⁠, respectively. The optimal solution, which is obtained by the Karush–Kuhn–Tucker (KKT) condition, is located at $x^{*} = {[0, 1, \dots, 0, 1]}^{T}$ for $f_{sph}$ and $f_{ell}$ ⁠, and $x^{*} \approx {[0.37, 1, \dots, 0.37, 1]}^{T}$ for $f_{rotell}$ ⁠. That is, the even-numbered coordinates of the optimum are on the boundary, while the others are not. The number of active constraints at the optimum is $n / 2$ ⁠.

We use the following two matrices for the transformation matrix $P$ in (21): $P_{rot} = Q_{θ^{'}}$ and $P_{illrot} = Q_{θ^{'}}^{T} D Q_{θ^{'}}$ with $θ^{'} = π / 4$ ⁠, where $Q_{θ^{'}}$ is as defined above, with $θ \neq θ^{'}$ ⁠, and $D = diag (1, 10, \dots, 1, 10)$ is a diagonal matrix. The transformation matrix $P_{rot}$ only rotates the search space, while $P_{illrot}$ transforms the rectangular feasible domain into an $n$ -parallelotope shape. Hereunder, problem (P1) is denoted by $Box$ ⁠, the linearly constrained problem (P2) using $P_{rot}$ is denoted by $rotBox$ ⁠, and (P2) using $P_{illrot}$ is denoted by $illrotBox$ ⁠. Their feasible domains and the levelsets of the objective function are illustrated in Figure 2.

For the box-constrained optimization problem (P1), the initial mean vector is $m^{(0)} = \frac{UB + LB}{2} + U {(- 1, 1)}^{n}$ ⁠, and the initial covariance matrix is $C^{(0)} = I_{n}$ ⁠. For the linearly constrained optimization problem (P2), $m^{(0)}$ and $C^{(0)}$ are transformed in the manner described in Section 6.1. The search space dimension is $n \in {20, 50}$ ⁠, the initial step size is $σ^{(0)} = \frac{1}{n} \sum_{i = 1}^{n} \frac{{[UB]}_{i} - {[LB]}_{i}}{4} = 1.25$ ⁠, and the other parameters are set to their default values, as defined in Table 1.

6.3 Results and Discussion

Figure 3 presents the performance of ARCH, the resampling technique with a maximum resampling number of 500, and AP-BCH (Hansen et al., 2009) on

f_{sph}

⁠,

f_{ell}

⁠, and

f_{rotell}

under three different coordinate systems (⁠

Box

⁠,

rotBox

⁠, and

illrotBox

⁠) on

n = 20

and

n = 50

dimensions. As AP-BCH is a box CHT, the results are illustrated only for

Box

⁠. The optimization progress is measured by the Mahalanobis distance between the mean vector and optimal solution

∥ m - x^{*} ∥_{H}^{2} = {(m - x^{*})}^{T} H (m - x^{*})

given the Hessian matrix

H \in R^{n \times n}

of the objective function. Figure 4 presents the results of typical runs of the CMA-ES on unconstrained problems, resampling, AP-BCH, and ARCH on

f_{sph}

and

f_{ell}

under the

Box

constraint on

n = 20

dimensions. For AP-BCH and ARCH, the ratio of the number of constraints satisfied by the mean vector

m

⁠,

r_{feas} = \frac{1}{n} |{i \in [[1, n]] | {[LB]}_{i} ⩽ {[m]}_{i} ⩽ {[UB]}_{i}}|,

(22)

is illustrated for discussion.

Figure 4:

View large Download slide

Typical single runs of CMA-ES (first row) on unconstrained $f_{sph}$ and $f_{ell}$ ⁠, resampling (second row), AP-BCH (third row), and ARCH (fourth row) on $f_{sph}$ and $f_{ell}$ with $Box$ ⁠. The figures indicate the Mahalanobis distance $∥ m - x^{*} ∥_{H}^{2}$ between the mean vector $m$ and optimal solution $x^{*}$ ⁠, given the Hessian matrix $H$ of the objective $f$ ⁠; the step size $σ$ ⁠; the eigenvalues $eig (\sqrt{C})$ of the square root of $C$ ⁠; the coordinates of $m$ ⁠; the ratio $r_{feas}$ of the constraints satisfied by $m$ ⁠; the coefficient $α$ and the parameter $d_{m}$ used for updating $α$ versus the number of iterations.

Firstly, we focus on the results of ARCH and the resampling technique in Figure 3. Although the underlying CMA-ES is not mathematically proven to be invariant to affine transformation of the search space, we observe in Figure 3 that the lines of ARCH for $Box$ ⁠, $rotBox$ ⁠, and $illrotBox$ overlap one another, owing to the invariance of ARCH to the affine transformation of the search space. For resampling as well, we observe the same. This indicates that the performance comparison can be conducted on the most convenient case, namely the $Box$ constraint case. In the following, we focus on analyzing the results on $Box$ ⁠, which can be generalized to general cases.

Next, we investigate the behavior of ARCH. We observe in Figure 4 that the parameter $d_{m}$ is maintained at approximately 1. This is the desired behavior, as we design the adaptation of $α$ to maintain $d_{m} \approx 1$ ⁠. Moreover, we can observe from the results that the behaviors of the CMA-ES on the unconstrained problem and the behaviors of the CMA-ES with ARCH and AP-BCH are similar, in that the covariance matrix tends to be proportional to the inverse of the Hessian matrix of the objective function. However, in Figure 3, we observe the difference between these two algorithms in terms of the adaptation speed of the covariance matrix. On $f_{ell}$ and $f_{rotell}$ ⁠, the CMA-ES with ARCH adapts the covariance matrix significantly faster than that with AP-BCH. By adapting the coefficient $α$ ⁠, ARCH appears to resemble the selection of candidate solutions on an unconstrained problem better than AP-BCH.

The resampling technique exhibits different behavior. The eigenvalues of the covariance matrix are divided into two, with the smaller values corresponding to the axes where the constraints are active at the optimum, and the greater values corresponding to the axes where the constraints are inactive at the optimum. If we do not use any information from the infeasible domain and generate only feasible candidate solutions, this is a reasonable approach to cause the distribution to be narrow in the directions of the active constraints. This concept has indeed been employed in Arnold and Hansen (2012). However, this strategy significantly reduces the speed of the approach of the mean vector to the optimum on the boundary.

7 Experiments on CEC 2006 Constrained Optimization Testbed

The aim of our second numerical experiments is twofold.

Firstly, we demonstrate that ARCH can be applied to nonlinearly as well as linearly constrained problems. To evaluate the efficacy, we compare ARCH with the active-set ES (Arnold, 2017) that can deal with explicit (a priori) and nonlinear constraints, because there is currently only one CHT designed for such a constraint. This is the $(1 + 1)$ -ES-based approach, and the covariance matrix adaptation is not incorporated.

Secondly, we show that, even if the infeasible solutions can be evaluated on the objective, ARCH, which only evaluates the solutions in the feasible domain, exhibits advantages in solving constrained optimization with explicit (a priori) constraints. For this purpose, we compare ARCH with the AL (Atamna et al., 2016) and MCR (de Paula Garcia et al., 2017), which are designed for simulation-based constraints under the relaxed assumption that any infeasible solutions can be evaluated on the objective.

7.1 CEC 2006 testbed

To compare the CHTs, we use the CEC 2006 testbed (Liang et al., 2006), which consists of 24 constrained problems including linear/nonlinear and equality/inequality constraints, where the equality constraints $h_{j} (x) = 0$ are transformed into inequalities, as described below Eq. (1), with $ɛ_{eq} = 10^{- 4}$ ⁠. We reuse the CEC 2006 testbed since this testbed is widely used in existing works, which makes it easy to compare algorithms. All optimization variables are bounded as ${[LB]}_{i} ⩽ {[x]}_{i} ⩽ {[UB]}_{i}$ (⁠ $i \in {1, \dots, n}$ ⁠). The bound constraints are transformed into a set of linear inequality constraints, as in Eq. (20). Note that ARCH is independent of the manner in which they are transformed, but certain existing CHTs, such as MCR, are dependent on this. Table 3 summarizes the features of the constrained problems. Our implementation of the CEC 2006 testbed is available in the repository (https://github.com/naoking158/ARCH).

Table 3:

Details of 24 constrained problems of CEC 2006 testbed: $n$ is the search space dimension; LI, NI, LE, and NE are the number of linear/nonlinear inequality constraints and linear/nonlinear equality constraints, respectively; and $m_{act}$ is the number of active constraints at the optimum.

Prob.	$n$	Type of $f$	LI	NI	LE	NE	$m_{act}$
g01	13	quadratic	9	0	0	0	6
g02	20	nonlinear	0	2	0	0	1
g03	10	polynomial	0	0	0	1	1
g04	5	quadratic	0	6	0	0	2
g05	4	cubic	2	0	0	3	3
g06	2	cubic	0	2	0	0	2
g07	10	quadratic	3	5	0	0	6
g08	2	nonlinear	0	2	0	0	0
g09	7	polynomial	0	4	0	0	2
g10	8	linear	3	3	0	0	6
g11	2	quadratic	0	0	0	1	1
g12	3	quadratic	0	1	0	0	0
g13	5	nonlinear	0	0	0	3	3
g14	10	nonlinear	0	0	3	0	3
g15	3	quadratic	0	0	1	1	2
g16	5	nonlinear	4	34	0	0	4
g17	6	nonlinear	0	0	0	4	4
g18	9	quadratic	0	13	0	0	6
g19	15	nonlinear	0	5	0	0	0
g20	24	linear	0	6	2	12	16
g21	7	linear	0	1	0	5	6
g22	22	linear	0	1	8	11	19
g23	9	linear	0	2	3	1	6
g24	2	linear	0	2	0	0	2

Prob.	$n$	Type of $f$	LI	NI	LE	NE	$m_{act}$
g01	13	quadratic	9	0	0	0	6
g02	20	nonlinear	0	2	0	0	1
g03	10	polynomial	0	0	0	1	1
g04	5	quadratic	0	6	0	0	2
g05	4	cubic	2	0	0	3	3
g06	2	cubic	0	2	0	0	2
g07	10	quadratic	3	5	0	0	6
g08	2	nonlinear	0	2	0	0	0
g09	7	polynomial	0	4	0	0	2
g10	8	linear	3	3	0	0	6
g11	2	quadratic	0	0	0	1	1
g12	3	quadratic	0	1	0	0	0
g13	5	nonlinear	0	0	0	3	3
g14	10	nonlinear	0	0	3	0	3
g15	3	quadratic	0	0	1	1	2
g16	5	nonlinear	4	34	0	0	4
g17	6	nonlinear	0	0	0	4	4
g18	9	quadratic	0	13	0	0	6
g19	15	nonlinear	0	5	0	0	0
g20	24	linear	0	6	2	12	16
g21	7	linear	0	1	0	5	6
g22	22	linear	0	1	8	11	19
g23	9	linear	0	2	3	1	6
g24	2	linear	0	2	0	0	2

7.2 ARCH vs Active-Set ES

In this section, we compare ARCH with the active-set ES (Arnold, 2017).

7.2.1 Settings

We follow the experimental setup of Arnold (2017). Note that the test problems used in Arnold (2017), known as the Michalewicz/Schoenauer (1996) test set are equivalent to the first 11 constrained problems in the CEC 2006 testbed. The initial mean vector $m^{(0)}$ is sampled from the uniform distribution $U [LB, UB]$ ⁠, and then projected onto the boundary of the feasible domain by the repair operator defined by each CHT if the sampled point is infeasible. The initial step size and covariance matrix are set to $σ^{(0)} = 0.2 min {UB - LB}$ and $C^{(0)} = I_{n}$ ⁠, respectively. A run is terminated if the algorithm reaches 1200 iterations or locates a feasible candidate solution $x$ with the objective value $f (x) < f^{*} + ɛ | f^{*} |$ ⁠, where $f^{*}$ is the optimal objective function value, as reported in Liang et al. (2006), and $ɛ \in {10^{- 4}, 10^{- 8}}$ is referred to as the target accuracy. The other parameters of the CMA-ES are set to the default values listed in Table 1.

7.2.2 Results and Discussion

We conduct 100 runs of each algorithm for each problem. The results are summarized in Table 4. The success rate is defined as the number of runs in which the algorithm can reach the target within 1200 iterations, divided by the total number of runs (100).

Table 4:

Median number of function evaluations (NFES), iterations (# Iterations) and success rate (SR) among 100 independent runs for each problem and each algorithm. The results of the active-set ES were obtained from Arnold (2017).

	Target Accuracy	Active-set ES		ARCH
Prob.	$ɛ$	NFES (# Iterations)	SR	NFES	# Iterations	SR
$g 01$	$10^{- 4}$	30	61%	154	$14$	96%
	$10^{- 8}$	30	62%	154	$14$	96%
$g 02$	$10^{- 4}$	–	0%	–	–	0%
	$10^{- 8}$	–	0%	–	–	0%
$g 03$	$10^{- 4}$	463	100%	900	$90$	100%
	$10^{- 8}$	863	100%	1720	$172$	100%
$g 04$	$10^{- 4}$	22	100%	176	22	100%
	$10^{- 8}$	24	100%	176	$22$	100%
$g 05$	$10^{- 4}$	$36$	100%	624	78	97%
	$10^{- 8}$	$82$	100%	1160	145	41%
$g 06$	$10^{- 4}$	5	100%	6	$1$	100%
	$10^{- 8}$	5	100%	6	$1$	100%
$g 07$	$10^{- 4}$	325	100%	1635	$163.5$	100%
	$10^{- 8}$	557	100%	2705	$270.5$	100%
$g 08$	$10^{- 4}$	107	38%	510	$85$	57%
	$10^{- 8}$	210	44%	636	$106$	57%
$g 09$	$10^{- 4}$	307	100%	846	$94$	100%
	$10^{- 8}$	582	100%	1620	$180$	100%
$g 10$	$10^{- 4}$	117	100%	580	$58$	100%
	$10^{- 8}$	$236$	100%	2985	298.5	100%
$g 11$	$10^{- 4}$	25	100%	60	$10$	100%
	$10^{- 8}$	73	100%	258	$43$	100%

	Target Accuracy	Active-set ES		ARCH
Prob.	$ɛ$	NFES (# Iterations)	SR	NFES	# Iterations	SR
$g 01$	$10^{- 4}$	30	61%	154	$14$	96%
	$10^{- 8}$	30	62%	154	$14$	96%
$g 02$	$10^{- 4}$	–	0%	–	–	0%
	$10^{- 8}$	–	0%	–	–	0%
$g 03$	$10^{- 4}$	463	100%	900	$90$	100%
	$10^{- 8}$	863	100%	1720	$172$	100%
$g 04$	$10^{- 4}$	22	100%	176	22	100%
	$10^{- 8}$	24	100%	176	$22$	100%
$g 05$	$10^{- 4}$	$36$	100%	624	78	97%
	$10^{- 8}$	$82$	100%	1160	145	41%
$g 06$	$10^{- 4}$	5	100%	6	$1$	100%
	$10^{- 8}$	5	100%	6	$1$	100%
$g 07$	$10^{- 4}$	325	100%	1635	$163.5$	100%
	$10^{- 8}$	557	100%	2705	$270.5$	100%
$g 08$	$10^{- 4}$	107	38%	510	$85$	57%
	$10^{- 8}$	210	44%	636	$106$	57%
$g 09$	$10^{- 4}$	307	100%	846	$94$	100%
	$10^{- 8}$	582	100%	1620	$180$	100%
$g 10$	$10^{- 4}$	117	100%	580	$58$	100%
	$10^{- 8}$	$236$	100%	2985	298.5	100%
$g 11$	$10^{- 4}$	25	100%	60	$10$	100%
	$10^{- 8}$	73	100%	258	$43$	100%

Comparing the median number of $f$ -calls for reaching the same accuracy, $ɛ$ ⁠, in Table 4, we can observe that the active-set ES achieves the target accuracy with a lower number of $f$ -calls. However, the median number of iterations is lower for ARCH, as it is the population-based approach. As the number of candidate solutions per iteration is increased, we expect the number of iterations to decrease. This is an advantage of ARCH when the candidate solutions can be evaluated on $f$ in parallel.

Focusing on problems g01 and g08, ARCH exhibits higher success rates than the active-set ES. This is possibly because ARCH is combined with the $(μ, λ)$ type ES, rather than the (1+1)-ES. As these problems have local minima, ARCH can reach the target at a higher rate than the active-set ES, which is (1+1)-ES based. None of the algorithms can solve problem g02, because this problem has a strongly multimodal landscape.

Another advantage that does not appear in this experiment is that we incorporate the covariance matrix adaptation, which is desirable for ill-conditioned and non-separable problems.

7.3 ARCH vs AL and MCR

In this section, we compare ARCH with the AL (Atamna et al., 2016) and MCR (de Paula Garcia et al., 2017) techniques.

7.3.1 Settings

In this experiment, we use all 24 problems in the CEC 2006 testbed. Following the default setup in Liang et al. (2006), we regard each run as successful if a feasible candidate solution $x$ with the objective function value $f (x) - f^{*} ⩽ 10^{- 4}$ is located within $5 \times 10^{5}$ $f$ -calls, where $f^{*}$ is the optimal objective value.

Because we consider the constraints to be explicit, it is a natural approach to prepare a set of feasible solutions and begin the optimization process with one feasible solution as the initial search point. The set of feasible solutions is generated as follows. We run the CMA-ES with the following loss function

L : R^{n} \to R

⁠,

L (x) = \sum_{j = 1}^{m} R_{g_{j}} (x), R_{g_{j}} (x) = \sum_{k = 1}^{λ} 1_{{g_{j}^{+} (x_{k}) < g_{j}^{+} (x)}} + \frac{1}{2} \sum_{k = 1}^{λ} 1_{{g_{j}^{+} (x_{k}) = g_{j}^{+} (x)}},

(23)

where

g_{j}^{+} (x) = max (0, g_{j} (x))

⁠. If the candidate solutions are all feasible, they receive the same ranking and the distribution parameter update becomes an unbiased random walk. Then, the CMA-ES tends to produce feasible solutions extensively in a connected feasible subset. To obtain diverse feasible solutions, we run the CMA-ES until

10 n

feasible solutions are generated, and repeat it 50 times with the following initial parameters:

σ^{(0)} = exp (\frac{1}{n} \sum_{i = 1}^{n} ln (\frac{{[UB]}_{i} - {[LB]}_{i}}{5}))

⁠,

C^{(0)} = diag {(\frac{UB - LB}{5 σ^{(0)}})}^{2}

⁠, and

m^{(0)} \sim U [LB, UB]

⁠.

Because multimodal functions exist in the CEC 2006 testbed, we use the BIPOP restart strategy (Hansen, 2009), which updates the population size and initial step size for every restart. Instead of randomly sampling an initial mean vector $m^{(0)}$ from the entire search space, we randomly sample one feasible solution from the above prepared feasible set in each restart. The initial step size and covariance matrix are set to $σ^{(0)} = exp (\frac{1}{n} \sum_{i = 1}^{n} ln (\frac{{[UB]}_{i} - {[LB]}_{i}}{5}))$ and $C^{(0)} = diag {(\frac{UB - LB}{5 σ^{(0)}})}^{2}$ ⁠, respectively, with the step size scaled by the BIPOP strategy in each restart. The other parameters of the CMA-ES are set to the default values listed in Table 1, and we follow the restart condition described in Hansen (2009).

7.3.2 Results and Discussion

All results are summarized in Table 5. The statistical significance is tested using the two-sided Mann–Whitney rank test, with a significance level of $5 % / m$ ⁠, where $m = 42$ is the number of tests (Bonferroni correction).^⁶

Table 5:

Median number of function evaluations (NFES) and success rate (SR) among 25 independent runs for each problem and each algorithm. The subscript of the NFES indicates the average number of restarts for the successful runs. The markers $*$ (AL vs ARCH), $+$ (AL vs MCR), and $\circ$ (MCR vs ARCH) indicate the statistical significance according to the two-sided Mann–Whitney rank test with a significance level of $5 / m %$ ⁠, where $m = 42$ is the number of tests (Bonferroni correction).

	AL		MCR		ARCH
Prob.	NFES	SR	NFES	SR	NFES	SR
g01	–	0%	$320245_{(4.92)}$	100%	$737_{(0.38)}^{\circ}$	100%
g02	–	0%	–	0%	–	0%
g03	–	0%	$10670_{(0.33)}$	100%	$1030_{(0.75)}^{\circ}$	100%
g04	–	0%	$7728_{(0.25)}$	100%	$120_{(0.21)}^{\circ}$	100%
g05	$22056_{(10.00)}^{+}$	100%	$269447_{(47.96)}$	8%	$737_{(1.46)}^{\circ *}$	100%
g06	$8981_{(43.50)}$	4%	$1332_{(0.00)}^{+}$	100%	$90_{(3.62)}^{\circ *}$	100%
g07	$14415_{(2.42)}$	100%	$19860_{(0.17)}$	100%	$2782_{(3.88)}^{\circ *}$	100%
g08	$468_{(0.46)}$	100%	$162_{(0.17)}$	100%	$390_{(3.08)}$	100%
g09	$2079_{(0.00)}^{+}$	100%	$5328_{(0.00)}$	100%	$1791_{(2.42)}^{\circ}$	100%
g10	–	0%	$41780_{(1.12)}$	100%	$2990_{(0.33)}^{\circ}$	100%
g11	$294_{(0.00)}^{+}$	100%	$3204_{(0.88)}$	100%	$54_{(0.17)}^{\circ *}$	100%
g12	$15601_{(6.08)}$	100%	$2870_{(1.38)}^{+}$	100%	$902_{(6.21)}^{\circ *}$	100%
g13	$14208_{(2.12)}$	100%	–	0%	$2686_{(4.96)}^{*}$	100%
g14	$418341_{(16.17)}$	12%	$16830_{(0.46)}^{+}$	100%	$2172_{(2.50)}^{\circ *}$	100%
g15	$2317_{(0.42)}^{+}$	100%	$72239_{(11.17)}$	100%	$175_{(0.92)}^{\circ *}$	100%
g16	$5816_{(0.54)}$	100%	$8664_{(0.46)}$	100%	–	0%
g17	$29400_{(5.92)}$	100%	–	0%	$1814_{(4.17)}^{*}$	100%
g18	$11490_{(2.75)}$	100%	$8810_{(1.58)}$	100%	$3939_{(6.83)}^{\circ *}$	100%
g19	–	0%	$127308_{(1.38)}$	100%	$5772_{(0.08)}^{\circ}$	100%
g20	–	0%	–	0%	–	0%
g21	–	0%	–	0%	–	0%
g22	–	0%	–	0%	–	0%
g23	–	0%	$96078_{(34.38)}$	4%	$9386_{(4.62)}^{\circ}$	100%
g24	$7422_{(3.33)}$	100%	$600_{(0.21)}^{+}$	100%	$54_{(0.38)}^{\circ *}$	100%

	AL		MCR		ARCH
Prob.	NFES	SR	NFES	SR	NFES	SR
g01	–	0%	$320245_{(4.92)}$	100%	$737_{(0.38)}^{\circ}$	100%
g02	–	0%	–	0%	–	0%
g03	–	0%	$10670_{(0.33)}$	100%	$1030_{(0.75)}^{\circ}$	100%
g04	–	0%	$7728_{(0.25)}$	100%	$120_{(0.21)}^{\circ}$	100%
g05	$22056_{(10.00)}^{+}$	100%	$269447_{(47.96)}$	8%	$737_{(1.46)}^{\circ *}$	100%
g06	$8981_{(43.50)}$	4%	$1332_{(0.00)}^{+}$	100%	$90_{(3.62)}^{\circ *}$	100%
g07	$14415_{(2.42)}$	100%	$19860_{(0.17)}$	100%	$2782_{(3.88)}^{\circ *}$	100%
g08	$468_{(0.46)}$	100%	$162_{(0.17)}$	100%	$390_{(3.08)}$	100%
g09	$2079_{(0.00)}^{+}$	100%	$5328_{(0.00)}$	100%	$1791_{(2.42)}^{\circ}$	100%
g10	–	0%	$41780_{(1.12)}$	100%	$2990_{(0.33)}^{\circ}$	100%
g11	$294_{(0.00)}^{+}$	100%	$3204_{(0.88)}$	100%	$54_{(0.17)}^{\circ *}$	100%
g12	$15601_{(6.08)}$	100%	$2870_{(1.38)}^{+}$	100%	$902_{(6.21)}^{\circ *}$	100%
g13	$14208_{(2.12)}$	100%	–	0%	$2686_{(4.96)}^{*}$	100%
g14	$418341_{(16.17)}$	12%	$16830_{(0.46)}^{+}$	100%	$2172_{(2.50)}^{\circ *}$	100%
g15	$2317_{(0.42)}^{+}$	100%	$72239_{(11.17)}$	100%	$175_{(0.92)}^{\circ *}$	100%
g16	$5816_{(0.54)}$	100%	$8664_{(0.46)}$	100%	–	0%
g17	$29400_{(5.92)}$	100%	–	0%	$1814_{(4.17)}^{*}$	100%
g18	$11490_{(2.75)}$	100%	$8810_{(1.58)}$	100%	$3939_{(6.83)}^{\circ *}$	100%
g19	–	0%	$127308_{(1.38)}$	100%	$5772_{(0.08)}^{\circ}$	100%
g20	–	0%	–	0%	–	0%
g21	–	0%	–	0%	–	0%
g22	–	0%	–	0%	–	0%
g23	–	0%	$96078_{(34.38)}$	4%	$9386_{(4.62)}^{\circ}$	100%
g24	$7422_{(3.33)}$	100%	$600_{(0.21)}^{+}$	100%	$54_{(0.38)}^{\circ *}$	100%

The MCR and AL are CHTs that are often used for simulation-based constraints. They do not exploit the fact that the constraints are explicit, but assume that $f$ is defined on the infeasible domain. With the exception of problems g08 and g16, ARCH overwhelmingly outperforms the MCR and AL, which do not explicitly utilize the fact that $g$ -calls are computationally cheaper than $f$ -calls, and perform as many $f$ -calls as $g$ -calls. This indicates that ARCH can efficiently exploit the fact that the constraints are explicit, even if the objective function is defined on the infeasible domain.

ARCH does not reach the target for problems g02, g16, and g20 to g22. For g02, g20, g21, and g22, no algorithms can reach the target value. Problem g02 has a strongly multimodal landscape, and it is difficult to locate the global optimum among many local optima. Problems g20 to g22 have many nonlinear equality constraints, and it is difficult to grasp the global landscape without using the function values outside of the feasible domain. In particular, a feasible solution has not yet been determined for problem g20 (Liang et al., 2006). Problem g16 is a unique problem that can be solved by the AL and MCR, but not ARCH. In ARCH, the internal optimizer (SLSQP) used in the repair operation fails to locate the solution. This may occur if there are too many complex constraints.

8 Conclusions

A constrained continuous optimization problem has been addressed in this study, in which the constraints are assumed to be explicitly written as mathematical expressions and their evaluation time is negligible compared to that of the objective function. We do not assume that the objective function values are defined outside of the feasible domain. Our proposed CHT for the CMA-ES, known as ARCH, can handle explicit and nonlinear constraints. ARCH is designed to be invariant to any element-wise transformation of the objective and constraint functions, as well as to any affine transformation of the search space coordinate, so as to preserve the invariance properties of the underlying CMA-ES. The invariance properties are mathematically proven and empirically validated. To the best of the authors' knowledge, this is the first approach for explicit constraints that is invariant to those transformations.

Two sets of experiments revealed the effectiveness of ARCH. The first demonstrated that the convergence speed of ARCH improved compared with BCH, in the especially ill-conditioned and nonseparable objective function while it is almost the same as in the well-conditioned problem. These results could be attributed to the fact that adapting the coefficient $α$ allows for faster adaptation of the covariance matrix $C$ ⁠, and the affine invariance allows for a search behavior that resembles one on the well-conditioned function. The second experiment revealed that ARCH is overwhelmingly more efficient than the CHTs that do not exploit the explicit constraints. This implies that the explicit constraints should be exploited in the constraint handling, which ARCH effectively achieves. However, note that the CHTs used in the comparison are limited to those available to the CMA-ES. Further comparison with state-of-the-art methods is necessary to show the usefulness of ARCH among existing approaches. On the other hand, the comparison to the active-set (1+1)-ES showed a possibility to further improve the efficacy of ARCH in terms of the number of the objective function calls. Incorporating the idea of the active-set ES into the repair operator in ARCH is a possible direction of a future work.

ARCH is designed to solve the constrained continuous optimization problems where the constraints are all explicit. Box-constrained optimization problems are such examples, where ARCH exhibits superior performance to an existing box constraint handling technique. However, it is often the case that there exist explicit constraints and simulation-based constraints at the same time, where the explicit constraints are prerequisites to a simulator that computes the objective function value and the simulation-based constraint violation values. ARCH is expected to be inefficient for such problems as it requires a substantial amount of constraint evaluations. To deal with the combination of these two types of constraints, one might need to combine an explicit constraint handling such as ARCH and a simulation-based constraint handling such as the multiple constraint ranking (MCR) technique. The invariance properties of ARCH will play an important role to preserve the invariance properties of the entire search algorithm when two CHTs are combined. This is an important direction of future work.

Notes

¹

Our implementation of ARCH is available in a GitHub repository (https://github.com/naoking158/ARCH).

²

Note that we introduce $γ_{σ}$ and $γ_{c}$ that does not appear in a standard formulation in order to treat the initialization effect of the evolution paths and write it short. See Akimoto and Hansen (2020) for more detail.

³

This transformation is applied after the transformation from the equality constraint to the inequality one described in the introduction is performed.

⁴

The invariance of CMA-ES can be proven by using $∥ {(\sqrt{C^{(t)}})}^{- 1} p_{c} ∥$ instead of $∥ p_{σ} ∥$ of the step size update formula in Section 2.1 (note the meaning changes because what is accumulated is not z). However, the step size formula using $p_{σ}$ is generally employed, and no adverse effect of this difference has been empirically observed. Therefore, we use $p_{σ}$ as well in this article.

⁵

We observed in preliminary experiments that the best result was obtained by using Eq. (13).

⁶

Since the tests have been performed excluding the results where no runs reached the target, the number of tests was 42, not $24 \times 3 = 72$ ⁠.

Acknowledgments

This work was supported by the JSPS KAKENHI under Grant Numbers 19H04179 and 19J21892.

References

Akimoto

,

Y.

,

Auger

,

A.

, and

Hansen

,

N.

(

2020

).

Quality gain analysis of the weighted recombination evolution strategy on general convex quadratic functions.

Theoretical Computer Science,

832

:

42

–

67

. Theory of Evolutionary Computation.

Google Scholar

Akimoto

,

Y.

, and

Hansen

,

N.

(

2020

).

Diagonal acceleration for covariance matrix adaptation evolution strategies.

Evolutionary Computation,

28

(

3

):

405

–

435

. PMID: 31120772.

Google Scholar

Arnold

,

D. V

. (

2016

).

An active-set evolution strategy for optimization with known constraints

. In

International Conference on Parallel Problem Solving from Nature

, pp.

192

–

202

.

Google Scholar

Crossref

Arnold

,

D. V

. (

2017

).

Reconsidering constraint release for active-set evolution strategies

. In

Proceedings of the Genetic and Evolutionary Computation Conference

, pp.

665

–

672

.

Google Scholar

Crossref

Arnold

,

D. V.

, and

Hansen

,

N

. (

2012

).

A (⁠

1 + 1

⁠)-CMA-ES for constrained optimisation

. In

Proceedings of the Fourteenth International Conference on Genetic and Evolutionary Computation Conference (GECCO)

, pp.

297

–

304

.

Google Scholar

Crossref

Arnold

,

D. V.

, and

Porter

,

J

. (

2015

).

Towards an augmented Lagrangian constraint handling approach for the (1+1)-ES

. In

Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation (GECCO)

, pp.

249

–

256

.

Google Scholar

Crossref

Atamna

,

A.

,

Auger

,

A.

, and

Hansen

,

N

. (

2016

).

Augmented Lagrangian constraint handling for CMA-ES—Case of a single linear constraint

. In

Proceedings of the 14th International Conference on Parallel Problem Solving from Nature

, pp.

181

–

191

.

Google Scholar

Crossref

Atamna

,

A.

,

Auger

,

A.

, and

Hansen

,

N.

(

2020

).

On invariance and linear convergence of evolution strategies with augmented Lagrangian constraint handling.

Theoretical Computer Science,

832

:

68

–

97

.

Google Scholar

Beyer

,

H.-G

. (

2001

).

The theory of evolution strategies

.

New York

:

Springer Science & Business Media

.

Google Scholar

Crossref

Chocat

,

R.

,

Brevault

,

L.

,

Balesdent

,

M.

, and

Defoort

,

S.

(

2015

).

Modified covariance matrix adaptation–evolution strategy algorithm for constrained optimization under uncertainty, application to rocket design.

International Journal for Simulation and Multidisciplinary Design Optimization,

6

:

A1

.

[PubMed]

Google Scholar

de Paula Garcia

,

R.

,

de Lima

,

B. S. L. P.

,

de Castro Lemonge

,

A. C.

, and

Jacob

,

B. P.

(

2017

).

A rank-based constraint handling technique for engineering design optimization problems solved by genetic algorithms.

Computers & Structures, 187:

77

–

87

.

Google Scholar

Eaton

,

M.

(

2007

).

Chapter 3: The normal distribution on a vector space

, series. Lecture Notes–Monograph Series. Institute of Mathematical Statistics,

53

.

Google Scholar

Hansen

,

N

. (

2009

).

Benchmarking a bi-population CMA-ES on the BBOB-2009 function testbed

. In

Workshop Proceedings of the GECCO Genetic and Evolutionary Computation Conference

, pp.

2389

–

2395

.

Google Scholar

Crossref

Hansen

,

N.

(

2016

).

The CMA evolution strategy: A tutorial

. arXiv:1604.00772

Google Scholar

Hansen

,

N.

,

Muller

,

S. D.

, and

Koumoutsakos

,

P

. (

2003

).

Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)

.

Evolutionary Computation

,

11

(

1

):

1

–

18

.

https://doi.org/10.1162/106365603321828970

[PubMed]

Google Scholar

Crossref

PubMed

Hansen

,

N.

,

Niederberger

,

A. S. P.

,

Guzzella

,

L.

, and

Koumoutsakos

,

P

. (

2009

).

A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion

.

IEEE Transactions on Evolutionary Computation

,

13

(

1

):

180

–

197

.

https://doi.org/10.1109/TEVC.2008.924423

Google Scholar

Crossref

Hansen

,

N.

, and

Ostermeier

,

A

. (

2001

).

Completely derandomized self-adaptation in evolution strategies

.

Evolutionary Computation

,

9

(

2

):

159

–

195

.

https://doi.org/10.1162/106365601750190398

[PubMed]

Google Scholar

Crossref

PubMed

Hansen

,

N.

,

Ros

,

R.

,

Mauny

,

N.

,

Schoenauer

,

M.

, and

Auger

,

A

. (

2011

).

Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems

.

Applied Soft Computing

,

11

(

8

):

5755

–

5769

.

https://doi.org/10.1016/j.asoc.2011.03.001

Google Scholar

Crossref

Hellwig

,

M.

, and

Beyer

,

H

. (

2018

).

A matrix adaptation evolution strategy for constrained real-parameter optimization

. In

2018 IEEE Congress on Evolutionary Computation

, pp.

1

–

8

.

Google Scholar

Crossref

Karafotias

,

G.

,

Hoogendoorn

,

M.

, and

Eiben

,

A. E

. (

2015

).

Parameter control in evolutionary algorithms: Trends and challenges

.

IEEE Transactions on Evolutionary Computation

,

19

(

2

):

167

–

187

.

https://doi.org/10.1109/TEVC.2014.2308294

Google Scholar

Crossref

Kraft

,

D.

(

1988

). A software package for sequential quadratic programming. Technical Report DFVLR-FB 88-28, DLR German Aerospace Center---Institute for Flight Mechanics, Koln, Germany.

Krause

,

O.

, and

Glasmachers

,

T

. (

2015

).

A CMA-ES with multiplicative covariance matrix updates

. In

Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation (GECCO)

, pp.

281

–

288

.

Google Scholar

Crossref

Le Digabel

,

S.

, and

Wild

,

S. M.

(

2015

).

A taxonomy of constraints in simulation-based optimization

. arXiv:1505.07881

Google Scholar

Liang

,

J.

,

Runarsson

,

T. P.

,

Mezura-Montes

,

E.

,

Clerc

,

M.

,

Suganthan

,

P. N.

,

Coello

,

C. C.

, and

Deb

,

K

. (

2006

).

Problem definitions and evaluation criteria for the CEC 2006 special session on constrained real-parameter optimization

.

Journal of Applied Mechanics

,

41

(

8

):

8

–

31

.

Google Scholar

Michalewicz

,

Z.

, and

Schoenauer

,

M

. (

1996

).

Evolutionary algorithms for constrained parameter optimization problems

.

Evolutionary Computation

,

4

(

1

):

1

–

32

.

https://doi.org/10.1162/evco.1996.4.1.1

Google Scholar

Crossref

Oyman

,

A. I.

,

Deb

,

K.

, and

Beyer

,

H

. (

1999

).

An alternative constraint handling method for evolution strategies

. In

Proceedings of the 1999 Congress on Evolutionary Computation

, Vol.

1

, pp.

612

–

619

.

https://doi.org/10.1109/CEC.1999.781989

Google Scholar

Crossref

Runarsson

,

T. P.

, and

Yao

,

X

. (

2000

).

Stochastic ranking for constrained evolutionary optimization

.

IEEE Transactions on Evolutionary Computation

,

4

(

3

):

284

–

294

.

https://doi.org/10.1109/4235.873238

Google Scholar

Crossref

Sakamoto

,

N.

, and

Akimoto

,

Y

. (

2017

).

Modified box constraint handling for the covariance matrix adaptation evolution strategy

. In

Proceedings of the Genetic and Evolutionary Computation Conference (GECCO), Companion

, pp.

183

–

184

.

Google Scholar

Crossref

Sakamoto

,

N.

, and

Akimoto

,

Y

. (

2019

).

Adaptive ranking based constraint handling for explicitly constrained black-box optimization

. In

Proceedings of the Genetic and Evolutionary Computation Conference (GECCO)

, p.

700

.

Google Scholar

Crossref

Spettel

,

P.

,

Beyer

,

H.

, and

Hellwig

,

M

. (

2019

).

A covariance matrix self-adaptation evolution strategy for optimization under linear constraints

.

IEEE Transactions on Evolutionary Computation

,

23

(

3

):

514

–

524

.

https://doi.org/10.1109/TEVC.2018.2871944

Google Scholar

Crossref

Spettel

,

P.

, and

Beyer

,

H.-G

. (

2019

).

A multi-recombinative active matrix adaptation evolution strategy for constrained optimization

.

Soft Computing

,

23

(

16

):

6847

–

6869

.

https://doi.org/10.1007/s00500-018-03736-z

Google Scholar

Crossref

2022

Massachusetts Institute of Technology

Adaptive Ranking-Based Constraint Handling for Explicitly Constrained Black-Box Optimization

Abstract

1 Introduction

2 CMA-ES and Related CHTs

2.1 CMA-ES

2.2 CHTs for ESs

2.2.1 Resampling and Death Penalty

2.2.2 Penalty Function Methods

2.2.3 Ranking-Based Methods

2.2.4 Active Constraint Handling

2.2.5 Other Explicit CHTs

2.3 Formal Classification of CHTs

3 Desired Invariance Properties for Constrained Optimization

3.1 Element-Wise Increasing Transformation of Functions

3.2 Affine Transformation of Search Space Coordinates

4 ARCH

4.1 Repair Operator

4.2 Total Ranking

4.3 Adaptation of Ranking Coefficient

5 Invariance Properties of ARCH

5.1 ARCH on Inner Product Space

5.2 Invariance to Affine Transformation of Search Space Coordinates

5.3 Invariance to Element-Wise Increasing Transformation

5.4 Invariance Properties of Entire Algorithm

6 Experiments on Linearly Constrained Quadratic Problems

6.1 Linearly Constrained Quadratic Minimization Problems

6.2 Settings

6.3 Results and Discussion

7 Experiments on CEC 2006 Constrained Optimization Testbed

7.1 CEC 2006 testbed

7.2 ARCH vs Active-Set ES

7.2.1 Settings

7.2.2 Results and Discussion

7.3 ARCH vs AL and MCR

7.3.1 Settings

7.3.2 Results and Discussion

8 Conclusions

Notes

Acknowledgments

References

Email alerts

Cited By

Related Articles

A product of The MIT Press

MIT Press Direct

Information

MIT Press

Contact Us

This Feature Is Available To Subscribers Only