Cross-Entropy Method for Design and Optimization of Pixelated Metasurfaces

Electromagnetic metasurfaces are planar two-dimensional metamaterials, typically of subwavelength thickness. Unit cell elements of different shapes have been widely explored, including electric and magnetic dipoles, patches, arbitrary geometries and pixelated surfaces. Although pixelated metasurfaces have a great advantage of geometric versatility, their design and analysis requires algorithmic approach. One of the techniques for their design is via evolutionary simulation-driven optimization. Since full-wave electromagnetic simulations are time-consuming, optimization methods with fast convergence properties are preferable. In this paper, we demonstrate the application of the cross-entropy optimization method to design of artiﬁcial magnetic conductors (AMCs) and thin printed phase shifters. Single-frequency AMCs at 10 GHz (X band) and dual-frequency AMCs at 8 and 12 GHz (X and Ku band) were produced that are more manufacturing-friendly, and thus cost effective, than previously reported AMCs. We also show that phase-shifting unit cells with transmission magnitudes over 0.9 (linear) can be designed using the proposed optimization technique. Other potential of these unit cells are in phase-correcting and beam-steering metasurfaces.


I. INTRODUCTION
E LECTROMAGNETIC surfaces, also called metasurfaces, are planar structures of finite thickness that are composed of subwavelength periodic or nonperiodic unit cells. A unit cell usually consists of one or more layers of dielectric materials with metallic patterns printed on their opposite sides. The prominent feature of such structures is a simultaneous control of the magnitude, phase and polarisation of electromagnetic (EM) waves. Interesting properties such as single-or dual-frequency narrowband or broadband antireflection, absorption and polarisation conversion [1], [2] can be achieved through manipulation of their reflection and transmission characteristics. The applications of EM surfaces in antenna and microwave engineering include absorbers [3], multiband spatial filters [4], conformal antenna radomes [5], reflectarrays [6] and transmitarrays [7], EM shieleding [8] and radar cross section reduction [9], to name a few.
Often, frequency-selective surfaces (FSSs), partially-reflective surfaces, high-impedance surfaces and artificial magnetic conductors (AMCs) are realised using metallic squares, loops, cross-type elements and their combinations [10]- [13]. The required reflection and transmission characteristics of an EM surface can be obtained by finding the appropriate geometrical parameters of the element in a unit cell. When designing either a band-pass or band-stop FSS, the choice of the element may be of utmost importance [10]. An alternative approach is to use a pixelated unit cell with discretised conductive layers and obtain the desired performance by optimizing the pixel pattern. The advantage of pixelated unit cells is the versatility of the geometry that is not limited to any canonical shape. Optimization of pixelated EM surfaces has been previously conducted by encoding a pattern of "metal-air" inclusions into a binary string. Using genetic algorithms (GAs) with their inherent binary representation of design variables were the first methods applied to such optimization. Using GAs, a variety of printed pixelated structures have been designed, such as a microwave absorber composed of multiple dielectric layers with FSS screens [3], a partially-reflective surface of an RCA for high-Q resonant cavity [14], periodic metamaterials [15] and a double-sided AMC [16]. Apart from GAs, other methods have also been implemented for the optimization of pixelated printed EM surfaces, such as particle swarm optimization, wind driven optimization and simulated annealing [17]- [19].
Recently, a new optimization method called the crossentropy (CE) method has been introduced to EM community [20]- [22] for direct optimization. It showed simple handling of continuous and mixed-variable variables, versatile sampling under design constraints and comparatively fast convergence. Therefore, it is worth investigating the performance of the CE method in application to the design of binary EM problems, such as pixelated metasurfaces. Previously, the possibility of using CE for binary EM problems has been explored in array synthesis [23] but no results have been demonstrated for simulation-driven optimization.
In this paper, we propose an optimization approach and implement the CE on the design of thin composite EM metasurfaces. Pixelated patterns printed on both sides of a thin dielectric material have been optimized targeting various frequency responses. By controlling the phase and magnitude of the EM wave reflected from both sides of the metasurface and the transmission properties, AMCs and planar phase shifters have been designed. This paper is organized as follows. In Section II, a method of analysis is described, along with a technique to improve the reliability of EM surface design. In Section III, an optimization strategy using the CE method is given in detail. Section IV-VI provide the results of optimized singlefrequency and dual-frequency AMCs and phase-correcting metasurfaces. Section VII concludes the paper.

II. METHOD OF ANALYSIS
In order to efficiently optimize an EM structure, an accurate simulation model with minimal computational burden has to be established. The method of analysis and some simplifications that have been introduced in this work in order to reduce the overall optimization time are described below.

