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

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.


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 nonunique conductivities sharing the same boundary measurements [9].

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 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.

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.

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 ; . . .

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: 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 l 0 = 2p/k 0 ), and s and x represent the scalar material properties.In the context of acoustic pressure waves in isotropic non-viscous fluids, s is the inverse of density and x the inverse of compressibility, and p is the amplitude of the pressure wave.For anti-plane shear waves in isotropic solids, s is the shear modulus and x 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 s is the inverse of the relative permittivity and x the relative permeability.The relative permittivity will be denoted by e 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 l = 500 nm the side of a triangle mesh has length l/6 in the freespace; and l/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 s can be written as a 2 Â 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 s 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 g 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 g 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: 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 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 impedanced 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 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 (x):ℝ P ↦ ℝ + and f 2 (x):ℝ P ↦ ℝ + .Vector x contains the parameter values K 1 , K 2 , … , K i , … , K P .

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 ≥ 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(.) 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(.) to be minimized.In the case of bi-objective optimization, there are two functions f 1 (.) and f 2 (.) to be minimized.Both GWO and NSGA-II are agent-based algorithms i.
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 x q (n)=[K 1 ,K 2 ,...,K P ] T .In Section 3.1, we give a background on a singleobjective 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].

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 a, b, and d All other agents are the v wolves.
x a (n), x b (n), and x d (n) denote the position of the leaders a, b and d respectively, at iteration n.
The position of any wolf at iteration n+1 is calculated as: It results from the equal contribution of the a, b and d wolves.These contributions are computed at each iteration iter as follows, for any leader l either a, b or d: with: d l ðnÞ ¼ j c x l ðnÞ À x q ðnÞj; j•j denoting absolute value.
The vectors b and c are calculated as b = 2a • 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 b i of b is provided by: the same whatever i.
The component d i l ðnÞ of d l (n) is provided by: where r 1 and r 2 are two random values between 0 and 1; x i q ðnÞ is the i th component of the q th agent at iteration n; x i l is the i th component of leader l.
The component y i l ðnÞ of y l (n) is provided by: 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 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: 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 a ð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 a (n max ), but also the convergence curve.

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: and 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 0 < 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 0 for 'max' For each solution of index q with 1< q < Q 0 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.

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.

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.

Contribution of a leader
We remind that in the continuous case, the contribution y i l ðnÞ of a leader l is computed as in equation (8).In this case y i l ðnÞ ¼ x i l À b i b i l ðnÞ depends on the product b i d i l ðnÞ which decreases to 0 simultaneously with a, reaching 0 when n ¼ n max : So y i l ðnÞ is a real value which gets closer to x i l across the iterations.In the following we propose to compute a contribution y i l ðnÞ which is a real value between 0 and 2. We propose, as expression of the contribution of leader l: where b i is defined as in equation ( 9) and d i l ð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 i l ðnÞ are between À8 and 8.So the values of x i l À b i d i l ðnÞ can be out of the interval [0:2].The 'modulo' operator, denoted by mod is meant to enforce the contributions y i l ðnÞ to remain into the interval [0:2].It is defined as follows: whatever the values z 1 ∈ ℝ and z 2 ∈ ℝ Ã : where |z 1 | denotes the largest value in ℤ which is smaller than the scalar z 1 ∈ ℝ.

