# Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility

Ronald Aznavourian, Guillaume Demesy, Sébastien Guenneau, Julien Marot

► To cite this version:

Ronald Aznavourian, Guillaume Demesy, Sébastien Guenneau, Julien Marot. Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility. EPJ Applied Metamaterials, 2024, 11, pp.11. 10.1051/epjam/2023003. hal-04568155

HAL Id: hal-04568155

https://hal.science/hal-04568155

Submitted on 3 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility

Ronald Aznavourian, Guillaume Demesy, Sébastien Guenneau, Julien Marot

► To cite this version:

Ronald Aznavourian, Guillaume Demesy, Sébastien Guenneau, Julien Marot. Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility. EPJ Applied Metamaterials, 2024, 11, pp.11. 10.1051/epjam/2023003. hal-04568155

HAL Id: hal-04568155

https://hal.science/hal-04568155v1

Submitted on 3 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility

Ronald Aznavourian<fnref n="1" />,*, Guillaume Demesy<fnref n="1" />, Sébastien Guenneau<fnref n="2" />,<fnref n="3" /> and Julien Marot<fnref n="1" /> ID

Received: 20 July 2023 / Accepted: 27 November 2023

**Abstract.** We revisit the design of cloaks, without resorting to any geometric transform. Cancellation techniques and anomalous resonances have been applied for this purpose. Instead of a deductive reasoning, we propose a novel mono-objective optimization algorithm, namely a ternary grey wolf algorithm, and we adapt a bi-objective optimization algorithm. Firstly, the proposed chaotic ternary grey wolf algorithm searches three-valued spaces for all permittivity values in the cloak while minimizing the summation of a protection criterion and an invisibility criterion. Secondly, a bi-objective genetic algorithm is adapted to find pairs of optimal values of invisibility and protection.

**Keywords:** Mono-objective optimization / bi-objective optimization / cloak design / protection / invisibility

## 1 Introduction

Since the publication of works by the teams of Leonhardt [1] and Pendry [2] in the same issue of the Science magazine over 16–18 years ago, cloaking has become a mature research area optics. It is by now well known that one can design invisibility cloaks via geometric transforms that either lead to the anisotropic heterogeneous material parameters (e.g. rank-2 tensors of permittivity and permeability in optics), see [2], or spatially varying, yet scalar valued, parameters [1]. The latter is achieved through conformal maps, hence constrained to the 2D case, and besides from that the cloak is of infinite extent. There is yet a third route to cloaking, that relaxes the severe material constraints in [2] (notably some infinite anisotropy on the inner boundary of cloak, rooted in the blow up of a point onto a ball of finite extent known as invisibility region, that can only be achieved over a narrow frequency bandwidth in practice [3]), and avoids the infinite extent of the cloak in [1]: so-called carpet cloaking is a combination of the previous two approaches that is based on quasi-conformal grids [4]. This third route requires only some moderate anisotropy, but only achieves invisibility for an object placed on a mirror.

Some works on cloaking focus on the mathematical aspects that are connected to famous inverse problems in particular on electric impedance tomography [5,6] wherein

one wishes to uniquely determine the conductivity within a bounded region, by applying a known static voltage to the surface and recording the resulting current at the boundary (a Dirichlet-to-Neumann map). The Dirichlet-to-Neumann map determines the conductivity [7], but this can only happen if the conductivity is scalar-valued, positive and finite. However, if some of these conditions are not met (e.g. the conductivity is matrix valued) electric impedance tomography fails [8]. This has been exploited to create non-unique conductivities sharing the same boundary measurements [9].

### 1.1 Main contributions

In the present work, we would like to revisit the design of cloaks, without resorting to any geometric transform. Not surprisingly, there is prior work that explored this alternative route, notably through scattering cancellation techniques [10,11] and anomalous resonances [12,13]. However, our rational for the design of cloaks is not based on a deductive reasoning, but rather on some optimization algorithm, and more precisely on some nature-inspired optimizer known as the Grey Wolf Optimizer (GWO). Here again, one may point out former work on design of invisibility cloaks [14,15] based on topology optimization [16]. In [17], some mono-objective genetic optimization algorithm has been applied to estimate the best value, in terms of bias with respect to free-space propagation conditions, of six parameters which define the desired

<footnote n="1">¹ Aix Marseille Univ, CNRS, Centrale Med, Institut Fresnel, Marseille, France</footnote>

<footnote n="2">² UMI 2004 Abraham de Moivre-CNRS, Imperial College London, SW7 2AZ, London, UK</footnote>

<footnote n="3">³ The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, UK</footnote>

<footnote>* e-mail: ronald.aznavourian@fresnel.fr</footnote>
# Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility

Ronald Aznavourian<fnref n="1" />,*, Guillaume Demesy<fnref n="1" />, Sébastien Guenneau<fnref n="2" />,<fnref n="3" /> and Julien Marot<fnref n="1" /> ID

Received: 20 July 2023 / Accepted: 27 November 2023

**Abstract.** We revisit the design of cloaks, without resorting to any geometric transform. Cancellation techniques and anomalous resonances have been applied for this purpose. Instead of a deductive reasoning, we propose a novel mono-objective optimization algorithm, namely a ternary grey wolf algorithm, and we adapt a bi-objective optimization algorithm. Firstly, the proposed chaotic ternary grey wolf algorithm searches three-valued spaces for all permittivity values in the cloak while minimizing the summation of a protection criterion and an invisibility criterion. Secondly, a bi-objective genetic algorithm is adapted to find pairs of optimal values of invisibility and protection.

**Keywords:** Mono-objective optimization / bi-objective optimization / cloak design / protection / invisibility

## 1 Introduction

Since the publication of works by the teams of Leonhardt [1] and Pendry [2] in the same issue of the Science magazine over 16–18 years ago, cloaking has become a mature research area optics. It is by now well known that one can design invisibility cloaks via geometric transforms that either lead to the anisotropic heterogeneous material parameters (e.g. rank-2 tensors of permittivity and permeability in optics), see [2], or spatially varying, yet scalar valued, parameters [1]. The latter is achieved through conformal maps, hence constrained to the 2D case, and besides from that the cloak is of infinite extent. There is yet a third route to cloaking, that relaxes the severe material constraints in [2] (notably some infinite anisotropy on the inner boundary of cloak, rooted in the blow up of a point onto a ball of finite extent known as invisibility region, that can only be achieved over a narrow frequency bandwidth in practice [3]), and avoids the infinite extent of the cloak in [1]: so-called carpet cloaking is a combination of the previous two approaches that is based on quasi-conformal grids [4]. This third route requires only some moderate anisotropy, but only achieves invisibility for an object placed on a mirror.

Some works on cloaking focus on the mathematical aspects that are connected to famous inverse problems in particular on electric impedance tomography [5,6] wherein

one wishes to uniquely determine the conductivity within a bounded region, by applying a known static voltage to the surface and recording the resulting current at the boundary (a Dirichlet-to-Neumann map). The Dirichlet-to-Neumann map determines the conductivity [7], but this can only happen if the conductivity is scalar-valued, positive and finite. However, if some of these conditions are not met (e.g. the conductivity is matrix valued) electric impedance tomography fails [8]. This has been exploited to create non-unique conductivities sharing the same boundary measurements [9].

### 1.1 Main contributions

In the present work, we would like to revisit the design of cloaks, without resorting to any geometric transform. Not surprisingly, there is prior work that explored this alternative route, notably through scattering cancellation techniques [10,11] and anomalous resonances [12,13]. However, our rational for the design of cloaks is not based on a deductive reasoning, but rather on some optimization algorithm, and more precisely on some nature-inspired optimizer known as the Grey Wolf Optimizer (GWO). Here again, one may point out former work on design of invisibility cloaks [14,15] based on topology optimization [16]. In [17], some mono-objective genetic optimization algorithm has been applied to estimate the best value, in terms of bias with respect to free-space propagation conditions, of six parameters which define the desired

<footnote n="1">¹ Aix Marseille Univ, CNRS, Centrale Med, Institut Fresnel, Marseille, France</footnote>

<footnote n="2">² UMI 2004 Abraham de Moivre-CNRS, Imperial College London, SW7 2AZ, London, UK</footnote>

<footnote n="3">³ The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, UK</footnote>

<footnote>* e-mail: ronald.aznavourian@fresnel.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the scattering problem: Perfectly Matched Layers (PMLs, light grey) surround the free-space region $S_2$ (cyan) that contains the source, the cloak (blue) and the protected region $S_1$ (orange).</figcaption></figure>

cloak. However, we stress that the nature of the optimization algorithms we shall use here is radically different: We aim at estimating a very elevated number of parameters, compared to some previous works such as [17], and we aim at minimizing two contradictory criteria instead of one.

## 1.2 Layout of the paper

In Section 2, we introduce the cloak design problem, pointing out the need for the minimization of two antagonist criteria: an invisibility criterion and a protection criterion. In Section 3, we provide a state-of-the-art of mono-objective and bi-objective optimization algorithms. We focus on the mono-objective grey wolf optimizer, and on the bi-objective non dominated sorting genetic algorithm. In Section 4, we propose a novel ternary version of the grey wolf algorithm, which is dedicated to search spaces with three values. We name it chaotic ternary grey wolf optimizer (CTGWO). In Section 5, we present the results obtained on cloak design with the mono-objective approach, including CTGWO, and with the bi-objective approach involving Non-dominated Sorting Genetic Algorithm - II (NSGA-II). In Section 6 we discuss the results obtained. We point out the superiority of CTGWO over comparative algorithms in the mono-objective approach; and we emphasize the interest of the bi-objective approach for an end-user. Conclusions are drawn in Section 7.

## 1.3 Notations

The following notations are used throughout the paper: manifolds are denoted by blackboard bold, A, matrices by boldface uppercase roman, A. Vectors are denoted by boldface lowercase roman, a, and scalars by lowercase or uppercase roman, a, b or A. The P scalar components of a vector a are accessed via $a^1, a^2, ..., a^P$, such that $a = [a^1, a^2, ..., a^P]^T$ with superscript T denoting the transpose. The interval of real values between scalars a and b is denoted by $[a : b]$ with square brackets. A set of values is denoted by $\{a,...,b\}$ with curly brackets.

The symbol ◦ denotes the Hadamard (also called component-wise) product of two vectors: $a \circ b$ is a vector whose K components are equal to $a_1b_1, a_2b_2, ..., a_Kb_K$. The symbol $|a|$ means element-wise absolute value and is a vector whose K components are equal to $|a_1|, |a_2|, ..., |a_K|$.

## 2 Cloak design problem: definition of experimental conditions, protection and invisibility criteria

We consider the 2D scattering problem sketched in Figure 1a. A point source is located at $(x_s, y_s)$ and radiates from freespace (in blue color) in the vicinity of the yellow zone $S_1$ which is the area to be protected. This yellow region is cloaked by the green rectangular zone. For computational purposes, the blue freespace zone is surrounded by Perfectly Matched Layers (PMLs) that model an unbounded medium [18]. The scalar scattering problem amounts to finding the total scalar field u such that:

$$ \mathrm{div}[\sigma \mathbf{grad} \, u] + k_0^2 \chi u = \delta_S, \quad (1) $$

where $k_0$ is the freespace wavenumber (associated with the unbounded region of space outside of the object and its surrounding cloak and corresponding to a freespace wavelength of $\lambda_0 = 2\pi/k_0$), and $\sigma$ and $\chi$ represent the scalar material properties. In the context of acoustic pressure waves in isotropic non-viscous fluids, $\sigma$ is the inverse of density and $\chi$ the inverse of compressibility, and $p$ is the amplitude of the pressure wave. For anti-plane shear waves in isotropic solids, $\sigma$ is the shear modulus and $\chi$ the density, whereas $u$ stands for the component of the displacement field perpendicular to the (xy)-plane. Finally, for electromagnetic waves in transverse magnetic (TM) polarization $\sigma$ is the inverse of the relative permittivity and $\chi$ the relative permeability. The relative permittivity will be denoted by $\epsilon$ in the rest of the paper. This is when $u$ represents the component of the magnetic field perpendicular to the (xy)-plane. The roles of permittivity and permeability are interchanged in transverse electric (TE) polarization, whereby the electric field is perpendicular to (xy)-plane. All these physical setups are equivalent from a mathematical standpoint. In the rest of the paper, the considered wavelength is $\lambda = 500$ nm the side of a triangle mesh has length $\lambda/6$ in the freespace; and $\lambda/12$ in the cloak and the protected zone.

In what follows, we focus on the TM case. In this scalar Helmholtz equation (1), one can still consider a certain form of anisotropy [19], provided that $\sigma$ can be written as a $2 \times 2$ matrix.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the scattering problem: Perfectly Matched Layers (PMLs, light grey) surround the free-space region $S_2$ (cyan) that contains the source, the cloak (blue) and the protected region $S_1$ (orange).</figcaption></figure>

cloak. However, we stress that the nature of the optimization algorithms we shall use here is radically different: We aim at estimating a very elevated number of parameters, compared to some previous works such as [17], and we aim at minimizing two contradictory criteria instead of one.

## 1.2 Layout of the paper

In Section 2, we introduce the cloak design problem, pointing out the need for the minimization of two antagonist criteria: an invisibility criterion and a protection criterion. In Section 3, we provide a state-of-the-art of mono-objective and bi-objective optimization algorithms. We focus on the mono-objective grey wolf optimizer, and on the bi-objective non dominated sorting genetic algorithm. In Section 4, we propose a novel ternary version of the grey wolf algorithm, which is dedicated to search spaces with three values. We name it chaotic ternary grey wolf optimizer (CTGWO). In Section 5, we present the results obtained on cloak design with the mono-objective approach, including CTGWO, and with the bi-objective approach involving Non-dominated Sorting Genetic Algorithm - II (NSGA-II). In Section 6 we discuss the results obtained. We point out the superiority of CTGWO over comparative algorithms in the mono-objective approach; and we emphasize the interest of the bi-objective approach for an end-user. Conclusions are drawn in Section 7.

## 1.3 Notations

The following notations are used throughout the paper: manifolds are denoted by blackboard bold, A, matrices by boldface uppercase roman, A. Vectors are denoted by boldface lowercase roman, a, and scalars by lowercase or uppercase roman, a, b or A. The P scalar components of a vector a are accessed via $a^1, a^2, ..., a^P$, such that $a = [a^1, a^2, ..., a^P]^T$ with superscript T denoting the transpose. The interval of real values between scalars a and b is denoted by $[a : b]$ with square brackets. A set of values is denoted by $\{a,...,b\}$ with curly brackets.

The symbol ◦ denotes the Hadamard (also called component-wise) product of two vectors: $a \circ b$ is a vector whose K components are equal to $a_1b_1, a_2b_2, ..., a_Kb_K$. The symbol $|a|$ means element-wise absolute value and is a vector whose K components are equal to $|a_1|, |a_2|, ..., |a_K|$.

## 2 Cloak design problem: definition of experimental conditions, protection and invisibility criteria

We consider the 2D scattering problem sketched in Figure 1a. A point source is located at $(x_s, y_s)$ and radiates from freespace (in blue color) in the vicinity of the yellow zone $S_1$ which is the area to be protected. This yellow region is cloaked by the green rectangular zone. For computational purposes, the blue freespace zone is surrounded by Perfectly Matched Layers (PMLs) that model an unbounded medium [18]. The scalar scattering problem amounts to finding the total scalar field u such that:

$$ \mathrm{div}[\sigma \operatorname{grad} u] + k_0^2 \chi u = \delta_S, \quad (1) $$