A. PERIODIC ANALYSIS AND SYMMETRY PLANES
A two-dimensional periodic structure is created by translating a unit cell in x and y directions, as shown in Fig. 1. According to Floquet-Bloch theory, it is sufficient to analyse the fields in a single unit cell in order to predict the characteristics of the infinitely large periodic EM surface. Periodic boundary conditions should be applied around the unit cell to imitate the infinite dimensions of the structure. The waveguide simulation method in conjunction with timedomain analysis in CST MWS has been used to predict the reflection and transmission characteristics of all the cases presented in the following discussion. Illumination by a normal plane wave from two opposite ports with the E- field aligned with x-axis was considered, and perfect electric and perfect magnetic boundary conditions were assigned to model an ideal parallel-plate waveguide. Taking xz-plane as a plane of incidence, the applied plane wave is TE-polarized.
A square unit cell, shown in Fig. 1, consists of a dielectric material of thickness d with two pixelated patterns printed on its opposite sides. The printed surface on both sides of a unit cell have equal side length a and level of pixelation n (number of columns/rows). Each pattern (both top and bottom) is represented by a matrix of n 2 elements that can take value either "1" or "0", where "1" corresponds to "metal" and "0" to "air".
The design variables are expressed as: where k = 1 and 2 correspond to the top and bottom pattern, respectively. Therefore, a unit cell is represented by a binary vector X = [X 1 , X 2 ] of 2 * n 2 elements. The optimization problem is to find the binary combination in X that corresponds to the unit cell with the desired transmission and reflection characteristics.
To further reduce the computational cost of the optimization problem, we impose symmetry within the unit cell. As shown in Fig. 2, a quarter of a unit cell can be transformed to a full pattern using three types of symmetries, i.e. reflectional, rotational and translational. Although there are other possible ways to reduce the number of variables in a unit cell, a 90degree rotational symmetry not only reduces the number of design variables by a factor of four, but also ensures polarisation independence of the structure. Thus, in the designs  presented further, a 90-degree rotational symmetry in the unit cell will be applied.

B. VERTICES-REMOVAL TECHNIQUE
Pixelated surfaces are known for the fabrication issue related to the undesired interconnections of the diagonal pixels. There are two types of possible connections between square pixels: vertex to vertex and edge to edge. The problem occurs because in numerical analysis, pixels with a vertex junction are often treated as not electrically connected. In manufactured prototypes, these vertices can sometimes become physically connected, allowing current to flow between the pixels. This causes disagreement between predicted and measured results.
A number of techniques have been proposed to overcome with this issue. In [24], [25], a vertex breaking technique separating the vertices of diagonal pixels by the finite distance of 0.04 mm has been applied. The disadvantage of this technique is that geometry modification is required for every generated pattern, which with the numerous simulations required by the optimization, might substantially increase the computational time. A geometry refinement technique proposed in [26], [27] eliminates the designs with diagonal connections of square pixels by analysing the geometry of every solution. However, by disregarding many potential solutions, this approach considerably reduces the search space available for the optimizer. Another approach called the nonuniform overlapping scheme, which ensures electrical connection of diagonal pixels by overlapping them, has been proposed in [28] and [29]. The idea is to intentionally overlap the metallic pixels in the simulation model by shifting the diagonally adjacent pixels in the y direction. This ensures the existence of the electrical connection in both numerical analysis and the fabricated prototype. Again, this scheme introduces an additional step to the overall optimization flow.
In this paper, the described problem is eliminated by using octagonal pixels instead of rectangular pixels. This can be easily achieved by truncating each square pixel. This approach has the advantage of not adding any complementary steps to the optimization procedure, as the conversion is only done once even before the initiation of the optimization. The proposed solution comes at a cost of an increased mesh size of the structure in the numerical analysis. To be specific, the number of meshcells increases from 14 000 with square pixels to 22 500 with octagonal pixels due to the need for a finer resolution. However, due to the short evaluation time (approximately 60 sec for each unit cell), the overall optimization time is still acceptable.