Weighted contribution of the leaders
We denote by y i ðnÞ the weighted contribution of four leaders a, b, d, and r: In equation ( 12), leader r is a wolf which is selected at random among the wolf pack.
We notice that y i ðnÞ gets closer to when the iteration index increases.
We can now derive a ternary rule which updates the position of any wolf in a ternary search space.

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: 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 ' u and ' d : Function ' u separates the uppermost part of the map from the rest of the map; and function ' 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 ' u (resp.below ' 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: 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 ' p ð0; c 1 ; c 2 ; c 3 ; c 4 ; aÞ ≥ 0 and ' p ð2; c 1 ; c 2 ; c 3 ; c 4 ; aÞ 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: 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: and in both equations ( 15) and ( 16).So: and The functions ' u (y, a) and ' d (y, a) defined in equations ( 15) and ( 16) can then be used in equation ( 13) to get the 'Power' transform map.

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 ' u and ' 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 y i ð n Þ ¼ 0 ðresp: y i ð n Þ ¼ 2Þ; -Whatever the value of a, either the values 0, 1, or 2 may be selected when y i ðnÞ= 1.
Indeed, as defined in equation ( 10) the values of y i 1 ð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.

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: 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 h which depends on the score of the agent: for the worst half of the agents (associated with the largest scores), h q ¼ 2 3 ; for the best half of the agents (associated with the smallest scores), h q ¼ 3 2 : Secondly, inserting a chaotic sequence G ðc 1 q ; nÞ enhances the variability of the behavior of the agents.c 1 q 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 1 q or each agent permits to emphasize diversity in the displacement of the agents.Meanwhile, we choose a sequence with one attractor to ensure that G ðc 1 q ; n max Þ is the same whatever q.

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 G ðc 1 q ; nÞ are positive, in the interval [0:1]; -The last value should be G ðc 1 q ; n max Þ ¼ 0 whatever the initial value G c 1 q ; 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 G 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: The number of attractors for this sequence depends on the value of k. Figure 3 shows chaotic sequences with one (Fig. 3a), two (Fig. 3b), or an infinite number of attractors (Fig. 3c).We chose k = 2.8 in Figure 3a, k = 3.2 in Figure 3b, k = 3.99 in Figure 3c.In Figures 3a-3c, each 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 k = 2.8.
For any agent q and for n = 1,...,n max , we define the ad hoc chaotic sequence G ðc 1 q ; nÞ as follows, based on the logistic sequence of equation ( 20): Hence, for any agent qG ðc 1 q ; n max Þ ¼ 0. We notice that the initial term c 1 q is defined as follows: 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 1 q 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 h q ¼ 3 2 and rather elevated values of a, and the second family with h q ¼ 2 3 and smaller values of a, which tend more rapidly towards 0.
5 Results: mono-objective and bi-objective approaches

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 In Section 5.2 we display the results obtained with a monoobjective approach.The criterion which is minimized is f (x):ℝ P ↦ ℝ + : fðxÞ ¼ 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 e is {7,10,12}.CTGWO In Section 5.3, we display the results obtained with a biobjective approach.The couple of criteria which are jointly minimized are f 1 (x) :ℝ P ↦ ℝ + , and f 2 (x):ℝ P ↦ ℝ + .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.

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 x a ðn max Þ .For instance, the solution provided by CTGWO is x a (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.

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 ≃ 0.003 and p c = 0.9; mutation and crossover distribution indices are h m = 20 and h c = 15.

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.

Bi-objective ternary optimization with 1000 iterations, 200 research agents
In this subsubsection, 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.

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.

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.

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 l 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 l 0 : for both C 1 (protection) and C 2 (invisibility) these criteria are 20 times smaller for the optimized cloaks: 1.5 Â 10 −4 vs 2.2 Â 10 −3 for C 1 ; and 6 Â 10 −5 vs 1.2 Â 10 −3 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 l/l 0 ∈ [0.95 : 1.03]), as opposed to the protection which turns out to be more broadband (the cloak optimized with C-NSGA-II+ lead to a better protection than any random realizations for l/l 0 ∈ [0.82 : 1.22]).We note that the optimization method operating over a continuous space lead to a more broadband response (see red and purple curves in Fig. 15).

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 Â 10 −4 and score, that is, f 1 þf 2 2 ¼ 5:128 Â 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 Â 10 −4 , at the expense of a mean value ¼ 8:063 Â 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Â10 −4 , with an associated mean value ¼ 4:021 Â 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 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    In the continuous case, still with 2.5 ×10 4 trials of both f 1 and f 2 , we reach a best protection criterion 7.845 Â 10 À5 , a best invisibility criterion 3.274 Â 10 À5 , and a best mean value 1.416 Â 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 Â 10 À5 ), best invisibility (2.59173 Â 10 −5 ), and best mean (1.10285 Â 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.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 s in equation ( 1) should be frequency dependent, and x is set to 1.However, the transverse electric polarization would be also worth investigating, in that case s is set to 1 in equation ( 1), and x plays the role of the permittivity.11, 11 (2024) ; 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 b is a vector whose K components are equal to a 1 b 1 , a 2 , b 2 , … , a K b 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 |.

Fig. 1 .
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).
for the Pareto front and Figure8for 2D plots of cloak design and wave fields.

Fig. 7 .
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 Table3.

Fig. 8 .
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 Table3for numerical values.

Fig. 9 .
Fig. 9. Pareto front for NSGA-II Table with 250 iterations, 100 research agents, with four selected solutions out of a total of 175.

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 Figure7for the Pareto front and Figure8for 2D plots of cloak design and wave fields.

Fig. 10 .
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.

Table 6 .
Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak.n max =1000,Q = 200 .

Table 5 .
Comparison of a NSGA-II bi-objective method in protection and invisibility for a half-rectangular cloak.n max = 250, Q = 100.