where $k_0$ is the freespace wavenumber (associated with the unbounded region of space outside of the object and its surrounding cloak and corresponding to a freespace wavelength of $\lambda_0 = 2\pi/k_0$), and $\sigma$ and $\chi$ represent the scalar material properties. In the context of acoustic pressure waves in isotropic non-viscous fluids, $\sigma$ is the inverse of density and $\chi$ the inverse of compressibility, and $p$ is the amplitude of the pressure wave. For anti-plane shear waves in isotropic solids, $\sigma$ is the shear modulus and $\chi$ the density, whereas $u$ stands for the component of the displacement field perpendicular to the (xy)-plane. Finally, for electromagnetic waves in transverse magnetic (TM) polarization $\sigma$ is the inverse of the relative permittivity and $\chi$ the relative permeability. The relative permittivity will be denoted by $\epsilon$ in the rest of the paper. This is when $u$ represents the component of the magnetic field perpendicular to the (xy)-plane. The roles of permittivity and permeability are interchanged in transverse electric (TE) polarization, whereby the electric field is perpendicular to (xy)-plane. All these physical setups are equivalent from a mathematical standpoint. In the rest of the paper, the considered wavelength is $\lambda = 500$ nm the side of a triangle mesh has length $\lambda/6$ in the freespace; and $\lambda/12$ in the cloak and the protected zone.

In what follows, we focus on the TM case. In this scalar Helmholtz equation (1), one can still consider a certain form of anisotropy [19], provided that $\sigma$ can be written as a $2 \times 2$ matrix.
Indeed, even if isotropic cloak and protected area only are considered here, the implementation of the PMLS relies on anisotropy (and absorption), see [18]. As in standard topology optimization, our design space is constructed on the mesh of the cloak shown in Figure 1. Each triangular element constitutes a voxel that can be filled with a particular material whose isotropic physical properties are represented by the scalar quantity $\sigma$. The practical Finite Elements discretization and implementation of the problem follow quite closely those described in reference [20] and its associated ONELAB open source tutorials found in reference [21]. In short, a first dummy run allows to initialize a table of discrete values of $\gamma$ for a particular mesh (one entry of the table corresponds to one triangle of the mesh in the design region), and this table is used to control discrete $\gamma$ values throughout the optimization procedure using the GetDP, GmshRead and ScalarField functions [22]. Then, equation (1) is solved using second order Lagrange elements in a standard manner.

Finally, we are in a position to define the optimization criteria. The first one, $C_1$ is a protection criterion, the integral over the region to be protected of the square norm of the field:

$$C_1 = \frac{1}{|S_1|} \int_{S_1} |u|^2 dS, \qquad (2)$$

where $|S_1|$ is the area of $S_1$ in Figure 1.

The second one, $C_2$ is an invisibility criterion, the integral over the freespace region of interest of the square norm of the difference between the field and a reference field denoted $u_0$

$$C_2 = \frac{1}{|S_2|} \int_{S_2} |u - u_0|^2 dS, \qquad (3)$$

where $|S_2|$ is the area of $S_2$ in Figure 1 and $u_0$ is the freespace solution to the problem (i.e. the scalar Green function of the problem without any cloak or region to protect).

Our goal is to minimize $C_1$ and $C_2$ to achieve both protection and invisibility. Intuitively, we expect there should be some trade-off between protection (i.e. a vanishing field $u$ inside the yellow region) and invisibility (i.e. $u$ as close as possible to $u_0$ in the blue region). For instance, surrounding the yellow region by an infinite conducting boundary would ensure a vanishing field $u$ inside the orange region, but $u$ and $u_0$ would be very different in freespace green region $S_2$ due to a large scattering. On the other hand, if we consider freespace in the green region, then $S_1$, the cloak and $S_2$ are impregnated matched at the interfaces between all three regions and thus $u = u_0$ in $S_2$ (the cloak and $S_1$ being transparent). But in that case there is no protection at all in $S_1$. Depending upon their need, cloak designers might just wish to give more weight to criterion $C_1$ or $C_2$ and $C_1$ and $C_2$ depend both on $P$ parameters, where $P$ is equal to the number of voxels (equivalently triangles) in the cloak. These parameters, denoted by $K^1, K^2, ..., K^i, ..., K^P$ take their values in a

so-called ‘search space’, which can be either discrete or continuous. In this paper we will firstly consider the realistic case where three possible permittivity values can be associated with each triangle. These values correspond to three different materials and yield the search space $\{7,10,12\}$ Equivalently, once these values are set, we may access them through the index values $\{0,1,2\}$. This will compound our ‘ternary search space’. Secondly, we will perform a study which is less realistic: in our simulations the permittivity for each voxel may be any real value in $[7:12]$ (up to seven decimal digits).

In a nutshell, we notice that we face a bi-objective problem with either ternary, or continuous search spaces. In the rest of the paper, the criteria $C_1$ and $C_2$ will also be denoted by $f_1(\mathbf{x}):\mathbb{R}^P \to \mathbb{R}_+$ and $f_2(\mathbf{x}):\mathbb{R}^P \to \mathbb{R}_+$. Vector $\mathbf{x}$ contains the parameter values $K^1, K^2, ..., K^i, ..., K^P$.

## 3 Global single-objective and bi-objective optimization methods

Optimization algorithms are meant to retrieve the location of the minimum value reached with a set of parameters. In the case of single-objective optimization, only one function is considered for optimization; in the case of bi-objective optimization, two functions should be minimized simultaneously.

Mono-objective problems are encountered in wave propagation phenomena, for instance for the structural optimization of a two-dimensional photonic crystal [23], or to perform the inversion of ellipsometric data [24], thanks to good convergence properties of bio-inspired optimization methods. Multi-objective problems are faced when a tradeoff should be found between two antagonistic objectives, in particular for material design as in [25], but also when considering electromagnetic wave propagation, for instance when designing optical diffracting devices [26].

We assume that $P$ parameters should be estimated: $K^1, K^2, ..., K^i, ..., K^P$ where $P \ge 1$. We remind that, as explained in Section 2, in the considered problem, $P$ is equal to the number of voxels in the cloak. The following notations will be used:

- $P$ is the number of expected parameters, which are indexed with $i$.

- $n$ denotes one iteration and $n_{max}$ the maximum allowed number of iterations.

- $f(\cdot)$ is a function to be optimized, also called the criterion. It depends on the $P$ parameters mentioned above. In this paper, unless specified, minimization problems are considered.

In the case of a single-objective optimization, there is one function $f(\cdot)$ to be minimized. In the case of bi-objective optimization, there are two functions $f_1(\cdot)$ and $f_2(\cdot)$ to be minimized.

Both GWO and NSGA-II are agent-based algorithms $i$. $-\mathbf{x}_q(n)$ is a vector corresponding to an agent $q = 1,...,Q$ at iteration $n$. It takes the form of a vector with a $P$-tuple of tested values $\mathbf{x}_q(n)=[K^1,K^2,...,K^P]^T$.
Indeed, even if isotropic cloak and protected area only are considered here, the implementation of the PMLS relies on anisotropy (and absorption), see [18]. As in standard topology optimization, our design space is constructed on the mesh of the cloak shown in Figure 1. Each triangular element constitutes a voxel that can be filled with a particular material whose isotropic physical properties are represented by the scalar quantity $\sigma$. The practical Finite Elements discretization and implementation of the problem follow quite closely those described in reference [20] and its associated ONELAB open source tutorials found in reference [21]. In short, a first dummy run allows to initialize a table of discrete values of $\gamma$ for a particular mesh (one entry of the table corresponds to one triangle of the mesh in the design region), and this table is used to control discrete $\gamma$ values throughout the optimization procedure using the GetDP, GmshRead and ScalarField functions [22]. Then, equation (1) is solved using second order Lagrange elements in a standard manner.

Finally, we are in a position to define the optimization criteria. The first one, $C_1$ is a protection criterion, the integral over the region to be protected of the square norm of the field:

$$C_1 = \frac{1}{|S_1|} \int_{S_1} |u|^2 dS, \qquad (2)$$

where $|S_1|$ is the area of $S_1$ in Figure 1.

The second one, $C_2$ is an invisibility criterion, the integral over the freespace region of interest of the square norm of the difference between the field and a reference field denoted $u_0$

$$C_2 = \frac{1}{|S_2|} \int_{S_2} |u - u_0|^2 dS, \qquad (3)$$

where $|S_2|$ is the area of $S_2$ in Figure 1 and $u_0$ is the freespace solution to the problem (i.e. the scalar Green function of the problem without any cloak or region to protect).

Our goal is to minimize $C_1$ and $C_2$ to achieve both protection and invisibility. Intuitively, we expect there should be some trade-off between protection (i.e. a vanishing field $u$ inside the yellow region) and invisibility (i.e. $u$ as close as possible to $u_0$ in the blue region). For instance, surrounding the yellow region by an infinite conducting boundary would ensure a vanishing field $u$ inside the orange region, but $u$ and $u_0$ would be very different in freespace green region $S_2$ due to a large scattering. On the other hand, if we consider freespace in the green region, then $S_1$, the cloak and $S_2$ are impregnated matched at the interfaces between all three regions and thus $u = u_0$ in $S_2$ (the cloak and $S_1$ being transparent). But in that case there is no protection at all in $S_1$. Depending upon their need, cloak designers might just wish to give more weight to criterion $C_1$ or $C_2$ and $C_1$ and $C_2$ depend both on $P$ parameters, where $P$ is equal to the number of voxels (equivalently triangles) in the cloak. These parameters, denoted by $K^1, K^2, ..., K^i, ..., K^P$ take their values in a

so-called ‘search space’, which can be either discrete or continuous. In this paper we will firstly consider the realistic case where three possible permittivity values can be associated with each triangle. These values correspond to three different materials and yield the search space $\{7,10,12\}$ Equivalently, once these values are set, we may access them through the index values $\{0,1,2\}$. This will compound our ‘ternary search space’. Secondly, we will perform a study which is less realistic: in our simulations the permittivity for each voxel may be any real value in $[7:12]$ (up to seven decimal digits).

In a nutshell, we notice that we face a bi-objective problem with either ternary, or continuous search spaces. In the rest of the paper, the criteria $C_1$ and $C_2$ will also be denoted by $f_1(\mathbf{x}):\mathbb{R}^P \to \mathbb{R}_+$ and $f_2(\mathbf{x}):\mathbb{R}^P \to \mathbb{R}_+$. Vector $\mathbf{x}$ contains the parameter values $K^1, K^2, ..., K^i, ..., K^P$.

## 3 Global single-objective and bi-objective optimization methods

Optimization algorithms are meant to retrieve the location of the minimum value reached with a set of parameters. In the case of single-objective optimization, only one function is considered for optimization; in the case of bi-objective optimization, two functions should be minimized simultaneously.

Mono-objective problems are encountered in wave propagation phenomena, for instance for the structural optimization of a two-dimensional photonic crystal [23], or to perform the inversion of ellipsometric data [24], thanks to good convergence properties of bio-inspired optimization methods. Multi-objective problems are faced when a tradeoff should be found between two antagonistic objectives, in particular for material design as in [25], but also when considering electromagnetic wave propagation, for instance when designing optical diffracting devices [26].

We assume that $P$ parameters should be estimated: $K^1, K^2, ..., K^i, ..., K^P$ where $P \ge 1$. We remind that, as explained in Section 2, in the considered problem, $P$ is equal to the number of voxels in the cloak. The following notations will be used:

- $P$ is the number of expected parameters, which are indexed with $i$.

- $n$ denotes one iteration and $n_{max}$ the maximum allowed number of iterations.

- $f(\cdot)$ is a function to be optimized, also called the criterion. It depends on the $P$ parameters mentioned above. In this paper, unless specified, minimization problems are considered.

In the case of a single-objective optimization, there is one function $f(\cdot)$ to be minimized. In the case of bi-objective optimization, there are two functions $f_1(\cdot)$ and $f_2(\cdot)$ to be minimized.

Both GWO and NSGA-II are agent-based algorithms $i$. $-\mathbf{x}_q(n)$ is a vector corresponding to an agent $q = 1,...,Q$ at iteration $n$. It takes the form of a vector with a $P$-tuple of tested values $\mathbf{x}_q(n)=[K^1,K^2,...,K^P]^T$.
In Section 3.1, we give a background on a single-objective optimization: the grey wolf algorithm (GWO). In Section 3.2, we give a background on a bi-objective optimization algorithm, namely non-dominated sorting genetic algorithm (NSGA) [27] and a fast version (NSGA-II) [28].

## 3.1 Single-objective optimization: Grey Wolf Optimizer

The GWO is a nature-inspired optimizer based on the observation of the social life of grey wolves in nature [29]. In this method an agent is called a wolf. It simulates the common behaviour and hunting strategies of grey wolves in their environment. The seminal GWO searches a continuous space [29]. Among the search agents, there are three leaders α, β, and δ All other agents are the ω wolves.

$x_α(n)$, $x_β(n)$, and $x_δ(n)$ denote the position of the leaders α, β and δ respectively, at iteration n.

The position of any wolf at iteration n+1 is calculated as:

$$x_q(n+1) = \frac{1}{3}(y_α(n) + y_β(n) + y_δ(n)). \quad (4)$$

It results from the equal contribution of the α, β and δ wolves. These contributions are computed at each iteration iter as follows, for any leader l either α, β or δ:

$$y_l(n) = x_l(n) - b \circ d_l(n), \quad (5)$$

with: $d_l(n) = |c \circ x_l(n) - x_q(n)|, |c|$ denoting absolute value.

The vectors b and c are calculated as $b = 2a \circ r_1 - a$ and $c = 2r_2$. In these expressions, vectors $r_1, r_2$ have random components between 0 and 1.

For the i-th parameter ($i=1,...,P$):

The component $d_l^i(n)$ of $d_l(n)$ is provided by:

$$b^i = 2ar_1 - a, \quad (6)$$

the same whatever i.

The component $d_l^i(n)$ of $d_l(n)$ is provided by:

$$d_l^i(n) = |2r_2x_l^i - x_q^i(n)|, \quad (7)$$

where $r_1$ and $r_2$ are two random values between 0 and 1; $x_q^i(n)$ is the i-th component of the q-th agent at iteration n; $x_l^i$ is the i-th component of leader l.

The component $y_l^i(n)$ of $y_l(n)$ is provided by:

$$y_l^i(n) = x_l^i - b^i d_l^i(n). \quad (8)$$

During the hunt, the wolves firstly diverge from each other to search for the prey, or equivalently to encircle it. Secondly, they converge to kill the prey. This is mathematically modeled through the deterministic vector a. The components of vector a are all equal to a, a scalar value which is a key parameter in the algorithm. When $a > 1$ the search agents are obliged to diverge from the prey: this is the exploration phase. Conversely, when $a \le 1$, the

search agents are obliged to attack towards the prey: this is the exploitation phase. In the vanilla version of GWO [29], the key parameter a decreased regularly from 2 to 0:

$$a = 2 \left( 1 - \frac{n}{n_{\max}} \right). \quad (9)$$

In more recent works, various expressions have been proposed for a such as a quadratic [30] or adaptive [31] function. Whatever the version [30, 31] the exploration phase lasts until $a=1$, then the exploitation phase lasts from $a=1$ to $a=0$. Storing all values of $f(x_α(n))$ across all iterations from 1 to $n_{max}$ yields a so-called convergence curve. The outcomes of a single-objective optimization method are mainly the solution $x_α(n_{max})$, but also the convergence curve.

## 3.2 Bi-objective optimization: non-dominated sorting genetic algorithm