III. OPTIMIZATION TECHNIQUE A. SAMPLING BINARY VARIABLES
The main principle of the CE method in search of the global optimal solution is to consecutively adapt the shape parameters of the sampling probability distributions [20], [22]. It is the nature of the design variables that guides the choice of the probability distribution used for sampling. For example, for optimization problems with design variables that have only two possible values (binary), a Bernoulli distribution can be used. Bernoulli distribution is a discrete probability distribution of a random variable that takes the value 1 with probability p and the value 0 with probability (1 − p). The probability mass function f (x; p) of a Bernoulli distribution over possible outcomes x is where x is either 1 or 0, meaning that f (x; p) = p when x = 1, and f (x; p) = (1 − p) when x = 0. Therefore, if encoding a pixelated unit cell into a binary string, the CE method with the Bernoulli sampling distribution can be used for its optimization. Using a 90-degree rotational symmetry in a unit cell with n = 6, the quarter of the pattern on top of the dielectric layer is represented by a probability distribution matrix P T : and the quarter of the pattern on the bottom of the dielectric layer is represented by P B : where the superscripts T and B refer to the top and bottom of the unit cell, respectively. The probability p ijk is the probability that the pixel located at row i = 1, ..., n/2, column j = 1, ..., n/2 and side k = {1; 2} is filled with metal. The quarter of a unit cell is, therefore, sampled using the sampling distribution parameter matrix P = P T P B with number of elements D = 2 × (n/2) 2 = 18 for n = 6, where n is the number of pixels in each row/column of a VOLUME 4, 2016 unit cell. The number of parameters D in the quarter matrix P is the dimensionality of the optimization problem. Fig. 3 illustrates that the difficulty of binary optimization problems increases exponentially with the number of dimensions. If D = 1, there are only two possible outcomes (straight line), that represent either "air" or "metal". If D = 2, the number of outcomes equals the number of vertices of a rectangle, i.e., 2 2 = 4. For D = 3, the number of outcomes is 2 3 = 8.
In the examples in this paper n = 6, thus, the binary optimization problem with 18 design variables has 2 D = 2 18 = 262 144 possible outcomes.

B. CE ALGORITHM FOR OPTIMIZATION OF BINARY PROBLEMS
The main idea of the CE optimizer is the minimization of the cross-entropy between two probability distributions: an empirical distribution describing the current elite subpopulation and a sampling distribution which is used to sample a new population. At each iteration, two important steps must be taken to implement this idea: • Generate a random population of solutions x t with N pop candidates from the sampling distribution f (x t ; p t ), where p t are the distributional parameters at the t-th iteration, and choose the N el best-performing candidates for the elite subpopulation by evaluating the fitness function. • Update the shape parameters p t+1 of f (x t+1 ; p t+1 ) by minimising its cross-entropy with the empirical distribution g(x t ; w t ) describing the current elite solutions using the maximum likelihood estimation of w t .
The flowchart of the proposed optimization scheme for the design of pixelated printed surfaces by the CE method with a Bernoulli sampling distribution is given in Fig. 4. The optimization begins by setting the parameters of the CE method, such as population size N pop , proportion of elite subpopulation ρ = N el /N pop and smoothing parameter α S . Also, it is required to set the parameters of the optimization problem, such as the number of rows N p and the number of columns M p in the probability matrices P T and P B . The initial unit cell patterns X 0 are generated by sampling from the uniform probability mass function P 0 with all probabilities equal 0.5. At each iteration t, the new population is generated by independently setting x t ijk = 1 with probability p and x t ijk = 0 with probability (1 − p), where "1" corresponds to perfect electric conductor and "0" corresponds to vacuum. The next step is to obtain the performance parameters via EM simulation and evaluate the fitness functions for each candidate. If the fitness function is to be maximized, the results are sorted in descending order, and if minimized, in ascending order. In either case, the first ρ% = N el ·100%/N pop constitute the elite subpopulation. These N el candidates are used to generate new probability mass functions P (t+1) , to increase the probability that the best-performing candidates occur in the next generation. This is achieved by maximum likelihood estimation, which for the Bernoulli distribution is simply the mean of the best-performing candidates [30]: An optional smoothing procedure described in [31] has been applied. The termination condition has been defined as the diversity of the elite subpopulation. If the variance of the elite candidates equals to or goes below the threshold, the optimization stops. For binary variables, the threshold is δ=0, which implies that optimization stops when the elite subpopulation has converged to the same solution.

IV. SINGLE-FREQUENCY AMC SURFACES
The first AMC design was optimized for X-band operation at 10 GHz. The side length of the unit cell was a=8 mm with each octagonal pixel having the height h=1.3 mm with one side of an octagon being 0.75h and another 0.15h. The dielectric material Rogers RO4003 with ε r = 3.55 and thickness d=1.58 mm was used to account for losses. The target phase and magnitude (on a linear scale) of the reflection and transmission coefficients at the required frequency are the following: The optimization goal was to minimise the objective function F.F., which is the sum of three terms: where FIGURE 4. The flowchart for optimization of binary problems using the CE method.
As the transmission coefficient magnitude |S 12 | of 0.3 is sufficient, all the designs with |S 12 | 0.3 have F 3 =0. The parameters of the CE method were set to the following values: population size N pop =30, elite subpopulation size N el =10 and smoothing parameter α S =0.6. The time required for a single simulation is between 30 and 60 s using a PC with Intel Core i7-4790 processor and 32 GB of memory. Five optimization runs were conducted to keep the balance between the statistical certainty of the results and the computational cost.
The best obtained solution has F.F.=6.8 and is represented by the following solution matrix: The top and bottom patterns of the surface with 3 × 3 unit cells, obtained after CE optimization for the goal given in (5), are shown in Fig. 5. To construct the geometry of a periodic EM surface from the optimized unit cell, the image theory has been applied. The mirror image of a current element flowing parallel (perpendicular) to the PEC boundary should be formed in the opposite (same) direction. The image of a current element flowing parallel (perpendicular) to a PMC boundary is formed in the same (opposite) direction [32]. The transmission and reflection results are given in Fig. 6, where S 11 and S 22 represent reflections from the top and the bottom surface, respectively. As required, ∠S 11 = 0, |S 12 | = 0.25, which is −12 dB, and ∠S 22 = 174 • at 10 GHz. Therefore, the top surface of the designed EM structure behaves as an AMC at the specified frequency. An in-phase reflection bandwidth, which is commonly defined as the frequency band with the phase of the reflection coefficient being within ±45 • , is from 9.41 GHz to 10.59 GHz.
It is also worth noting that, to obtain this AMC, a complete metallic backing was not required. As shown in Table 1, in comparison to the pixelated surfaces obtained through GA in [15], [16], the thickness of this AMC is merely the same but the level of discretisation is significantly lower, 6 × 6 pixels as opposed to 16 × 16, which is advantageous for fabrication purposes. The side length of the unit cell in [15] is 0.1λ 0 at 5 GHz, which is 6.6 mm, and thus the dimensions of each pixel are 0.4 × 0.4 mm. The table also shows whether the unit cells are polarisation-dependent (PD). Neither of the designs from Table 1 have been fabricated due to the focus of the work being on the optimization methodologies. However, practical considerations should be taken into account. The side length a = 8 mm and coarser discretisation in the unit cell have been chosen due to anticipated fabrication limits. The thickness t = 1.58 mm has been considered because the AMC surface of 100λ 0 × 100λ 0 has to be rigid for the measurement purposes. To demonstrate that AMC design with a = 0.1λ 0 is also VOLUME 4, 2016 possible, a unit cell with a = 3 mm and 36 pixels (each pixel having b=0.5 mm) on each side has been optimized. The 3×3 surface is shown in Fig. 7, and its reflection and transmission characteristics are given in Fig. 8. The desired specifications have been satisfied, as |S 12 | = 0.14, ∠S 11 = 2.2 • and S 22 = 175.5 • at 10 GHz. The best obtained solution is represented by the following solution matrix: The optimization of the single-frequency AMC only took from 8 to 15 iterations, which is three times less computational effort than that required by GA in [15]. The convergence curve of the optimization of the single-frequency AMC with a = 0.1λ 0 is shown in Fig. 9.  A more detailed visualisation of the fitness function improvement can be observed in Fig. 10, which compares the fitness function values calculated by (6) for every candidate at the first and the last iteration of the optimization for a = 0.27λ 0 . The candidates are sorted in descending order of their final F.F. value. Initial score distribution is random with the highest F.F. = 890. A significant overall reduction can be seen in the fitness function that indicates convergence. At the last iteration, the lowest F.F. = 6.8.