Non-dominated sorting genetic algorithms (NSGA and NSGA-II) are inspired by Darwin's rules of evolution. In this method an agent is called a chromosome. The interest of multi-objective optimization has already been emphasized in physical phenomenon modeling for instance in [32], where NSGA-II is used to model turbulence; or in the design of finite 3D periodic structures [33], etc. The main steps of the algorithm are as follows: generate a random population of chromosomes, calculate function values $f_1$ and $f_2$ for each chromosome, sort chromosomes in the population, choose parents in the next generation by tournament algorithm, generate children by crossover and mutation, extract a new generation through ranking, and repeat the process from parent choice. The expected outcome of NSGA-II is different from the outcome of GWO: For a single-objective method, we will represent the results as convergence curves, and the solution is a single set of values extracted from the search spaces. We will try to minimize our cost function so that our objective tends as much as possible towards zero. It is up to the user to determine a threshold value for our cost function, from which we will retrieve one optimal solution. For a bi-objective method, we will no longer have convergence curves, but Pareto fronts. The solution is composed of several sets of values extracted from the search space and located on the Pareto front. The principle of the Pareto front is that we will represent the value of the first cost function on the horizontal axis, and the value of the second cost function on the vertical axis. Thus, visually, we will see very quickly if a solution favors either $f_1$ or $f_2$. The main method to estimate the quality of a solution is the 'domination' principle. In fact, for each solution, we will calculate the distance between the point of origin of coordinates (0,0), and the considered solution, with its coordinates. Of course, the solution which has the smallest distance will 'dominate' the solution which has a larger distance. Considering two solutions $x_{q1}$ and $x_{q2}$, we can say that $x_{q1}$ 'dominates' $x_{q2}$, if the following condition is respected:

$$((f_1(x_{q1}) \le f_1(x_{q2})) \text{ and } (f_2(x_{q1}) \le f_2(x_{q2}))),$$
In Section 3.1, we give a background on a single-objective optimization: the grey wolf algorithm (GWO). In Section 3.2, we give a background on a bi-objective optimization algorithm, namely non-dominated sorting genetic algorithm (NSGA) [27] and a fast version (NSGA-II) [28].

## 3.1 Single-objective optimization: Grey Wolf Optimizer

The GWO is a nature-inspired optimizer based on the observation of the social life of grey wolves in nature [29]. In this method an agent is called a wolf. It simulates the common behaviour and hunting strategies of grey wolves in their environment. The seminal GWO searches a continuous space [29]. Among the search agents, there are three leaders $\alpha, \beta$, and $\delta$. All other agents are the $\omega$ wolves.

$\mathbf{x}_{\alpha}(n)$, $\mathbf{x}_{\beta}(n)$, and $\mathbf{x}_{\delta}(n)$ denote the position of the leaders $\alpha, \beta$ and $\delta$ respectively, at iteration $n$.

The position of any wolf at iteration $n+1$ is calculated as:

$$ \mathbf{x}_{q}(n+1) = \frac{1}{3}(\mathbf{y}_{\alpha}(n) + \mathbf{y}_{\beta}(n) + \mathbf{y}_{\delta}(n)). \quad (4) $$

It results from the equal contribution of the $\alpha, \beta$ and $\delta$ wolves. These contributions are computed at each iteration $i$ as follows, for any leader $l$ either $\alpha, \beta$ or $\delta$:

$$ \mathbf{y}_l(n) = \mathbf{x}_l(n) - \mathbf{b} \circ \mathbf{d}_l(n), \quad (5) $$

with: $\mathbf{d}_l(n) = |\mathbf{c} \circ \mathbf{x}_l(n) - \mathbf{x}_q(n)|, |\cdot|$ denoting absolute value.

The vectors $\mathbf{b}$ and $\mathbf{c}$ are calculated as $\mathbf{b} = 2\mathbf{a} \circ \mathbf{r}_1 - \mathbf{a}$ and $\mathbf{c} = 2\mathbf{r}_2$. In these expressions, vectors $\mathbf{r}_1, \mathbf{r}_2$ have random components between 0 and 1.

For the $i^{th}$ parameter ($i=1,...,P$):

The component $d_i^i$ of $\mathbf{d}_l(n)$ is provided by:

$$ b^i = 2ar_1 - a, \quad (6) $$

the same whatever $i$.

The component $d_i^i(n)$ of $\mathbf{d}_l(n)$ is provided by:

$$ d_l^i(n) = |2r_2x_l^i - x_q^i(n)|, \quad (7) $$

where $r_1$ and $r_2$ are two random values between 0 and 1; $x_q^i(n)$ is the $i^{th}$ component of the $q^{th}$ agent at iteration $n$; $x_l^i$ is the $i^{th}$ component of leader $l$.

The component $y_l^i(n)$ of $\mathbf{y}_l(n)$ is provided by:

$$ y_l^i(n) = x_l^i - b^i d_l^i(n). \quad (8) $$

During the hunt, the wolves firstly diverge from each other to search for the prey, or equivalently to encircle it. Secondly, they converge to kill the prey. This is mathematically modeled through the deterministic vector $\mathbf{a}$. The components of vector $\mathbf{a}$ are all equal to $a$, a scalar value which is a key parameter in the algorithm. When $a>1$ the search agents are obliged to diverge from the prey: this is the exploration phase. Conversely, when $a \le 1$, the

search agents are obliged to attack towards the prey: this is the exploitation phase. In the vanilla version of GWO [29], the key parameter $a$ decreased regularly from 2 to 0:

$$ a = 2 \left( 1 - \frac{n}{n_{\max}} \right). \quad (9) $$

In more recent works, various expressions have been proposed for $a$ such as a quadratic [30] or adaptive [31] function. Whatever the version [30, 31] the exploration phase lasts until $a=1$, then the exploitation phase lasts from $a=1$ to $a=0$. Storing all values of $f(\mathbf{x}_{\alpha}(n))$ across all iterations from 1 to $n_{\max}$ yields a so-called convergence curve. The outcomes of a single-objective optimization method are mainly the solution $\mathbf{x}_{\alpha}(n_{\max})$, but also the convergence curve.

## 3.2 Bi-objective optimization: non-dominated sorting genetic algorithm

Non-dominated sorting genetic algorithms (NSGA and NSGA-II) are inspired by Darwin's rules of evolution. In this method an agent is called a chromosome. The interest of multi-objective optimization has already been emphasized in physical phenomenon modeling for instance in [32], where NSGA-II is used to model turbulence; or in the design of finite 3D periodic structures [33], etc. The main steps of the algorithm are as follows: generate a random population of chromosomes, calculate function values $f_1$ and $f_2$ for each chromosome, sort chromosomes in the population, choose parents in the next generation by tournament algorithm, generate children by crossover and mutation, extract a new generation through ranking, and repeat the process from parent choice. The expected outcome of NSGA-II is different from the outcome of GWO: For a single-objective method, we will represent the results as convergence curves, and the solution is a single set of values extracted from the search spaces. We will try to minimize our cost function so that our objective tends as much as possible towards zero. It is up to the user to determine a threshold value for our cost function, from which we will retrieve one optimal solution. For a bi-objective method, we will no longer have convergence curves, but Pareto fronts. The solution is composed of several sets of values extracted from the search space and located on the Pareto front. The principle of the Pareto front is that we will represent the value of the first cost function on the horizontal axis, and the value of the second cost function on the vertical axis. Thus, visually, we will see very quickly if a solution favors either $f_1$ or $f_2$. The main method to estimate the quality of a solution is the 'domination' principle. In fact, for each solution, we will calculate the distance between the point of origin of coordinates (0,0), and the considered solution, with its coordinates. Of course, the solution which has the smallest distance will 'dominate' the solution which has a larger distance. Considering two solutions $\mathbf{x}_{q1}$ and $\mathbf{x}_{q2}$, we can say that $\mathbf{x}_{q1}$ 'dominates' $\mathbf{x}_{q2}$, if the following condition is respected:

$$ ((f_1(\mathbf{x}_{q1}) \le f_1(\mathbf{x}_{q2})) \text{ and } (f_2(\mathbf{x}_{q1}) \le f_2(\mathbf{x}_{q2}))), $$
and

$$((f_1(\mathbf{x}_{q1}) < f_1(\mathbf{x}_{q2})) \text{ or } (f_2(\mathbf{x}_{q1}) < f_2(\mathbf{x}_{q2}))).$$

It is up to the user to choose the best solution(s) to keep. Indeed, the user may very well want to favour one of the two objectives, or seek the best compromise between the two objectives.

Actually, in the last decades, different methods to determine the ‘domination’ principle have been proposed, and then, the ‘non-domination’ principle emerged, notably thanks to Srinivas and Deb [27] who proposed the NSGA algorithm [27] and then an improved version, called NSGA-II [28]. This fast sorting method by ‘non-domination’ has been widely spread by other algorithms as an efficient technique. The particularity of NSGA-II is to hierarchize the levels of ‘domination’, with a first Pareto front containing only the non-dominated solutions, a second Pareto front with the solutions dominated by one or two solutions, and finally, a third Pareto front with all the other solutions, those dominated by more than two solutions. For this last category, we compute ‘crowding distances’ for the solutions of this category, then we sort the set of results thus obtained. The ‘crowding distance’ is calculated criterion by criterion. For example, for the criterion represented on the horizontal axis, we will first determine the extreme solutions, which we will call ‘min’ and ‘max’ it being understood that ‘min’ will be the solution which will have the smallest value on the horizontal axis, and ‘max’ the solution which will have the largest value on the horizontal axis. Considering that we have $Q' < Q$ solutions on the Pareto front, we will then assign an index to each solution, with the index 1 for ‘min’. and the index $Q'$ for ‘max’. For each solution of index $q$ with $1 < q < Q'$ we will calculate the distance on the horizontal axis between the $(q-1)^{th}$ solution and the solution $(q+1)^{th}$ solution, and we will divide this distance by the distance on the horizontal axis between ‘max’ and ‘min’. Thus, it is this result that will be considered as the crowding distance of each solution.

The next step is to create a ‘descendant’. To do this, we first organize a selection tournament, i.e. we will randomly draw pairs of solutions in the set of solutions, and for each pair of solutions, we will determine which one can become ‘parent’, by comparing the hierarchical levels of the Pareto front. For example, if the first solution belongs to the second Pareto front and the second solution belongs to the third Pareto front, then the selection tournament will be won by the first solution because the second Pareto front contains solutions that are better than the third Pareto front. If two solutions belonging to the same Pareto front were drawn, then the solution with the smallest crowding distance would be selected as the ‘parent’. Then, for each pair of ‘parents’, we will generate a child who will have the values of the first ‘parent’ for some unknowns, and the values of the other “parent” for the other unknowns. This step is called ‘cross-over’. Finally, according to a percentage defined beforehand, the values of some unknowns of the child will be slightly modified. This last step is called the ‘mutation’. The whole process is repeated, as many times as there are iterations, and in the end, we obtain the set of ‘optimized’ solutions.

# 4 Chaotic ternary grey wolf algorithm

We propose a ternary version of GWO, which searches specifically ternary spaces with enhanced exploration abilities.

Our first motivation is that, in the considered application, the search space is associated to three values of epsilon. But this method could be applied to other situations and applications, involving for instance sensors with three possible states.

Our second motivation is to propose a method with enhanced exploration properties. Indeed metaheuristics with enhanced exploration properties are of great interest to cope with applications where the objective function exhibits an elevated number of local minima. We aim at achieving such enhanced exploration properties while proposing a ternary map which evolves across iterations, and inserting chaotic sequences in the update rules of the agents.

Our third motivation is to improve the diversity of the agents behavior. Indeed, GWO exhibits premature convergence due to poor diversity of the population of wolves. So we propose to divide the wolf pack into two groups: the first with enhanced ‘exploration’ abilities, and the second with ‘exploitation’ abilities.

The proposed chaotic ternary GWO is denoted by CTGWO. In Section 4.1, we derive the update rules which relies on specific transform functions depending on a parameter $a$. We wish to preserve the original philosophy of GWO: the number of leaders ruling the update of the agents is superior to 1, and the parameter $a$ permits to distinguish between an exploration phase at the beginning of the algorithm and an exploitation phase at the end. In Section 4.1, we just assume about parameter $a$ that it decreases from 2 to 0 across the iterations. Then in Section 4.2, we investigate a chaotic expression for $a$.

## 4.1 Ternary update rules based on dedicated transform maps

We propose here innovative update rules, dedicated to a ternary search space, performed with the help of an ad hoc transform function. Firstly propose a novel manner to compute the contribution of a leader, and the mean contribution of several leaders. Secondly, we propose an update rule depending on this mean contribution.

### 4.1.1 Contribution of a leader

We remind that in the continuous case, the contribution $y_l^i(n)$ of a leader $l$ is computed as in equation (8). In this case $y_l^i(n) = x_l^i - b^i b_l^i(n)$ depends on the product $b^i d_l^i(n)$ which decreases to 0 simultaneously with $a$, reaching 0 when $n = n_{max}$. So $y_l^i(n)$ is a real value which gets closer to $x_l^i$ across the iterations. In the following we propose to compute a contribution $y_l^i(n)$ which is a real value between 0 and 2. We propose, as expression of the contribution of leader $l$:

$$y_l^i(n) = (x_l^i - b^i d_l^i(n)) \mod 2, \quad (10)$$
and

$$((f_1(\mathbf{x}_{q1}) < f_1(\mathbf{x}_{q2})) \text{ or } (f_2(\mathbf{x}_{q1}) < f_2(\mathbf{x}_{q2}))).$$

It is up to the user to choose the best solution(s) to keep. Indeed, the user may very well want to favour one of the two objectives, or seek the best compromise between the two objectives.

Actually, in the last decades, different methods to determine the ‘domination’ principle have been proposed, and then, the ‘non-domination’ principle emerged, notably thanks to Srinivas and Deb [27] who proposed the NSGA algorithm [27] and then an improved version, called NSGA-II [28]. This fast sorting method by ‘non-domination’ has been widely spread by other algorithms as an efficient technique. The particularity of NSGA-II is to hierarchize the levels of ‘domination’, with a first Pareto front containing only the non-dominated solutions, a second Pareto front with the solutions dominated by one or two solutions, and finally, a third Pareto front with all the other solutions, those dominated by more than two solutions. For this last category, we compute ‘crowding distances’ for the solutions of this category, then we sort the set of results thus obtained. The ‘crowding distance’ is calculated criterion by criterion. For example, for the criterion represented on the horizontal axis, we will first determine the extreme solutions, which we will call ‘min’ and ‘max’ it being understood that ‘min’ will be the solution which will have the smallest value on the horizontal axis, and ‘max’ the solution which will have the largest value on the horizontal axis. Considering that we have $Q' < Q$ solutions on the Pareto front, we will then assign an index to each solution, with the index 1 for ‘min’. and the index $Q'$ for ‘max’. For each solution of index $q$ with $1 < q < Q'$ we will calculate the distance on the horizontal axis between the $(q-1)^{th}$ solution and the solution $(q+1)^{th}$ solution, and we will divide this distance by the distance on the horizontal axis between ‘max’ and ‘min’. Thus, it is this result that will be considered as the crowding distance of each solution.

The next step is to create a ‘descendant’. To do this, we first organize a selection tournament, i.e. we will randomly draw pairs of solutions in the set of solutions, and for each pair of solutions, we will determine which one can become ‘parent’, by comparing the hierarchical levels of the Pareto front. For example, if the first solution belongs to the second Pareto front and the second solution belongs to the third Pareto front, then the selection tournament will be won by the first solution because the second Pareto front contains solutions that are better than the third Pareto front. If two solutions belonging to the same Pareto front were drawn, then the solution with the smallest crowding distance would be selected as the ‘parent’. Then, for each pair of ‘parents’, we will generate a child who will have the values of the first ‘parent’ for some unknowns, and the values of the other “parent” for the other unknowns. This step is called ‘cross-over’. Finally, according to a percentage defined beforehand, the values of some unknowns of the child will be slightly modified. This last step is called the ‘mutation’. The whole process is repeated, as many times as there are iterations, and in the end, we obtain the set of ‘optimized’ solutions.

# 4 Chaotic ternary grey wolf algorithm

We propose a ternary version of GWO, which searches specifically ternary spaces with enhanced exploration abilities.

Our first motivation is that, in the considered application, the search space is associated to three values of epsilon. But this method could be applied to other situations and applications, involving for instance sensors with three possible states.

Our second motivation is to propose a method with enhanced exploration properties. Indeed metaheuristics with enhanced exploration properties are of great interest to cope with applications where the objective function exhibits an elevated number of local minima. We aim at achieving such enhanced exploration properties while proposing a ternary map which evolves across iterations, and inserting chaotic sequences in the update rules of the agents.

Our third motivation is to improve the diversity of the agents behavior. Indeed, GWO exhibits premature convergence due to poor diversity of the population of wolves. So we propose to divide the wolf pack into two groups: the first with enhanced ‘exploration’ abilities, and the second with ‘exploitation’ abilities.

The proposed chaotic ternary GWO is denoted by CTGWO. In Section 4.1, we derive the update rules which relies on specific transform functions depending on a parameter $a$. We wish to preserve the original philosophy of GWO: the number of leaders ruling the update of the agents is superior to 1, and the parameter $a$ permits to distinguish between an exploration phase at the beginning of the algorithm and an exploitation phase at the end. In Section 4.1, we just assume about parameter $a$ that it decreases from 2 to 0 across the iterations. Then in Section 4.2, we investigate a chaotic expression for $a$.

## 4.1 Ternary update rules based on dedicated transform maps

We propose here innovative update rules, dedicated to a ternary search space, performed with the help of an ad hoc transform function. Firstly propose a novel manner to compute the contribution of a leader, and the mean contribution of several leaders. Secondly, we propose an update rule depending on this mean contribution.

### 4.1.1 Contribution of a leader

We remind that in the continuous case, the contribution $y_l^i(n)$ of a leader $l$ is computed as in equation (8). In this case $y_l^i(n) = x_l^i - b^i b_l^i(n)$ depends on the product $b^i d_l^i(n)$ which decreases to 0 simultaneously with $a$, reaching 0 when $n = n_{max}$. So $y_l^i(n)$ is a real value which gets closer to $x_l^i$ across the iterations. In the following we propose to compute a contribution $y_l^i(n)$ which is a real value between 0 and 2. We propose, as expression of the contribution of leader $l$:

$$y_l^i(n) = (x_l^i - b^i d_l^i(n)) \mod 2, \quad (10)$$
where $b^i$ is defined as in equation (9) and $d_l^i(n)$ is defined as in equation (7). Based on the hypothesis that the maximum value of $a$ is 2, we deduce from equations (6) and (7) that the values of $b^i d_l^i(n)$ are between -8 and 8. So the values of $x_l^i - b^i d_l^i(n)$ can be out of the interval $[0:2]$.

The 'modulo' operator, denoted by $\mod$ is meant to enforce the contributions $y_l^i(n)$ to remain into the interval $[0:2]$. It is defined as follows: whatever the values $z_1 \in \mathbb{R}$ and $z_2 \in \mathbb{R}^*$:

$$z_1 \mod z_2 = \begin{cases} z_1 - z_2 \lfloor z_1/z_2 \rfloor & \text{if } z_1 \neq z_2 \\ z_2 & \text{if } z_1 = z_2, \text{ or } z_1 = 0 \end{cases} \quad (11)$$

where $|z_1|$ denotes the largest value in $\mathbb{Z}$ which is smaller than the scalar $z_1 \in \mathbb{R}$.

### 4.1.2 Weighted contribution of the leaders

We denote by $\overline{y^i(n)}$ the weighted contribution of four leaders $\alpha, \beta, \delta$, and $\rho$:

$$\overline{y^i(n)} = \begin{cases} \frac{1}{3}(y_{\alpha}^i(n) + y_{\beta}^i(n) + (1-a/2)y_{\delta}^i(n) + a/2y_{\rho}^i(n)) & \text{if } a > 1 \\ \frac{1}{3}(y_{\alpha}^i(n) + y_{\beta}^i(n) + y_{\delta}^i(n)) & \text{if } a \le 1. \end{cases} \quad (12)$$

In equation (12), leader $\rho$ is a wolf which is selected at random among the wolf pack.

We notice that $\overline{y^i(n)}$ gets closer to $\frac{1}{3}(y_{\alpha}^i(n) + y_{\beta}^i(n) + y_{\delta}^i(n))$ when the iteration index increases.

We can now derive a ternary rule which updates the position of any wolf in a ternary search space.

### 4.1.3 Ternary update rule

We propose an evolving map which enhances exploration at the beginning of the algorithm, and exploitation at the end of the algorithm. Compared to recent works about the binary GWO [34, 35], the new feature in the proposed transform map is of course the division of the map into three regions instead of 2, but also, the fact that the map evolves across iterations: we will introduce a term which is proportional to $a$ in the transform functions and which emphasizes exploration at the beginning of the algorithm and exploitation at the end of the algorithm.

The proposed novel process dedicated to ternary search spaces permits to select either value 0, 1 or 2. In dimension $i$ ($i = 1,...,P$) wolf $q$ is updated from iteration $n$ to iteration $n + 1$ as follows:

$$x_q^i(n+1) = \begin{cases} 0 & \text{if } r \ge \varphi^u(\overline{y^i(n)}, a) \\ 1 & \text{if } r < \varphi^u(\overline{y^i(n)}, a) \text{ and } r \ge \varphi^d(\overline{y^i(n)}, a) \\ 2 & \text{if } r < \varphi^d(\overline{y^i(n)}, a), \end{cases} \quad (13)$$

where the scalar $r$ is a random value between 0 and 1 and taken from a normal distribution. In equation (13) we introduce two functions, which are necessary to define the ternary map. These functions are denoted by $\varphi^u$ and $\varphi^d$:

$$\varphi^u : [0:2] \times \mathbb{R}_+ \rightarrow [0:1]; y \mapsto \varphi^u(y, a)$$

$$\varphi^d : [0:2] \times \mathbb{R}_+ \rightarrow [0:1]; y \mapsto \varphi^d(y, a).$$

Function $\varphi^u$ separates the uppermost part of the map from the rest of the map; and function $\varphi^d$ separates the lowermost part of the map from the rest of the map. The basic idea is that if the random number $r$ leads to the region in-between, the value 1 will be chosen as an updated value. If $r$ leads to the region which is above $\varphi^u$ (resp. below $\varphi^d$), the value 0 (resp. 2) will be selected. We will now detail the shape of the frontiers between regions 0, 1, and 2. We set a ternary map based on a 'power function'. The basic function we use is a power function applied to any scalar $y$ and depending on 5 parameters $c_1, c_2, c_3, c_4$, and $a$:

$$\varphi^p(y, c_1, c_2, c_3, c_4, a) = \left( \frac{(y-c_3)}{c_1} \right)^{c_2} + c_4 + \frac{a}{5} \left( \exp\left(-\frac{y^2}{2}\right) - \exp\left(-\frac{(y-2)^2}{2}\right) \right). \quad (14)$$

In equation (14), the first term depending on $c_1, c_2, c_3$ and $c_4$ gives the overall shape of the function. The second term depending on $a$ permits to get $\varphi^p(0, c_1, c_2, c_3, c_4, a) \ge 0$ and $\varphi^p(2, c_1, c_2, c_3, c_4, a) \le 1$. We use two versions of this function to define the ternary map. The first one, with $c_3 = 2$ and $c_4 = 1$; the second one, with $c_3 = 0$ and $c_4 = 0$:

$$\varphi^u(y, a) = \varphi^p(y, 2, 3, 2, 1, a) \quad (15)$$

$$\varphi^d(y, a) = \varphi^p(y, 2, 3, 0, 0, a). \quad (16)$$

The second term depending on $a$ (see Eq. (14)) permits to get values on $y = 0$ which are slightly larger than 0, and values on $y = 2$ which are slightly smaller than 1:

$$\varphi^P(0, c_1, c_2, c_3, c_4, a) = \frac{a}{5} \left(1 - \exp(-2)\right), \\ \text{and } \varphi^P(2, c_1, c_2, c_3, c_4, a) = 1 - \frac{a}{5} \left(1 - \exp(-2)\right)$$

in both equations (15) and (16).

So:

$$\varphi^u(0, a) = \varphi^d(0, a) = \frac{a}{5}(1 - \exp(-2)) \quad (17)$$

and

$$\varphi^u(2, a) = \varphi^d(2, a) = 1 - \frac{a}{5}(1 - \exp(-2)). \quad (18)$$
where $b^i$ is defined as in equation (9) and $d_l^i(n)$ is defined as in equation (7). Based on the hypothesis that the maximum value of $a$ is 2, we deduce from equations (6) and (7) that the values of $b^i d_l^i(n)$ are between -8 and 8. So the values of $x_l^i - b^i d_l^i(n)$ can be out of the interval $[0:2]$.

The ‘modulo’ operator, denoted by $\mod$ is meant to enforce the contributions $y_l^i(n)$ to remain into the interval $[0:2]$. It is defined as follows: whatever the values $z_1 \in \mathbb{R}$ and $z_2 \in \mathbb{R}^*$:

$$z_1 \mod z_2 = \begin{cases} z_1 - z_2 \lfloor z_1/z_2 \rfloor & \text{if } z_1 \neq z_2 \\ z_2 & \text{if } z_1 = z_2, \text{ or } z_1 = 0 \end{cases} \quad (11)$$

where $|z_1|$ denotes the largest value in $\mathbb{Z}$ which is smaller than the scalar $z_1 \in \mathbb{R}$.

### 4.1.2 Weighted contribution of the leaders

We denote by $\overline{y^i(n)}$ the weighted contribution of four leaders $\alpha, \beta, \delta$, and $\rho$:

$$\overline{y^i(n)} = \begin{cases} \frac{1}{3}(y_{\alpha}^i(n) + y_{\beta}^i(n) + (1-a/2)y_{\delta}^i(n) + a/2y_{\rho}^i(n)) & \text{if } a > 1 \\ \frac{1}{3}(y_{\alpha}^i(n) + y_{\beta}^i(n) + y_{\delta}^i(n)) & \text{if } a \le 1. \end{cases} \quad (12)$$

In equation (12), leader $\rho$ is a wolf which is selected at random among the wolf pack.

We notice that $\overline{y^i(n)}$ gets closer to $\frac{1}{3}(y_{\alpha}^i(n) + y_{\beta}^i(n) + y_{\delta}^i(n))$ when the iteration index increases.

We can now derive a ternary rule which updates the position of any wolf in a ternary search space.

### 4.1.3 Ternary update rule

We propose an evolving map which enhances exploration at the beginning of the algorithm, and exploitation at the end of the algorithm. Compared to recent works about the binary GWO [34, 35], the new feature in the proposed transform map is of course the division of the map into three regions instead of 2, but also, the fact that the map evolves across iterations: we will introduce a term which is proportional to $a$ in the transform functions and which emphasizes exploration at the beginning of the algorithm and exploitation at the end of the algorithm.

The proposed novel process dedicated to ternary search spaces permits to select either value 0, 1 or 2. In dimension $i$ ($i = 1,...,P$) wolf $q$ is updated from iteration $n$ to iteration $n + 1$ as follows:

$$x_q^i(n+1) = \begin{cases} 0 & \text{if } r \ge \varphi^u(\overline{y^i(n)}, a) \\ 1 & \text{if } r < \varphi^u(\overline{y^i(n)}, a) \text{ and } r \ge \varphi^d(\overline{y^i(n)}, a) \\ 2 & \text{if } r < \varphi^d(\overline{y^i(n)}, a), \end{cases} \quad (13)$$

where the scalar $r$ is a random value between 0 and 1 and taken from a normal distribution. In equation (13) we introduce two functions, which are necessary to define the ternary map. These functions are denoted by $\varphi^u$ and $\varphi^d$:

$$\varphi^u : [0:2] \times \mathbb{R}_+ \rightarrow [0:1]; y \mapsto \varphi^u(y, a)$$

$$\varphi^d : [0:2] \times \mathbb{R}_+ \rightarrow [0:1]; y \mapsto \varphi^d(y, a).$$

Function $\varphi^u$ separates the uppermost part of the map from the rest of the map; and function $\varphi^d$ separates the lowermost part of the map from the rest of the map. The basic idea is that if the random number $r$ leads to the region in-between, the value 1 will be chosen as an updated value. If $r$ leads to the region which is above $\varphi^u$ (resp. below $\varphi^d$), the value 0 (resp. 2) will be selected. We will now detail the shape of the frontiers between regions 0, 1, and 2. We set a ternary map based on a ‘power function’. The basic function we use is a power function applied to any scalar $y$ and depending on 5 parameters $c_1, c_2, c_3, c_4$, and $a$:

$$\varphi^p(y, c_1, c_2, c_3, c_4, a) = \left( \frac{(y-c_3)}{c_1} \right)^{c_2} + c_4 + \frac{a}{5} \left( \exp\left(-\frac{y^2}{2}\right) - \exp\left(-\frac{(y-2)^2}{2}\right) \right). \quad (14)$$

In equation (14), the first term depending on $c_1, c_2, c_3$ and $c_4$ gives the overall shape of the function. The second term depending on $a$ permits to get $\varphi^p(0, c_1, c_2, c_3, c_4, a) \ge 0$ and $\varphi^p(2, c_1, c_2, c_3, c_4, a) \le 1$. We use two versions of this function to define the ternary map. The first one, with $c_3 = 2$ and $c_4 = 1$; the second one, with $c_3 = 0$ and $c_4 = 0$:

$$\varphi^u(y, a) = \varphi^p(y, 2, 3, 2, 1, a) \quad (15)$$

$$\varphi^d(y, a) = \varphi^p(y, 2, 3, 0, 0, a). \quad (16)$$

The second term depending on $a$ (see Eq. (14)) permits to get values on $y = 0$ which are slightly larger than 0, and values on $y = 2$ which are slightly smaller than 1:

$$\varphi^P(0, c_1, c_2, c_3, c_4, a) = \frac{a}{5} \left(1 - \exp(-2)\right), \\ \text{and } \varphi^P(2, c_1, c_2, c_3, c_4, a) = 1 - \frac{a}{5} \left(1 - \exp(-2)\right)$$

in both equations (15) and (16).

So:

$$\varphi^u(0, a) = \varphi^d(0, a) = \frac{a}{5}(1 - \exp(-2)) \quad (17)$$

and

$$\varphi^u(2, a) = \varphi^d(2, a) = 1 - \frac{a}{5}(1 - \exp(-2)). \quad (18)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Ternary ‘Power’ selection map with different values of parameter a. (a) a = 2; (b) a = 1; (c) a = 0.</figcaption></figure>

The functions $\varphi^u(y, a)$ and $\varphi^d(y, a)$ defined in equations (15) and (16) can then be used in equation (13) to get the ‘Power’ transform map.

### 4.1.4 Representation of the ternary transform map and interpretation

Figure 2 shows the ‘Power’ ternary map. In each case the uppermost region maps for 0, the central region maps for 1, and the lowermost region maps for 2. It can be noticed that the shape of the functions $\varphi^u$ and $\varphi^d$ is consistent with the derivations in equations (17) and (18).

We can check that:

- The term which is proportional to $a$ in equation (14) yields a map which evolves across the iterations; this is an important difference compared to the binary map presented in [34];

- For a value of parameter $a$ which is 2, either the value 0 or 2 may be selected with an elevated probability when all leader contributions are equal to 0, or 2;

- For a value of parameter $a$ which is 0, the value 0 (resp. 2) is selected with probability 1 when $\overline{y^i(n)} = 0$ (resp. $\overline{y^i(n)} = 2$);

- Whatever the value of $a$, either the values 0, 1, or 2 may be selected when $\overline{y^i(n)} = 1$.

Indeed, as defined in equation (10) the values of $y_1^i(n)$ are real and between 0 and 2, whatever $i$ and 1. In Section 4.2 we embed chaotic sequences in the expression of $a$ to improve again the exploration abilities of our algorithm.

## 4.2 Embedding chaotic sequences in the ternary grey wolf optimizer

### 4.2.1 Chaotic expression of parameter ‘a’ for improved exploration abilities