V. DUAL-FREQUENCY AMC SURFACES
The second AMC design, with the unit cell dimensions a = 8 mm and d = 3.16 mm, was optimized to operate at 8 GHz and 12 GHz. The target characteristics were defined as in (5) at these two frequencies. The optimization goal was to minimise the objective function where F.F. 8GHz and F.F. 12GHz are calculated using (6) - (9). Out of five consecutive optimization runs, each with population size N pop =30, elite subpopulation size N el =10 and smoothing parameter α S =0.1, the two best solutions had almost the same F.F. value and design parameters. The solution shown in Fig. 11 with F.F.=41.3 has the following solution matrix: The second solution has F.F. = 35.2, and the only difference in the solution matrix is a single bit in the bottom pattern, which has all "1s". It suggests that for dual-band performance, a complete or nearly complete metallic backing is required.
The transmission and reflection results are given in Fig. 12. Phase ∠S 22 is 177 • and 175 • , and |S 12 | is -32 dB and - 40 dB at 8 and 12 GHz, respectively. Phase ∠S 11 is -25 • at 8 GHz, which is within the in-phase reflection bandwidth (±45 • ), and it passes through 0 • at 12 GHz. Therefore, the optimized surface behaves as an AMC at both target frequencies of 8 and 12 GHz. In comparison to the pixelated surface obtained through GA in [15], this AMC is electrically thinner, t = 0.09λ r versus 0.13λ r , and has a larger side length, a = 0.27λ 0 versus 0.11λ 0 .
As in the previous case, a comparison given in Table 2 shows that the thickness and the side length of this AMC are merely the same as in [16] but the level of discretisation of the optimized unit cell is significantly lower due to the fabrication considerations. The optimization run that produced the design in Fig. 11 continued for 47 iterations. Its convergence curve is given in Fig. 13. The best result was obtained after the 27th iteration, and the optimization was terminated after not producing any better solution for as long as 20 iterations.

VI. THIN PHASE-CORRECTING METASURFACES
Other applications of EM surfaces are found in the design of low-profile planar phase correcting structures and beamsteering metasurfaces [33], [34]. Using the proposed optimization methodology by means of the CE method, we designed phase-shifting structures made of a single-dielectric layer with pixelated surfaces printed on the two sides. The unit cells were optimized to provide a required phase shift without introducing significant losses to the transmitted  wave. Eleven optimization problems have been formulated to obtain phase shifters covering a range of phases from 180 • to 360 • with 10 • to 20 • steps. The same unit cell as in the design of AMCs, with 36 pixels on each side, has been used for realisation of the phase shifters. Motivated by practical considerations for operation at 20 GHz (λ 0 = 15 mm), the side length of the unit cell has been set to a = 0.3λ 0 = 5 mm. In the simulation model, the pixelated patterns are placed on a commerciallyavailable dielectric material Taconic TLY-5 with ε r = 2.2, t = 1.58 mm and loss tangent of 0.0009.
The optimization goal for all phase-shift values has been defined as: where VOLUME 4, 2016 and The coefficients in (11) have been chosen to reflect the relative importance of the phase and magnitude objectives and to compensate for the difference in their units.  Table 3, where the target transmission-phase ∠S 12 values and the ones obtained after optimization can be compared. For all solutions, the magnitude of the transmission phase |S 12 | is 0.9 (linear scale), which ensures transmission losses of less than -1 dB. The largest difference between the desired and obtained phase shift is 5.67 • , which is assumed to introduce only minor errors in application to phase correction. Fig. 14 demonstrates the phase shifts provided by the optimized designs given in Table 3, and the transmission phase and magnitude of a unit cell with a 200 • phase shift is shown in Fig. 15. It can be seen that a thin doublesided single-layer metasurface is capable of providing a 180 • phase-shifting range at 20 GHz. 19 19 .

VII. CONCLUSION
Optimization of pixelated metasurfaces is a binary multidimensional problem. Due to its fast convergence rate and inherent versatility in handling variables, we applied the CE method with a Bernoulli sampling distribution to optimize pixelated metasurfaces. Application of the Floquet-Bloch theory allowed us to minimize the design space to a single unit cell with periodic boundary conditions enforced in CST MWS. Rotational symmetry was assumed in the unit cell for polarisation independence, which allowed further reduction in the design space by a factor of four. Each side of the unit cell was discretised into 36 pixels, resulting in the total of 18 binary design variables. To eliminate the issue of poor correlation between simulated and measured results, we implemented an octagonal pixel shape. Two 10 GHz-AMCs were produced with a unit-cell side length of 0.1 λ 0 and 0.27λ 0 that have phase response of a perfect magnetic conductor and transmission magnitude −17 dB and −12 dB, respectively. Dual-frequency AMCs for 8 and 12 GHz were designed using the same method by changing the objective function. The results with the CE method were obtained after only around 500-1000 function evaluations, which is much shorter than typical thousands of evaluations required by genetic algorithms. The optimized phase-shifting unit cells with pixelated patterns are half the thickness of the metasurface with square patches in [34] and can be used to improve the beam-steering solution by decreasing the profile, weight and complexity of the turning metasurfaces. The proposed optimization methodology and the obtained unit cells can also be used in the design of transmitarrays [7]. A more precise phaseshifting element might be obtained by optimizing a unit cell with an increased level of discretisation. Other phase values can be obtained by defining the desired ∠S 12,obj in the optimization goal.