We modify the expression of parameter $a$ with respect to other versions of GWO, and propose:

$$a = 2 \left( 1 - \left( \frac{n}{n_{\max}} \right)^{(n_q(1+\Gamma(c_q^1, n)))} \right). \quad (19)$$

In equation (19), we notice that, for the first time in this paper, the expression of $a$ depends on the agent index $q$. Firstly, we reposition the worst agents closer to the three leaders, with a value of $\eta$ which depends on the score of the agent: for the worst half of the agents (associated with the largest scores), $\eta_q = \frac{2}{3}$; for the best half of the agents (associated with the smallest scores), $\eta_q = \frac{3}{2}$. Secondly, inserting a chaotic sequence $\Gamma(c_q^1, n)$ enhances the variability of the behavior of the agents. $c_q^1$ is the initial value of the chaotic sequence, which is different for every agent $q$. The principles of the proposed method is that the $a$ parameter of the GWO which rules the displacement of the agents is perturbed through the value of a chaotic sequence.

Choosing a different value of $c_q^1$ or each agent permits to emphasize diversity in the displacement of the agents. Meanwhile, we choose a sequence with one attractor to ensure that $\Gamma(c_q^1, n_{\max})$ is the same whatever $q$.

### 4.2.2 Construction of the chaotic sequences

To privilege exploration abilities, the behaviors of the agents should differ one from the other. To privilege exploitation abilities at the end of the algorithm, the agents behavior should get closer to each other while the iteration index increases. So we set the following constraints on the chaotic sequences:

- The values of $\Gamma(c_q^1, n)$ are positive, in the interval $[0:1]$;

- The last value should be $\Gamma(c_q^1, n_{\max}) = 0$ whatever the initial value $\Gamma(c_q^1, 1)$. In this way, at the last iteration $n = n_{\max}$, the behavior of all agents is the same.

To fulfill easily those constraints, we have to choose a sequence with one known attractor. We base our sequence $\Gamma$ on a logistic sequence, denoted by $c(n)$. Given an initial term $c(1)$, each subsequent term (for $n = 2,...,n_{\max}$) is defined as:

$$c(n + 1) = \kappa c(n)(1 - c(n)), \quad (20)$$

where $\kappa \in \mathbb{R}_+^*$. The number of attractors for this sequence depends on the value of $\kappa$. Figure 3 shows chaotic sequences with one (Fig. 3a), two (Fig. 3b), or an infinite number of attractors (Fig. 3c). We chose $\kappa = 2.8$ in Figure 3a, $\kappa = 3.2$ in Figure 3b, $\kappa = 3.99$ in Figure 3c. In Figures 3a–3c, each
<figure><img src="image_3.png" /><figcaption>Fig. 2. Ternary ‘Power’ selection map with different values of parameter a. (a) a = 2; (b) a = 1; (c) a = 0.</figcaption></figure>

The functions $\varphi^u(y, a)$ and $\varphi^d(y, a)$ defined in equations (15) and (16) can then be used in equation (13) to get the ‘Power’ transform map.

### 4.1.4 Representation of the ternary transform map and interpretation

Figure 2 shows the ‘Power’ ternary map. In each case the uppermost region maps for 0, the central region maps for 1, and the lowermost region maps for 2. It can be noticed that the shape of the functions $\varphi^u$ and $\varphi^d$ is consistent with the derivations in equations (17) and (18).

We can check that:

- The term which is proportional to $a$ in equation (14) yields a map which evolves across the iterations; this is an important difference compared to the binary map presented in [34];

- For a value of parameter $a$ which is 2, either the value 0 or 2 may be selected with an elevated probability when all leader contributions are equal to 0, or 2;

- For a value of parameter $a$ which is 0, the value 0 (resp. 2) is selected with probability 1 when $\overline{y^i(n)} = 0$ (resp. $\overline{y^i(n)} = 2$);

- Whatever the value of $a$, either the values 0, 1, or 2 may be selected when $\overline{y^i(n)} = 1$.

Indeed, as defined in equation (10) the values of $y_1^i(n)$ are real and between 0 and 2, whatever $i$ and 1. In Section 4.2 we embed chaotic sequences in the expression of $a$ to improve again the exploration abilities of our algorithm.

## 4.2 Embedding chaotic sequences in the ternary grey wolf optimizer

### 4.2.1 Chaotic expression of parameter ‘a’ for improved exploration abilities

We modify the expression of parameter $a$ with respect to other versions of GWO, and propose:

$$a = 2 \left( 1 - \left( \frac{n}{n_{\max}} \right)^{(n_q(1+\Gamma(c_q^1,n)))} \right). \quad (19)$$

In equation (19), we notice that, for the first time in this paper, the expression of $a$ depends on the agent index $q$. Firstly, we reposition the worst agents closer to the three leaders, with a value of $\eta$ which depends on the score of the agent: for the worst half of the agents (associated with the largest scores), $\eta_q = \frac{2}{3}$; for the best half of the agents (associated with the smallest scores), $\eta_q = \frac{3}{2}$. Secondly, inserting a chaotic sequence $\Gamma(c_q^1, n)$ enhances the variability of the behavior of the agents. $c_q^1$ is the initial value of the chaotic sequence, which is different for every agent $q$. The principles of the proposed method is that the $a$ parameter of the GWO which rules the displacement of the agents is perturbed through the value of a chaotic sequence.

Choosing a different value of $c_q^1$ or each agent permits to emphasize diversity in the displacement of the agents. Meanwhile, we choose a sequence with one attractor to ensure that $\Gamma(c_q^1, n_{\max})$ is the same whatever $q$.

### 4.2.2 Construction of the chaotic sequences

To privilege exploration abilities, the behaviors of the agents should differ one from the other. To privilege exploitation abilities at the end of the algorithm, the agents behavior should get closer to each other while the iteration index increases. So we set the following constraints on the chaotic sequences:

- The values of $\Gamma(c_q^1, n)$ are positive, in the interval $[0:1]$;

- The last value should be $\Gamma(c_q^1, n_{\max}) = 0$ whatever the initial value $\Gamma(c_q^1, 1)$. In this way, at the last iteration $n = n_{\max}$, the behavior of all agents is the same.

To fulfill easily those constraints, we have to choose a sequence with one known attractor. We base our sequence $\Gamma$ on a logistic sequence, denoted by $c(n)$. Given an initial term $c(1)$, each subsequent term (for $n = 2,...,n_{\max}$) is defined as:

$$c(n + 1) = \kappa c(n)(1 - c(n)), \quad (20)$$

where $\kappa \in \mathbb{R}_+^*$. The number of attractors for this sequence depends on the value of $\kappa$. Figure 3 shows chaotic sequences with one (Fig. 3a), two (Fig. 3b), or an infinite number of attractors (Fig. 3c). We chose $\kappa = 2.8$ in Figure 3a, $\kappa = 3.2$ in Figure 3b, $\kappa = 3.99$ in Figure 3c. In Figures 3a–3c, each
<figure><img src="image_4.png" /><figcaption>Fig. 3. Chaotic sequences with various values of chaos parameter $κ$. (a) $κ = 2.80$, one attractor; (b) $κ = 3.20$, two attractors; (c) $κ = 3.99$, an infinite number of attractors.</figcaption></figure>

plot with a given color corresponds to a different value for $c(1)$. There are ten chaotic sequences in each case. We choose a sequence such as in Figure 3a, where $c(n_{max})$ bears the same value, equal to 0.65 approximately, whatever the sequence. For this we set with $κ = 2.8$.

For any agent q and for $n = 1,...,n_{max}$, we define the ad hoc chaotic sequence $Γ(c_q^1, n)$ as follows, based on the logistic sequence of equation (20):

$$Γ(c_q^1, n) = 0.1(c(n) - c(n_{max})). \quad (21)$$

Hence, for any agent q$Γ(c_q^1, n_{max}) = 0$. We notice that the initial term $c_q^1$ is defined as follows:

$$c_q^1 = 0.1(c(1) - c(n_{max})). \quad (22)$$

We set the initial term $c(1)$ of the logistic sequence as a random value between 0 and 1, taken from a normal distribution. As $c(1)$ is a random value, $c_q^1$ is also random and a different sequence is generated for each agent.

Figure 4 shows the sequences of values for a for all agents. We clearly distinguish two families of agents: the first family with $η_q = \frac{3}{2}$ and rather elevated values of a, and the second family with $η_q = \frac{2}{3}$ and smaller values of a, which tend more rapidly towards 0.

# 5 Results: mono-objective and bi-objective approaches

## 5.1 Experimental conditions and metrics

In this section, the test environment is a server running Linux, equipped with 4 Intel(R) Xeon(R) CPU X7560 @ 2.27GHz (64 processors, Hyper-threading activated) and 1000 GB RAM. The software is Python.

We consider cloaks with $P = 329$ voxels (of triangular shape) and the expected outcomes of the algorithms are vectors $x$ containing permittivity values $K^1, K^2, ..., K^i, ..., K^P$.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Chaotic values of parameter a for all agents, and all iterations $n = 1,...,50$.</figcaption></figure>

In Section 5.2 we display the results obtained with a monoobjective approach. The criterion which is minimized is $f(x):\mathbb{R}^P \mapsto \mathbb{R}_+: f(x) = \frac{1}{2}(f_1(x) + f_2(x))$. We remind (see Sect. 2) that $f_1$ stands for protection, and $f_2$ stands for invisibility.

We remind that topological optimization limits the number of possible materials to two. Instead, and for the first time in this paper, we adapt our new variant CTGWO of the GWO which searches solutions among three materials. We compare the proposed CTGWO with the adaptive mixed GWO in discrete mode [31] (denoted by amixedGWO), and the vanilla continuous GWO [29]. The computational time for one trial of either $f_1$ or $f_2$ is 0.99 sec by trial. For GWO and amixedGWO, we use the expression of a presented in equation (9). For CTGWO, we use the expression of a presented in equation (19). The search space for the values of the permittivity $ε$ is $\{7,10,12\}$. CTGWO
<figure><img src="image_4.png" /><figcaption>Fig. 3. Chaotic sequences with various values of chaos parameter $κ$. (a) $κ = 2.80$, one attractor; (b) $κ = 3.20$, two attractors; (c) $κ = 3.99$, an infinite number of attractors.</figcaption></figure>

plot with a given color corresponds to a different value for $c(1)$. There are ten chaotic sequences in each case. We choose a sequence such as in Figure 3a, where $c(n_{max})$ bears the same value, equal to 0.65 approximately, whatever the sequence. For this we set with $κ = 2.8$.

For any agent q and for $n = 1,...,n_{max}$, we define the ad hoc chaotic sequence $Γ(c_q^1, n)$ as follows, based on the logistic sequence of equation (20):

$$Γ(c_q^1, n) = 0.1(c(n) - c(n_{max})). \quad (21)$$

Hence, for any agent q$Γ(c_q^1, n_{max}) = 0$. We notice that the initial term $c_q^1$ is defined as follows:

$$c_q^1 = 0.1(c(1) - c(n_{max})). \quad (22)$$

We set the initial term $c(1)$ of the logistic sequence as a random value between 0 and 1, taken from a normal distribution. As $c(1)$ is a random value, $c_q^1$ is also random and a different sequence is generated for each agent.

Figure 4 shows the sequences of values for a for all agents. We clearly distinguish two families of agents: the first family with $η_q = \frac{3}{2}$ and rather elevated values of a, and the second family with $η_q = \frac{2}{3}$ and smaller values of a, which tend more rapidly towards 0.

# 5 Results: mono-objective and bi-objective approaches

## 5.1 Experimental conditions and metrics

In this section, the test environment is a server running Linux, equipped with 4 Intel(R) Xeon(R) CPU X7560 @ 2.27GHz (64 processors, Hyper-threading activated) and 1000 GB RAM. The software is Python.

We consider cloaks with $P = 329$ voxels (of triangular shape) and the expected outcomes of the algorithms are vectors $x$ containing permittivity values $K^1, K^2, ..., K^i, ..., K^P$.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Chaotic values of parameter a for all agents, and all iterations $n = 1,...,50$.</figcaption></figure>

In Section 5.2 we display the results obtained with a monoobjective approach. The criterion which is minimized is $f(x):\mathbb{R}^P \mapsto \mathbb{R}_+: f(x) = \frac{1}{2}(f_1(x) + f_2(x))$. We remind (see Sect. 2) that $f_1$ stands for protection, and $f_2$ stands for invisibility.

We remind that topological optimization limits the number of possible materials to two. Instead, and for the first time in this paper, we adapt our new variant CTGWO of the GWO which searches solutions among three materials. We compare the proposed CTGWO with the adaptive mixed GWO in discrete mode [31] (denoted by amixedGWO), and the vanilla continuous GWO [29]. The computational time for one trial of either $f_1$ or $f_2$ is 0.99 sec by trial. For GWO and amixedGWO, we use the expression of a presented in equation (9). For CTGWO, we use the expression of a presented in equation (19). The search space for the values of the permittivity $ε$ is $\{7,10,12\}$. CTGWO
<table><caption>**Table 1.** Cloak design experimental conditions: mono-objective approach (GWO, amixedGWO, CTGWO); and bi-objective approach (NSGA-II).</caption><thead><tr><th>Approach</th><th>Search space</th><th>n<sub>max</sub></th><th>Q</th><th>Link</th></tr></thead><tbody><tr><td>Mono-objective</td><td>Ternary</td><td>250</td><td>100</td><td>5.2</td></tr><tr><td rowspan="4">Bi-objective</td><td>Ternary</td><td>250</td><td>100</td><td>5.3.1</td></tr><tr><td>Ternary</td><td>1000</td><td>200</td><td>5.3.2</td></tr><tr><td>Continuous</td><td>250</td><td>100</td><td>5.3.3</td></tr><tr><td>Continuous</td><td>1000</td><td>200</td><td>5.3.4</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Convergence curve of GWO, amixedGWO, and CTGWO algorithms.</figcaption></figure>

and amixedGWO access these values *via* indices retrieved from the discrete search space {1,1,2}. GWO access these values *via* rounded indices retrieved from the continuous search space [0:2]. We run the three algorithms with $Q = 100$ agents and $n_{max} = 250$ iterations, that is, $25 \times 10^3$ trials of the objective function. The agents are initialized with random integers between 0 and 2.

In Section 5.3, we display the results obtained with a biobjective approach. The couple of criteria which are jointly minimized are $f_1(\mathbf{x}):\mathbb{R}^P \mapsto \mathbb{R}_+$, and $f_2(\mathbf{x}):\mathbb{R}^P \mapsto \mathbb{R}_+$. We display the results of GWO, amixed GWO, and the proposed CTGWO; as well as the results obtained by NSGA-II in four situations. These five experimental conditions are summarized in Table 1. The ternary search space mentioned in Table 1 is {7,10,12}: NSGA-II accesses these values *via* rounded indices retrieved from the continuous search space [0:2]. The continuous search space mentioned in Table 1 is [7:12]. This last situation is prospective in the sense that we assume that we afford any material with any permittivity value between 7 and 12.

## 5.2 Mono-objective approach

In this section, we display the results obtained by CTGWO, amixedGWO, and GWO: the convergence curves in Figure 5; the scores reached by each method, and corresponding

protection and invisibility performances in Table 2; the cloak design and wave propagation field in Figure 6.

In Table 2, we display, for GWO, amixedGWO, and CTGWO, the score of wolf *a* at the last iteration $f(\mathbf{x}_\alpha(n_{max}))$. For instance, the solution provided by CTGWO is $\mathbf{x}_\alpha(n_{max})=[7,12,12,...,10]^T$. We notice that, in CTGWO, the power map is such that there is still a small probability that the number of the area in which we are going to place ourselves does not correspond to the value obtained, used as a position on the horizontal axis. The interest of such a mismatch is to avoid our algorithm to be locked in a local minimum. The chaotic sequences also play a role in the good behavior of CTGWO shown by the convergence curve in Figure 5. Figure 6f indicates that the protection abilities of the cloak designed by CTGWO are clearly better than for GWO and amixedGWO. Values in Table 2 confirm this impression: the protection criterion is two times smaller, while the mean criterion is also significantly smaller.

## 5.3 Bi-objective approach

In a second step, we used a bi-objective genetic algorithm as an optimization method. Thus, we no longer work on convergence curves, but on Pareto fronts, to balance the two minimized criteria.

The parameter values for NSGA-II are as follows: mutation and crossover probabilities are $p_m = 1/P \approx 0.003$ and $p_c = 0.9$; mutation and crossover distribution indices are $\eta_m = 20$ and $\eta_c = 15$.

### 5.3.1 Bi-objective ternary optimization with 250 iterations, 100 research agents

In this section, we display the results obtained by NSGA-II: the Pareto front in Figure 7; the protection and invisibility performances in Table 3; the cloak design and wave propagation field in Figure 8.

### 5.3.2 Bi-objective ternary optimization with 1000 iterations, 200 research agents

In this subsection, we display the results obtained by NSGA-II: the Pareto front in Figure 9; the protection and invisibility performances in Table 4; the cloak design and wave propagation field in Figure 10.

### 5.3.3 Bi-objective continuous optimization with 250 iterations, 100 research agents

In this section, we display the results obtained by NSGA-II: the Pareto front in Figure 11; the protection and invisibility performances in Table 5; the cloak design and wave propagation field in Figure 12.

### 5.3.4 Bi-objective continuous optimization with 1000 iterations, 200 research agents

In this section, we display the results obtained by NSGA-II: the Pareto front in Figure 13; the protection and invisibility performances in Table 6; the cloak design and wave propagation field in Figure 14.
<table><caption>**Table 1.** Cloak design experimental conditions: mono-objective approach (GWO, amixedGWO, CTGWO); and bi-objective approach (NSGA-II).</caption><thead><tr><th>Approach</th><th>Search space</th><th>n<sub>max</sub></th><th>Q</th><th>Link</th></tr></thead><tbody><tr><td>Mono-objective</td><td>Ternary</td><td>250</td><td>100</td><td>5.2</td></tr><tr><td rowspan="4">Bi-objective</td><td>Ternary</td><td>250</td><td>100</td><td>5.3.1</td></tr><tr><td>Ternary</td><td>1000</td><td>200</td><td>5.3.2</td></tr><tr><td>Continuous</td><td>250</td><td>100</td><td>5.3.3</td></tr><tr><td>Continuous</td><td>1000</td><td>200</td><td>5.3.4</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Convergence curve of GWO, amixedGWO, and CTGWO algorithms.</figcaption></figure>

and amixedGWO access these values *via* indices retrieved from the discrete search space {1,1,2}. GWO access these values *via* rounded indices retrieved from the continuous search space [0:2]. We run the three algorithms with $Q = 100$ agents and $n_{max} = 250$ iterations, that is, $25 \times 10^3$ trials of the objective function. The agents are initialized with random integers between 0 and 2.

In Section 5.3, we display the results obtained with a biobjective approach. The couple of criteria which are jointly minimized are $f_1(\mathbf{x}):\mathbb{R}^P \mapsto \mathbb{R}_+$, and $f_2(\mathbf{x}):\mathbb{R}^P \mapsto \mathbb{R}_+$. We display the results of GWO, amixed GWO, and the proposed CTGWO; as well as the results obtained by NSGA-II in four situations. These five experimental conditions are summarized in Table 1. The ternary search space mentioned in Table 1 is {7,10,12}: NSGA-II accesses these values *via* rounded indices retrieved from the continuous search space [0:2]. The continuous search space mentioned in Table 1 is [7:12]. This last situation is prospective in the sense that we assume that we afford any material with any permittivity value between 7 and 12.

## 5.2 Mono-objective approach

In this section, we display the results obtained by CTGWO, amixedGWO, and GWO: the convergence curves in Figure 5; the scores reached by each method, and corresponding

protection and invisibility performances in Table 2; the cloak design and wave propagation field in Figure 6.

In Table 2, we display, for GWO, amixedGWO, and CTGWO, the score of wolf *a* at the last iteration $f(\mathbf{x}_\alpha(n_{max}))$. For instance, the solution provided by CTGWO is $\mathbf{x}_\alpha(n_{max})=[7,12,12,...,10]^T$. We notice that, in CTGWO, the power map is such that there is still a small probability that the number of the area in which we are going to place ourselves does not correspond to the value obtained, used as a position on the horizontal axis. The interest of such a mismatch is to avoid our algorithm to be locked in a local minimum. The chaotic sequences also play a role in the good behavior of CTGWO shown by the convergence curve in Figure 5. Figure 6f indicates that the protection abilities of the cloak designed by CTGWO are clearly better than for GWO and amixedGWO. Values in Table 2 confirm this impression: the protection criterion is two times smaller, while the mean criterion is also significantly smaller.

## 5.3 Bi-objective approach

In a second step, we used a bi-objective genetic algorithm as an optimization method. Thus, we no longer work on convergence curves, but on Pareto fronts, to balance the two minimized criteria.

The parameter values for NSGA-II are as follows: mutation and crossover probabilities are $p_m = 1/P \approx 0.003$ and $p_c = 0.9$; mutation and crossover distribution indices are $\eta_m = 20$ and $\eta_c = 15$.

### 5.3.1 Bi-objective ternary optimization with 250 iterations, 100 research agents

In this section, we display the results obtained by NSGA-II: the Pareto front in Figure 7; the protection and invisibility performances in Table 3; the cloak design and wave propagation field in Figure 8.

### 5.3.2 Bi-objective ternary optimization with 1000 iterations, 200 research agents

In this subsection, we display the results obtained by NSGA-II: the Pareto front in Figure 9; the protection and invisibility performances in Table 4; the cloak design and wave propagation field in Figure 10.

### 5.3.3 Bi-objective continuous optimization with 250 iterations, 100 research agents

In this section, we display the results obtained by NSGA-II: the Pareto front in Figure 11; the protection and invisibility performances in Table 5; the cloak design and wave propagation field in Figure 12.

### 5.3.4 Bi-objective continuous optimization with 1000 iterations, 200 research agents

In this section, we display the results obtained by NSGA-II: the Pareto front in Figure 13; the protection and invisibility performances in Table 6; the cloak design and wave propagation field in Figure 14.
<table><caption>**Table 2.** Comparison of GWO, amixedGWO and CTGWO mono-objective methods in protection ($f_1$) and invisibility ($f_2$) for a rectangular and half-rectangular cloak. $n_{max}= 250$, $Q =100$.</caption><thead><tr><th></th><th>Protection f<sub>1</sub>(x<sub>α</sub>(n<sub>max</sub>))</th><th>Invisibility f<sub>2</sub>(x<sub>α</sub>(n<sub>max</sub>))</th><th>Score f(x<sub>α</sub>(n<sub>max</sub>))</th><th>Links to figures</th></tr></thead><tbody><tr><td>GWO</td><td>4.234 × 10<sup>-3</sup></td><td>1.740 × 10<sup>-3</sup></td><td>2.987 × 10<sup>-3</sup></td><td>Figure 6d</td></tr><tr><td>AmixedGWO</td><td>3.481 × 10<sup>-4</sup></td><td>9.572 × 10<sup>-4</sup></td><td>6.527 × 10<sup>-4</sup></td><td>Figure 6e</td></tr><tr><td>CTGWO</td><td>1.561 × 10<sup>-4</sup></td><td>8.696 × 10<sup>-4</sup></td><td>5.128 × 10<sup>-4</sup></td><td>Figure 6f</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Column 1: in (a) GWO optimized cloak, in (b) discrete GWO optimized cloak, and in (c), CTGWO optimized cloak; column 2: in (d) GWO result, in (e) Discrete GWO result, and in (f), CTGWO result.</figcaption></figure>
<table><caption>**Table 2.** Comparison of GWO, amixedGWO and CTGWO mono-objective methods in protection ($f_1$) and invisibility ($f_2$) for a rectangular and half-rectangular cloak. $n_{max}$= 250, $Q$ =100.</caption><thead><tr><th></th><th>Protection f<sub>1</sub>(x<sub>α</sub>(n<sub>max</sub>))</th><th>Invisibility f<sub>2</sub>(x<sub>α</sub>(n<sub>max</sub>))</th><th>Score f(x<sub>α</sub>(n<sub>max</sub>))</th><th>Links to figures</th></tr></thead><tbody><tr><td>GWO</td><td>4.234 × 10<sup>-3</sup></td><td>1.740 × 10<sup>-3</sup></td><td>2.987 × 10<sup>-3</sup></td><td>Figure 6d</td></tr><tr><td>AmixedGWO</td><td>3.481 × 10<sup>-4</sup></td><td>9.572 × 10<sup>-4</sup></td><td>6.527 × 10<sup>-4</sup></td><td>Figure 6e</td></tr><tr><td>CTGWO</td><td>1.561 × 10<sup>-4</sup></td><td>8.696 × 10<sup>-4</sup></td><td>5.128 × 10<sup>-4</sup></td><td>Figure 6f</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Column 1: in (a) GWO optimized cloak, in (b) discrete GWO optimized cloak, and in (c), CTGWO optimized cloak; column 2: in (d) GWO result, in (e) Discrete GWO result, and in (f), CTGWO result.</figcaption></figure>
<table><caption>**Table 3.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}=250, Q=100$, see Figure 7 for the Pareto front and Figure 8 for 2D plots of cloak design and wave fields.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>1.690 × 10<sup>-4</sup></td><td>8.922 × 10<sup>-4</sup></td><td>5.306 × 10<sup>-4</sup></td><td>Figure 8f</td></tr><tr><td>Best invisibility</td><td>1.077 × 10<sup>-3</sup></td><td>5.363 × 10<sup>-4</sup></td><td>8.064 × 10<sup>-4</sup></td><td>Figure 8d</td></tr><tr><td>Best compromise</td><td>4.643 × 10<sup>-4</sup></td><td>6.527 × 10<sup>-4</sup></td><td>5.585 × 10<sup>-4</sup></td><td>Figure 8e</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Left: Pareto front for NSGA-II Table with 250 iterations, 100 research agents and four selected solutions marked in red. Right: Magnified view of the four selected solutions, three of which appear in Table 3.</figcaption></figure>

## 5.4 Spectral tolerance and comparison with random cloaks

To further verify the robustness of our approach to a change in frequency, we compute the invisibility and protection criteria for several frequency values around the central freespace wavelength $\lambda_0$. We compare the results obtained on the optimized cloaks to the results obtained with random cloaks obtained by filling each voxel of the design space by an arbitrary integer (resp. real) value in {7,10,12} (resp. [7:12]). For the central freespace wavelength $\lambda_0$: for both $C_1$ (protection) and $C_2$ (invisibility) these criteria are 20 times smaller for the optimized cloaks: 1.5 × 10⁻⁴ vs 2.2 × 10⁻³ for $C_1$; and 6 × 10⁻⁵ vs 1.2 × 10⁻³ for $C_2$.

We infer from these results that it is worth optimizing the structure of the cloaks: this good behavior at the target frequency is obtained at the expense of slightly worse performances at other frequencies, in particular for the invisibility which appears to be quite resonant (optimized cloaks lead to a better invisibility than any random realizations for $\lambda/\lambda_0 \in [0.95:1.03]$), as opposed to the protection which turns out to be more broadband (the cloak optimized with C-NSGA-II+ leads to a better protection than any random realizations for $\lambda/\lambda_0 \in [0.82:1.22]$). We note that the optimization method operating over a continuous space leads to a more broadband response (see red and purple curves in Fig. 15).

## 6 Discussion

We distinguish two situations: the ternary case, where we afford three possible values for epsilon, and the continuous case, where we afford any real value between 7 and 12.

In the ternary case, when $n_{max}=250$ and $Q=100$, the mono-objective approach implemented with the proposed CTGWO yields the best results in terms of protection, that is, $1.56 \times 10^{-4}$ and score, that is, $f_1/f_2 = 5.128 \times 10^{-4}$. Figure 6f confirms this with a good protection behavior obtained with CTGWO. So in these conditions the proposed CTGWO algorithm yields the best trade-off between protection and invisibility.

NSGA-II yields the best result in terms of invisibility, that is, $5.363 \times 10^{-4}$, at the expense of a mean value $f_1/f_2 = 8.063 \times 10^{-4}$ which is larger than the score obtained by CTGWO. However the interest of NSGA-II is to enable an end-user to choose to favor one criterion (for instance invisibility) rather than the other. When $n_{max} = 1000$ and $Q=200$ affording therefore 8 times more trials of each objective function, the best protection reaches $1.233 \times 10^{-4}$, with an associated mean value $f_1/f_2 = 4.021 \times 10^{-4}$. In Figures 8 and 10, we can check the coherence of the shape of the wavefront and the score values obtained. A less realistic manner to improve the results is to enable the search for any real permittivity value between 7 and 12.
<table><caption>**Table 3.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}=250, Q=100$, see Figure 7 for the Pareto front and Figure 8 for 2D plots of cloak design and wave fields.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>1.690 × 10<sup>-4</sup></td><td>8.922 × 10<sup>-4</sup></td><td>5.306 × 10<sup>-4</sup></td><td>Figure 8f</td></tr><tr><td>Best invisibility</td><td>1.077 × 10<sup>-3</sup></td><td>5.363 × 10<sup>-4</sup></td><td>8.064 × 10<sup>-4</sup></td><td>Figure 8d</td></tr><tr><td>Best compromise</td><td>4.643 × 10<sup>-4</sup></td><td>6.527 × 10<sup>-4</sup></td><td>5.585 × 10<sup>-4</sup></td><td>Figure 8e</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Left: Pareto front for NSGA-II Table with 250 iterations, 100 research agents and four selected solutions marked in red. Right: Magnified view of the four selected solutions, three of which appear in Table 3.</figcaption></figure>

## 5.4 Spectral tolerance and comparison with random cloaks

To further verify the robustness of our approach to a change in frequency, we compute the invisibility and protection criteria for several frequency values around the central freespace wavelength $\lambda_0$. We compare the results obtained on the optimized cloaks to the results obtained with random cloaks obtained by filling each voxel of the design space by an arbitrary integer (resp. real) value in {7,10,12} (resp. [7:12]). For the central freespace wavelength $\lambda_0$: for both $C_1$ (protection) and $C_2$ (invisibility) these criteria are 20 times smaller for the optimized cloaks: 1.5 × 10⁻⁴ vs 2.2 × 10⁻³ for $C_1$; and 6 × 10⁻⁵ vs 1.2 × 10⁻³ for $C_2$.

We infer from these results that it is worth optimizing the structure of the cloaks: this good behavior at the target frequency is obtained at the expense of slightly worse performances at other frequencies, in particular for the invisibility which appears to be quite resonant (optimized cloaks lead to a better invisibility than any random realizations for $\lambda/\lambda_0 \in [0.95:1.03]$), as opposed to the protection which turns out to be more broadband (the cloak optimized with C-NSGA-II+ leads to a better protection than any random realizations for $\lambda/\lambda_0 \in [0.82:1.22]$). We note that the optimization method operating over a continuous space leads to a more broadband response (see red and purple curves in Fig. 15).

## 6 Discussion

We distinguish two situations: the ternary case, where we afford three possible values for epsilon, and the continuous case, where we afford any real value between 7 and 12.

In the ternary case, when $n_{max}=250$ and $Q=100$, the mono-objective approach implemented with the proposed CTGWO yields the best results in terms of protection, that is, $1.56 \times 10^{-4}$ and score, that is, $f_1/f_2 = 5.128 \times 10^{-4}$. Figure 6f confirms this with a good protection behavior obtained with CTGWO. So in these conditions the proposed CTGWO algorithm yields the best trade-off between protection and invisibility.

NSGA-II yields the best result in terms of invisibility, that is, $5.363 \times 10^{-4}$, at the expense of a mean value $f_1/f_2 = 8.063 \times 10^{-4}$ which is larger than the score obtained by CTGWO. However the interest of NSGA-II is to enable an end-user to choose to favor one criterion (for instance invisibility) rather than the other. When $n_{max} = 1000$ and $Q=200$ affording therefore 8 times more trials of each objective function, the best protection reaches $1.233 \times 10^{-4}$, with an associated mean value $f_1/f_2 = 4.021 \times 10^{-4}$. In Figures 8 and 10, we can check the coherence of the shape of the wavefront and the score values obtained. A less realistic manner to improve the results is to enable the search for any real permittivity value between 7 and 12.
Cloak design

<figure><img src="image_9.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

**Fig. 8.** Column 1: in (a) best invisibility cloak, in (b) best compromise between invisibility and protection cloak, and in (c), best protection cloak; column 2: in (d) best invisibility result, in (e) best compromise between invisibility and protection result, and in (f), best protection result. See Table 3 for numerical values.

In the continuous case, still with $2.5 \times 10^4$ trials of both $f_1$ and $f_2$, we reach a best protection criterion $7.845 \times 10^{-5}$, a best invisibility criterion $3.274 \times 10^{-5}$, and a best mean value $1.416 \times 10^{-4}$. The results are improved but the corresponding physical constraints are very much strengthened, as we assume we can choose any material for any voxel in the cloak. In the continuous case, with $n_{max}=1000$ and $Q = 200$, we reach very low best protection ($6.30584 \times 10^{-5}$), best invisibility ($2.59173 \times 10^{-5}$), and best mean

($1.10285 \times 10^{-4}$) values. The visual results in the continuous case (see Figs. 12 and 14) are very convincing, and clearly illustrate the difference between the best protection, the best invisibility, the best compromise, and the best protection cases. Finally, the tolerance with respect to the operating freespace incident wavelength showed a broadband behavior in terms of protection and a rather narrowband behavior for the invisibility criterion, as shown by the comparison to random cloaks in Figure 15.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Column 1: in (a) best invisibility cloak, in (b) best compromise between invisibility and protection cloak, and in (c), best protection cloak; column 2: in (d) best invisibility result, in (e) best compromise between invisibility and protection result, and in (f), best protection result. See Table 3 for numerical values.</figcaption></figure>

In the continuous case, still with $2.5 \times 10^4$ trials of both $f_1$ and $f_2$, we reach a best protection criterion $7.845 \times 10^{-5}$, a best invisibility criterion $3.274 \times 10^{-5}$, and a best mean value $1.416 \times 10^{-4}$. The results are improved but the corresponding physical constraints are very much strengthened, as we assume we can choose any material for any voxel in the cloak. In the continuous case, with $n_{max}=1000$ and $Q = 200$, we reach very low best protection ($6.30584 \times 10^{-5}$), best invisibility ($2.59173 \times 10^{-5}$), and best mean

($1.10285 \times 10^{-4}$) values. The visual results in the continuous case (see Figs. 12 and 14) are very convincing, and clearly illustrate the difference between the best protection, the best invisibility, the best compromise, and the best protection cases. Finally, the tolerance with respect to the operating freespace incident wavelength showed a broadband behavior in terms of protection and a rather narrowband behavior for the invisibility criterion, as shown by the comparison to random cloaks in Figure 15.
<figure><img src="image_27.png" /><figcaption>Fig. 9. Pareto front for NSGA-II Table with 250 iterations, 100 research agents, with four selected solutions out of a total of 175.</figcaption></figure>

<table><caption>**Table 4.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}$ =1000, $Q$ = 200. See Figure 7 for the Pareto front and Figure 8 for 2D plots of cloak design and wave fields.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>1.234 × 10<sup>-4</sup></td><td>6.808 × 10<sup>-4</sup></td><td>4.021 × 10<sup>-4</sup></td><td>Figure 10f</td></tr><tr><td>Best invisibility</td><td>1.006 × 10<sup>-3</sup></td><td>4.312 × 10<sup>-4</sup></td><td>7.188 × 10<sup>-4</sup></td><td>Figure 10d</td></tr><tr><td>Best compromise</td><td>4.060 × 10<sup>-4</sup></td><td>5.198 × 10<sup>-4</sup></td><td>4.629 × 10<sup>-4</sup></td><td>Figure 10e</td></tr></tbody></table>

<table><caption>**Table 5.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}$ = 250, $Q$ = 100.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>7.845 × 10<sup>-5</sup></td><td>4.302 × 10<sup>-4</sup></td><td>2.543 × 10<sup>-4</sup></td><td>Figure 12f</td></tr><tr><td>Best invisibility</td><td>4.419 × 10<sup>-4</sup></td><td>3.274 × 10<sup>-5</sup></td><td>2.373 × 10<sup>-4</sup></td><td>Figure 12d</td></tr><tr><td>Best compromise</td><td>1.666 × 10<sup>-4</sup></td><td>1.167 × 10<sup>-4</sup></td><td>1.416 × 10<sup>-4</sup></td><td>Figure 12e</td></tr></tbody></table>

<table><caption>**Table 6.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}$=1000, $Q$ = 200.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>6.306 × 10<sup>-5</sup></td><td>2.341 × 10<sup>-4</sup></td><td>1.486 × 10<sup>-4</sup></td><td>Figure 14f</td></tr><tr><td>Best invisibility</td><td>7.782 × 10<sup>-4</sup></td><td>2.592 × 10<sup>-5</sup></td><td>4.021 × 10<sup>-4</sup></td><td>Figure 14d</td></tr><tr><td>Best compromise</td><td>1.530 × 10<sup>-4</sup></td><td>6.741 × 10<sup>-5</sup></td><td>1.102 × 10<sup>-4</sup></td><td>Figure 14e</td></tr></tbody></table>

## 7 Conclusion

In this work, we address the issue of an electromagnetic cloak's design in the transverse magnetic TM polarization (whereby the magnetic field is perpendicular to the $(xy)$-plane containing the computational domain). This polarization has been chosen as the model can then be

adapted to plasmonics by assuming some Drude-like dispersion in the permittivity [36]. More precisely, this means that $\sigma$ in equation (1) should be frequency dependent, and $\chi$ is set to 1. However, the transverse electric polarization would be also worth investigating, in that case $\sigma$ is set to 1 in equation (1), and $\chi$ plays the role of the permittivity.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Pareto front for NSGA-II Table with 250 iterations, 100 research agents, with four selected solutions out of a total of 175.</figcaption></figure>

<table><caption>**Table 4.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}$ =1000, $Q$ = 200. See Figure 7 for the Pareto front and Figure 8 for 2D plots of cloak design and wave fields.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>1.234 × 10<sup>-4</sup></td><td>6.808 × 10<sup>-4</sup></td><td>4.021 × 10<sup>-4</sup></td><td>Figure 10f</td></tr><tr><td>Best invisibility</td><td>1.006 × 10<sup>-3</sup></td><td>4.312 × 10<sup>-4</sup></td><td>7.188 × 10<sup>-4</sup></td><td>Figure 10d</td></tr><tr><td>Best compromise</td><td>4.060 × 10<sup>-4</sup></td><td>5.198 × 10<sup>-4</sup></td><td>4.629 × 10<sup>-4</sup></td><td>Figure 10e</td></tr></tbody></table>

<table><caption>**Table 5.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}$ = 250, $Q$ = 100.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>7.845 × 10<sup>-5</sup></td><td>4.302 × 10<sup>-4</sup></td><td>2.543 × 10<sup>-4</sup></td><td>Figure 12f</td></tr><tr><td>Best invisibility</td><td>4.419 × 10<sup>-4</sup></td><td>3.274 × 10<sup>-5</sup></td><td>2.373 × 10<sup>-4</sup></td><td>Figure 12d</td></tr><tr><td>Best compromise</td><td>1.666 × 10<sup>-4</sup></td><td>1.167 × 10<sup>-4</sup></td><td>1.416 × 10<sup>-4</sup></td><td>Figure 12e</td></tr></tbody></table>

<table><caption>**Table 6.** Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak. $n_{max}$ =1000, $Q$ = 200.</caption><thead><tr><th></th><th>Protection f<sub>1</sub></th><th>Invisibility f<sub>2</sub></th><th>f<sub>1</sub>+f<sub>2</sub>/2</th><th>Links to figures</th></tr></thead><tbody><tr><td>Best protection</td><td>6.306 × 10<sup>-5</sup></td><td>2.341 × 10<sup>-4</sup></td><td>1.486 × 10<sup>-4</sup></td><td>Figure 14f</td></tr><tr><td>Best invisibility</td><td>7.782 × 10<sup>-4</sup></td><td>2.592 × 10<sup>-5</sup></td><td>4.021 × 10<sup>-4</sup></td><td>Figure 14d</td></tr><tr><td>Best compromise</td><td>1.530 × 10<sup>-4</sup></td><td>6.741 × 10<sup>-5</sup></td><td>1.102 × 10<sup>-4</sup></td><td>Figure 14e</td></tr></tbody></table>

## 7 Conclusion

In this work, we address the issue of an electromagnetic cloak's design in the transverse magnetic TM polarization (whereby the magnetic field is perpendicular to the $(xy)$-plane containing the computational domain). This polarization has been chosen as the model can then be

adapted to plasmonics by assuming some Drude-like dispersion in the permittivity [36]. More precisely, this means that $\sigma$ in equation (1) should be frequency dependent, and $\chi$ is set to 1. However, the transverse electric polarization would be also worth investigating, in that case $\sigma$ is set to 1 in equation (1), and $\chi$ plays the role of the permittivity.
<figure><img src="image_28.png" /><figcaption>**Fig. 10.** Column 1: in (a) best invisibility cloak, in (b) best compromise between invisibility and protection cloak, and in (c), best protection cloak; column 2: in (d) best invisibility result, in (e) best compromise between invisibility and protection result, and in (f), best protection result.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 11.** In red, the four selected solutions.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Column 1: in (a) best invisibility cloak, in (b) best compromise between invisibility and protection cloak, and in (c), best protection cloak; column 2: in (d) best invisibility result, in (e) best compromise between invisibility and protection result, and in (f), best protection result.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** In red, the four selected solutions.</figcaption></figure>
Cloak design

<figure><img src="image_30.png" /><figcaption>(a)</figcaption></figure>

Permittivity
$ε_r = 7.39$

y

<figure><img src="image_31.png" /><figcaption>(d)</figcaption></figure>

utot - real part
$0$

y

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /><figcaption>(b)</figcaption></figure>

Permittivity
$ε_r = 12.7$

y

<figure><img src="image_34.png" /><figcaption>(e)</figcaption></figure>

utot - real part
$0$

y

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /><figcaption>(c)</figcaption></figure>

Permittivity
$ε_r = 12.8$

y

<figure><img src="image_37.png" /><figcaption>(f)</figcaption></figure>

utot - real part
$0$

y

<figure><img src="image_38.png" /></figure>

**Fig. 12.** Column 1: in (a) best invisibility cloak, in (b) best compromise between invisibility and protection cloak, and in (c), best protection cloak; column 2: in (d) best invisibility result, in (e) best compromise between invisibility and protection result, and in (f), best protection result.

<figure><img src="image_39.png" /></figure>

**Fig. 13.** In red, the four selected solutions.
Cloak design

<figure><img src="image_13.png" /><figcaption>(a)</figcaption></figure>

Permittivity
$\epsilon_r = 7.39$

$\gamma$

utot - real part
$0.1$

$\gamma$

Wave propagation field

<figure><img src="image_14.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b)</figcaption></figure>

Permittivity
$\epsilon_r = 7.37$

$\gamma$

utot - real part
$0.1$

$\gamma$

<caption>(c)</caption>

Permittivity
$\epsilon_r = 7.41$

$\gamma$

utot - real part
$0.1$

$\gamma$

<caption>(f)</caption>

**Fig. 12.** Column 1: in (a) best invisibility cloak, in (b) best compromise between invisibility and protection cloak, and in (c), best protection cloak; column 2: in (d) best invisibility result, in (e) best compromise between invisibility and protection result, and in (f), best protection result.

<figure><img src="image_17.png" /></figure>

**Fig. 13.** In red, the four selected solutions.
Cloak design

<figure><img src="image_40.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_41.png" /></figure>

Wave propagation field

<figure><img src="image_42.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_49.png" /></figure>

**Fig. 14.** Column 1: in (a) Best invisibility cloak, in (b) Best compromise between invisibility and protection cloak, and in (c), Best protection cloak; Column 2: in (d) Best invisibility result, in (e) Best compromise between invisibility and protection result, and in (f), Best protection result.
Cloak design

<figure><img src="image_18.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_19.png" /></figure>

Wave propagation field

<figure><img src="image_20.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_27.png" /></figure>

**Fig. 14.** Column 1: in (a) Best invisibility cloak, in (b) Best compromise between invisibility and protection cloak, and in (c), Best protection cloak; Column 2: in (d) Best invisibility result, in (e) Best compromise between invisibility and protection result, and in (f), Best protection result.
<figure><img src="image_50.png" /><figcaption>**Fig. 15.** Spectral tolerance of the best candidates in terms of protection ($C_1$, left panel) and invisibility ($C_2$, right panel). The wavelength in abscissa is normalized by $λ_0$, the freespace wavelength targeted for the optimization process. The blue curves correspond to the CTGWO method, the orange ones to T-NSGA-II, the green ones to T-NSGA-II+, the red ones to C-NSGA-II and the purple ones to C-NSGA-II+. The five grey (resp. black) curves in each plot represent the spectral behavior of $C_1$ and $C_2$ for random cloaks obtained by filling each voxel of the design space by an arbitrary integer (resp. real) value in {7, 10, 12} (resp. [7:12]).</figcaption></figure>

Our objective is here to achieve the best compromise between protection and invisibility for TM waves. In other words, we are looking for the best trade-off between protection while considering the wave amplitude inside the cloak, and invisibility while considering the wave behavior outside of the cloak. This is a large scale bi-objective optimization problem. We propose two approaches: in the first one we transform this problem into a mono-objective optimization problem and seek for the best mean value of protection and invisibility criteria. In the second one we look for the best protection, the best invisibility, and the best trade-off with a bi-objective optimization algorithm. GWO is a well known mono-objective optimization algorithm which reaches the desired solution with a reduced number of iterations. We propose a novel mono-objective version of GWO, namely the chaotic ternary GWO, with three main innovations: ad hoc update rules to face ternary search spaces, evolving map and chaotic sequences to improve exploration abilities, and division of the pack into two groups to improve diversity. We apply this algorithm, and the comparative GWO and amixedGWO (both mono-objective) as well as the NSGA-II (bi-objective) to solve the considered cloak design problem. In the considered cloaking application, and with the help of $2.5 \times 10^4$ evaluations of these criteria, the proposed CTGWO algorithm yields the best mean value, that is, $5.128 \times 10^{-4}$, hence the best ‘trade-off’, between protection and invisibility.

It surpasses GWO ($29.871 \times 10^{-4}$), amixedGWO ($6.526 \times 10^{-4}$); and the best trade-off provided by NSGA-II ($5.306 \times 10^{-4}$). A possible prospect could consist in altering the ternary map in CTGWO to favor one particular

material, among the three which are chosen to build the cloak. This could help in implementing the cloak with a preferred material.

## Funding

This research was funded by BTU Trust and ANRT, grant number [CONVENTION CIFRE N° 2020/0078].

## Conflicts of interest

No author associated with this paper has disclosed any potential or pertinent conflict which may be perceived to have impending conflict with this work.

## Data availability statement

This article has no associated available data.

## Author contribution statement

**Ronald Aznavourian:** Conceptualization, Software, Validation, Writing. **Guillaume Demesy:** Software, Validation, Formal analysis.

**Sébastien Guenneau:** Conceptualization, Validation, Supervision, Final preparing. **Julien Marot:** Conceptualization, Software, Editing, Writing, Supervision, Final preparing.

## References

1. U. Leonhardt, T. Philbin, *Geometry and light: the science of invisibility* (Courier Corporation, 2010)
<figure><img src="image_28.png" /><figcaption>**Fig. 15.** Spectral tolerance of the best candidates in terms of protection ($C_1$, left panel) and invisibility ($C_2$, right panel). The wavelength in abscissa is normalized by $λ_0$, the freespace wavelength targeted for the optimization process. The blue curves correspond to the CTGWO method, the orange ones to T-NSGA-II, the green ones to T-NSGA-II+, the red ones to C-NSGA-II and the purple ones to C-NSGA-II+. The five grey (resp. black) curves in each plot represent the spectral behavior of $C_1$ and $C_2$ for random cloaks obtained by filling each voxel of the design space by an arbitrary integer (resp. real) value in {7, 10, 12} (resp. [7:12]).</figcaption></figure>

Our objective is here to achieve the best compromise between protection and invisibility for TM waves. In other words, we are looking for the best trade-off between protection while considering the wave amplitude inside the cloak, and invisibility while considering the wave behavior outside of the cloak. This is a large scale bi-objective optimization problem. We propose two approaches: in the first one we transform this problem into a mono-objective optimization problem and seek for the best mean value of protection and invisibility criteria. In the second one we look for the best protection, the best invisibility, and the best trade-off with a bi-objective optimization algorithm. GWO is a well known mono-objective optimization algorithm which reaches the desired solution with a reduced number of iterations. We propose a novel mono-objective version of GWO, namely the chaotic ternary GWO, with three main innovations: ad hoc update rules to face ternary search spaces, evolving map and chaotic sequences to improve exploration abilities, and division of the pack into two groups to improve diversity. We apply this algorithm, and the comparative GWO and amixedGWO (both mono-objective) as well as the NSGA-II (bi-objective) to solve the considered cloak design problem. In the considered cloaking application, and with the help of $2.5 \times 10^4$ evaluations of these criteria, the proposed CTGWO algorithm yields the best mean value, that is, $5.128 \times 10^{-4}$, hence the best ‘trade-off’, between protection and invisibility.

It surpasses GWO ($29.871 \times 10^{-4}$), amixedGWO ($6.526 \times 10^{-4}$); and the best trade-off provided by NSGA-II ($5.306 \times 10^{-4}$). A possible prospect could consist in altering the ternary map in CTGWO to favor one particular

material, among the three which are chosen to build the cloak. This could help in implementing the cloak with a preferred material.

## Funding

This research was funded by BTU Trust and ANRT, grant number [CONVENTION CIFRE N° 2020/0078].

## Conflicts of interest

No author associated with this paper has disclosed any potential or pertinent conflict which may be perceived to have impending conflict with this work.

## Data availability statement

This article has no associated available data.

## Author contribution statement

**Ronald Aznavourian:** Conceptualization, Software, Validation, Writing. **Guillaume Demesy:** Software, Validation, Formal analysis.

**Sébastien Guenneau:** Conceptualization, Validation, Supervision, Final preparing. **Julien Marot:** Conceptualization, Software, Editing, Writing, Supervision, Final preparing.

## References

1. U. Leonhardt, T. Philbin, *Geometry and light: the science of invisibility* (Courier Corporation, 2010)
2. J.B. Pendry, D. Schurig, D.R. Smith, Controlling electromagnetic fields, Science **312**, 1780 (2006)

3. M. Cassier, G.W. Milton, Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking, J. Math. Phy, **58**, 071504 (2017)

4. J. Li, J.B. Pendry, Hiding under the carpet: a new strategy for cloaking, Phys. Rev. Lett. **101**, 203901 (2008)

5. R. Kohn, M. Vogelius, Determining conductivity by boundary measurements, Commun. Pure Appl. Math. **37**, 289 (1984)

6. J.M. Lee, G. Uhlmann, Determining anisotropic real-analytic conductivities by boundary measurements, Commun. Pure Appl. Math. **42**, 1097 (1989)

7. R.V. Kohn, H. Shen, M.S. Vogelius, M.I. Weinstein, Cloaking via change of variables in electric impedance tomography, Inverse Probl. **24**, 015016 (2008)

8. J. Sylvester, G. Uhlmann, A global uniqueness theorem for an inverse boundary value problem, Ann. Math. **125**, 153 (1987)

9. A. Greenleaf, M. Lassas, G. Uhlmann, Anisotropic conductivities that cannot be detected by EIT, Physiol. Meas. **24**, 413 (2003)

10. A. Alù, N. Engheta, Achieving transparency with plasmonic and metamaterial coatings. Phys. Rev. E **72**, 016623 (2005)

11. P.-Y. Chen, J. Soric, A. Alù, Invisibility and cloaking based on scattering cancellation, Adv. Mater. **24**, OP281 (2012)

12. G.W. Milton, N.-A.P. Nicorovici, R.C. McPhedran, V.A. Podolski, A proof of superlensing in the quasistatic regime, and limitations of superlenses in this regime due to anomalous localized resonance, Proc. Roy. Soc. A: Math. Phys. Eng. Sci. **461**, 3999 (2005)

13. G.W. Milton, M. Briane, J.R. Willis, On cloaking for elasticity and physical equations with a transformation invariant form, N. J. Phys. **8**, 248 (2006)

14. M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. **69**, 635 (1999)

15. B. Vial, Y. Hao, Yang, Topology optimized all-dielectric cloak: design, performances and modal picture of the invisibility effect, Opt. Express **23**, 23551 (2015)

16. J. Andkjær, O. Sigmund, Topology optimized low-contrast all-dielectric optical cloak, Appl. Phys. Lett. **98**, 021112 (2011)

17. L. Pomot, C. Payan, M. Remillieux, S. Guenneau, Acoustic cloaking: geometric transform, homogenization and a genetic algorithm, Wave Motion **92**, 102413 (2020)

18. F.L. Teixeira, W.C. Chew, General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media, IEEE Microwave Guided Wave Lett. **8**, 223 (1998)

19. F. Zolla, G. Renversez, A. Nicolet, B. Kuhlme, S. Guenneau, D. Felbacq, A. Argyros, S. Leon-Saval, *Foundations of photonic crystal fibres*, 2nd ed. (Imperial College Press, 2012)

20. E. Kuci, F. Henrotte, P. Duysinx, C. Geuzaine, Design sensitivity analysis for shape optimization based on the Lie derivative, Comput. Methods Appl. Mech. Eng. **317**, 702 (2017)

21. C. Geuzaine, ONELAB. http://onelab.info/, (accessed November 14, 2022)

22. GetDP documentation. http://getdp.info/#Documentation (accessed November 14, 2022)

23. J. Smajic, C. Hafner, D. Erni, Optimization of photonic crystal structures, J. Opt. Soc. Am. A **21**, 2223 (2004)

24. T. Ross, G. Cormier, Particle swarm optimization for ellipsometric data inversion of samples having an arbitrary number of layers, J. Opt. Soc. Am. A **27**, 319 (2010)

25. P. Zhang, Y. Qian, Q. Qian, Multi-objective optimization for materials design with improved NSGA-II, Mater. Today Commun. **28**, 102709 (2021)

26. W. Costa, H. Camporez, M. Pontes, M. Segatto, H. Rocha, J. Silva, M. Hinrichs, A. Paraskevopoulos, V. Jungnickel, R. Freund, Increasing the power and spectral efficiencies of an OFDM-based VLC system through multi-objective optimization, J. Opt. Soc. Am. A **40**, 1268 (2023)

27. N. Srinivas, K. Deb, Multiobjective optimization using nondomi- nated sorting in genetic algorithms, Evol. Comput. **2**, 221 (1994)

28. K. Deb, A. Pratap, S. Agarwal, T.A.M.T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Trans. Evol. Comput. **6**, 182 (2002)

29. S. Mirjalili, S. Mohammad Mirjalili, A, Lewis, Grey wolf optimizer, Adv. Eng. Softw. **69**, 46 (2014)

30. N. Mittal, U. Singh, B. Singh Sohi, Modified grey wolf optimizer for global engineering optimization, Appl. Comput. Intell. Soft Comput. **2016**, 7950348 (2016)

31. B. Martin, J. Marot, S. Bourennane, Mixed grey wolf optimizer for the joint denoising and unmixing of multispectral images, Appl. Soft Comput. **74**, 385 (2019)

32. F. Waschkowski, Y. Zhao, R. Sandberg, J. Klewicki, Multi-objective CFD-driven development of coupled turbulence closure models, J. Comput. Phys. **452**, 110922 (2022)

33. Y. Chen, S. Zhou, Q. Li, Multiobjective topology optimization for finite periodic structures, Comput. Struct. **88**, 806 (2010)

34. E. Emary, H.M. Zawbaa, A. Ella Hassanien, Binary grey wolf optimization approaches for feature selection, Neurocomputing **172**, 371 (2016)

35. P. Hu, J.-S. Pan, S.-C. Chu, Improved Binary Grey Wolf Optimizer and Its application for feature selection, Knowledge Based Syst. **195**, 105746 (2020)

36. M. Kadic, S. Guenneau, S. Enoch, P.A. Huidobro, L. Martin- Moreno, F.J. García-Vidal, J. Renger, R. Quidant, Transformation plasmonics, Nanophotonics **1**, 51 (2012)

Cite this article as: Ronald Aznavourian, Guillaume Demesy, Sébastien Guenneau, Julien Marot, Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility, EPJ Appl. Metamat. **11**, 11 (2024)
2. J.B. Pendry, D. Schurig, D.R. Smith, Controlling electromagnetic fields, Science **312**, 1780 (2006)

3. M. Cassier, G.W. Milton, Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking, J. Math. Phy, **58**, 071504 (2017)

4. J. Li, J.B. Pendry, Hiding under the carpet: a new strategy for cloaking, Phys. Rev. Lett. **101**, 203901 (2008)

5. R. Kohn, M. Vogelius, Determining conductivity by boundary measurements, Commun. Pure Appl. Math. **37**, 289 (1984)

6. J.M. Lee, G. Uhlmann, Determining anisotropic real-analytic conductivities by boundary measurements, Commun. Pure Appl. Math. **42**, 1097 (1989)

7. R.V. Kohn, H. Shen, M.S. Vogelius, M.I. Weinstein, Cloaking via change of variables in electric impedance tomography, Inverse Probl. **24**, 015016 (2008)

8. J. Sylvester, G. Uhlmann, A global uniqueness theorem for an inverse boundary value problem, Ann. Math. **125**, 153 (1987)

9. A. Greenleaf, M. Lassas, G. Uhlmann, Anisotropic conductivities that cannot be detected by EIT, Physiol. Meas. **24**, 413 (2003)

10. A. Alù, N. Engheta, Achieving transparency with plasmonic and metamaterial coatings. Phys. Rev. E **72**, 016623 (2005)

11. P.-Y. Chen, J. Soric, A. Alù, Invisibility and cloaking based on scattering cancellation, Adv. Mater. **24**, OP281 (2012)

12. G.W. Milton, N.-A.P. Nicorovici, R.C. McPhedran, V.A. Podolski, A proof of superlensing in the quasistatic regime, and limitations of superlenses in this regime due to anomalous localized resonance, Proc. Roy. Soc. A: Math. Phys. Eng. Sci. **461**, 3999 (2005)

13. G.W. Milton, M. Briane, J.R. Willis, On cloaking for elasticity and physical equations with a transformation invariant form, N. J. Phys. **8**, 248 (2006)

14. M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. **69**, 635 (1999)

15. B. Vial, Y. Hao, Yang, Topology optimized all-dielectric cloak: design, performances and modal picture of the invisibility effect, Opt. Express **23**, 23551 (2015)

16. J. Andkjær, O. Sigmund, Topology optimized low-contrast all-dielectric optical cloak, Appl. Phys. Lett. **98**, 021112 (2011)

17. L. Pomot, C. Payan, M. Remillieux, S. Guenneau, Acoustic cloaking: geometric transform, homogenization and a genetic algorithm, Wave Motion **92**, 102413 (2020)

18. F.L. Teixeira, W.C. Chew, General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media, IEEE Microwave Guided Wave Lett. **8**, 223 (1998)

19. F. Zolla, G. Renversez, A. Nicolet, B. Kuhlme, S. Guenneau, D. Felbacq, A. Argyros, S. Leon-Saval, *Foundations of photonic crystal fibres*, 2nd ed. (Imperial College Press, 2012)

20. E. Kuci, F. Henrotte, P. Duysinx, C. Geuzaine, Design sensitivity analysis for shape optimization based on the Lie derivative, Comput. Methods Appl. Mech. Eng. **317**, 702 (2017)

21. C. Geuzaine, ONELAB. http://onelab.info/, (accessed November 14, 2022)

22. GetDP documentation. http://getdp.info/#Documentation (accessed November 14, 2022)

23. J. Smajic, C. Hafner, D. Erni, Optimization of photonic crystal structures, J. Opt. Soc. Am. A **21**, 2223 (2004)

24. T. Ross, G. Cormier, Particle swarm optimization for ellipsometric data inversion of samples having an arbitrary number of layers, J. Opt. Soc. Am. A **27**, 319 (2010)

25. P. Zhang, Y. Qian, Q. Qian, Multi-objective optimization for materials design with improved NSGA-II, Mater. Today Commun. **28**, 102709 (2021)

26. W. Costa, H. Camporez, M. Pontes, M. Segatto, H. Rocha, J. Silva, M. Hinrichs, A. Paraskevopoulos, V. Jungnickel, R. Freund, Increasing the power and spectral efficiencies of an OFDM-based VLC system through multi-objective optimization, J. Opt. Soc. Am. A **40**, 1268 (2023)

27. N. Srinivas, K. Deb, Multiobjective optimization using nondomi- nated sorting in genetic algorithms, Evol. Comput. **2**, 221 (1994)

28. K. Deb, A. Pratap, S. Agarwal, T.A.M.T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Trans. Evol. Comput. **6**, 182 (2002)

29. S. Mirjalili, S. Mohammad Mirjalili, A, Lewis, Grey wolf optimizer, Adv. Eng. Softw. **69**, 46 (2014)

30. N. Mittal, U. Singh, B. Singh Sohi, Modified grey wolf optimizer for global engineering optimization, Appl. Comput. Intell. Soft Comput. **2016**, 7950348 (2016)

31. B. Martin, J. Marot, S. Bourennane, Mixed grey wolf optimizer for the joint denoising and unmixing of multispectral images, Appl. Soft Comput. **74**, 385 (2019)

32. F. Waschkowski, Y. Zhao, R. Sandberg, J. Klewicki, Multi-objective CFD-driven development of coupled turbulence closure models, J. Comput. Phys. **452**, 110922 (2022)

33. Y. Chen, S. Zhou, Q. Li, Multiobjective topology optimization for finite periodic structures, Comput. Struct. **88**, 806 (2010)

34. E. Emary, H.M. Zawbaa, A. Ella Hassanien, Binary grey wolf optimization approaches for feature selection, Neurocomputing **172**, 371 (2016)

35. P. Hu, J.-S. Pan, S.-C. Chu, Improved Binary Grey Wolf Optimizer and Its application for feature selection, Knowledge Based Syst. **195**, 105746 (2020)

36. M. Kadic, S. Guenneau, S. Enoch, P.A. Huidobro, L. Martin- Moreno, F.J. García-Vidal, J. Renger, R. Quidant, Transformation plasmonics, Nanophotonics **1**, 51 (2012)

Cite this article as: Ronald Aznavourian, Guillaume Demesy, Sébastien Guenneau, Julien Marot, Electromagnetic cloak design with mono-objective and bi-objective optimizers: seeking the best tradeoff between protection and invisibility, EPJ Appl. Metamat. **11**, 11 (2024)