IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO.

, 2021 1

A Combined Machine-Learning /
Optimization-Based Approach for Inverse Design of
Nonuniform Bianisotropic Metasurfaces
Parinaz Naseri, Stewart Pearson, Zhengzheng Wang, and Sean V. Hum
Submitted to IEEE Transactions on Antennas & Propagation,
Machine Learning in Antenna Design, Modeling, and Measurements Special Issue
arXiv:2105.14133v1 [[Link]] 28 May 2021

Abstract—Electromagnetic metasurface design based on far- Kem , is also needed for useful wave transformations [3]–
field constraints without the complete knowledge of the fields [4]. Since these macroscopic parameters describe the elec-
on both sides of the metasurface is typically a time consuming tromagnetic properties of the same EMMS, they are equiv-
and iterative process, which relies heavily on heuristics and
ad hoc methods. This paper proposes an end-to-end systematic alent and interchangeable. Here, we use the surface electric
and efficient approach where the designer inputs high-level far- impedance (Z̄¯se ), surface magnetic admittance (Ȳ¯sm ), and
field constraints such as nulls, sidelobe levels, and main beam electro-magnetic coupling coefficients (K̄ ¯ ) to describe the
em
level(s); and a 3-layer nonuniform passive, lossless, omega-type surface properties of the EMMS.
bianisotropic electromagnetic metasurface design to satisfy them
The macroscopic and microscopic steps of the inverse
is returned. The surface parameters to realize the far-field criteria
are found using the alternating direction method of multipliers problem are shown with steps 1 and 2 in Figure 1, respectively.
on a homogenized model derived from the method of moments. In practical and general problems where the description of the
This model incorporates edge effects of the finite surface and field on one side is incomplete or high level constraints on
mutual coupling in the inhomogenous impedance sheet. Opti- the desired radiation pattern are given, the surface properties
mization through the physical unit cell space integrated with
need to be obtained without the knowledge of the tangential
machine learning-based surrogate models is used to realize the
desired surface parameters from physical meta-atom (or unit electric and magnetic fields on one side of the surface. Hence,
cell) designs. Two passive lossless examples with different feeding the first step of the inverse problem can no longer be solved
systems and far-field constraints are shown to demonstrate the analytically and requires optimization to search for the surface
effectiveness of this method. parameters [8]–[13]. Moreover, the design of the right unit
Index Terms—Electromagnetic metasurfaces, inverse design, cell structure based on the macroscopic surface parameters is
metasurface synthesis, machine learning, deep learning, surrogate a one-to-many mapping where multiple solutions may exist in
models, optimization. a high-dimensional design space. This makes the search for
the right design difficult and has resulted in this step being
I. I NTRODUCTION implemented by mostly ad hoc or empirical approaches in
the literature. So far, this approach relies on time-consuming

E LECTROMAGNETIC metasurfaces (EMMSs) are thin

uniform or nonuniform 2D arrangements of sub-
wavelength unit cells that are composed of patterned metallic
and resource-demanding cycles of optimization and full-wave
simulations, that has been the bottleneck of the EMMS design
process.
scatterers and/or dielectric substrates. These special surfaces
provide the ability to manipulate electromagnetic waves in Far-Field Constraints: Step 1 Step 2
extraordinary ways. Examples include spectrum filtering, wave • Main beam direction(s) 𝑍𝑠𝑒 𝑥 , 𝑌𝑠𝑚 𝑥 , Physical Unit cell
manipulation, and polarization conversion [1]. The inverse • Sidelobe levels 𝑎𝑛𝑑 𝐾𝑒𝑚 (𝑥) Structures
• Nulls locations
design of an EMMS generally involves two steps: mapping
the constraints on the fields of one side of the surface
while knowing the fields on the other side of the EMMS to Figure 1. The steps to inverse design a metasurface based on the
macroscopic properties; and then, achieving these properties constraints on the desired radiated far-field, based on the the surface
impedance approach..
microscopically using physical unit cells. The EMMS macro-
scopic electromagnetic properties can be described in terms of Interesting power-conserving wave transformations are pos-
the scattering parameters, surface susceptibilities [2], surface sible with nonuniform surfaces where the surface parameters
impedance/admittance [3]–[5], or equivalent permittivity and are a function of position on the EMMS surface. These
permeability tensors [6]–[7]. It is shown that bianisotropy, parameters translate to a nonuniform array of unit cells that
which can be modeled by a magneto-electric coupling term are composed of at least three layers of asymmetric patterned
metallic scatterers [3] in the physical solution domain. It is
P. Naseri, S. Pearson, Z. Wang, and S. V. Hum are with the Ed- worth noting that the surface properties of each unit cell
ward S. Rogers Sr. Department of Electrical and Computer Engineering,
10 King’s College Road, toronto, Ontario, Canada, M5S3G4, email: pari- are characterized by simulating them in periodic boundary
[Link]@[Link] conditions. Therefore, when placed in a nonuniform array, the
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 2

unit cells’ effective properties can vary from their designed

values. This change is attributed to mutual coupling between Far-field
each unit cell with its neighbors. Neglecting to account for constraints
this coupling results in errors in the wave transformation. The
designer is then required to simulate and tune each of the ADMM-based convex
many variables such as scatterers’ dimensions to obtain the macroscopic optimization
desired response. This ad hoc approach requires many cycles using MoM model
of optimization through expensive full-wave simulations and
is very time- and resource-demanding. For transverse electric
(TE) incident waves, the use of vertical baffles to implement Optimized
perfect conductor walls between neighbor unit cells has been (𝑍𝑠𝑒 , 𝑌𝑠𝑚 , 𝐾𝑒𝑚 ) or [𝑆]
proposed to suppress the mutual coupling [14], but it compli-
cates the fabrication and only works for TE excitation.
Exploration in the
The second step involves optimizing the unit cell struc-
physical solution
ture based on the required macroscopic properties. Designers ML-Based
space
mostly rely on empirical approaches including many cycles of PSO Microscopic
optimization and simulations to obtain the right choices for the Optimization
unit cells’ scatterer shapes, dimensions, substrate thicknesses
ML-based
and permittivity(ies). Searching for the optimized structure
is a challenging problem since the solution space includes surrogate models
a wide variety of scatterer shapes and corresponding feature
dimensions. This solution space might include many local
Optimized
minima, which make the optimization over this space diffi-
physical EMMS
cult. Moreover, to find the right choice, each candidate must
be simulated and evaluated against the required properties.
Figure 2. Proposed systematic approach for inverse design of a
Inspired by the revolution that data-driven machine-learning
nonuniform bianisotropic metasurface based on the constraints on
methods have made in material informatics such as discovery the radiated far-field pattern.
of new quantum materials, pharmaceuticals, and other com-
pounds [15], deep machine learning [16]–[21] and statistical
learning [22] methods can help to build surrogate models that optimization step, the alternating direction method of multi-
can provide fast predictions of the properties of each unit pliers (ADMM)-based convex optimization using the periodic
cell. Moreover, several machine-learning techniques have been method of moments (MoM) is implemented to deterministi-
proposed to tackle the challenges of the second step [23]–[33]. cally obtain the surface parameters. Then, in the microscopic
Some of these proposed methods deal with the inverse design optimization step, the desired surface parameters are realized
of a uniform EMMS where the impact of mutual coupling is with physical unit cells using particle swarm optimization
less of an issue [23]–[29], while the rest optimize over a simple (PSO) integrated with machine learning (ML) surrogate mod-
solution space that is composed of only one scatterer shape els. Mutual coupling between the unit cells and edge effects
[30]–[33]. It is worth noting that dielectric optical EMMSs that are of main concerns in the design of such finite sur-
can be designed using global optimization methods due to the faces are taken into account to minimize time-consuming and
analytical relation between the scatterers’ properties and the resource-demanding simulation-based optimization of the full
EMMS’s scattering parameters [34]. However, due to the lack array as well as avoiding using more manufacturing-intensive
of such relations in EMMSs composed of metallic scatterers, structures such as baffles (vias) [14].
the inverse design of a heterogeneous bianisotropic EMMS An electromagnetic-aware efficient method is used to gen-
is more challenging. Nonetheless, solving this problem more erate training data for ML surrogate models. These models
efficiently has led to the proposal of systematic approaches that successfully and efficiently predict the magnitude as well as
provide both the optimized surface properties and the actual the phase of the scattering properties and replace the time-
physical unit cells to synthesize a stack of two EMMSs that consuming full-wave simulations that are required for the
generates two pencil beams at two frequency bands [35]. optimization over the physical solution space. A new method
Here, to expand the applications of the EMMSs, arbi- to represent the 3-layer unit cells as the inputs of the surrogate
trary far-field constraints are given as the inputs and the models is implemented that reflects both the categorical and
macroscopic and microscopic steps of the inverse problem continuous nature of the solution space composed of different
of the EMMSs are combined together into a unified design types of scatterers with various dimensions that are separated
approach that helps to implement a nonuniform bianisotropic with different substrate thicknesses.
metasurface to satisfy them. Figure 2 shows the flowchart This paper is organized as follows. Section II details the
of the proposed systematic approach. This approach allows macroscopic step to convert the high-level far-field constraints
a designer to place practical constraints on the radiated far- to surface parameters. In Section III, we explain the different
zone fields emanating from the physical implementation of the parts of the microscopic step involving the ML-based surro-
EMMS in an efficient and systematic way. In the macroscopic gate models integrated in the particle swarm optimization to
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 3

obtain the optimized physical EMMS. Two design examples Since Kem = Kme , we will simply use Kem to denote the
to demonstrate the effectiveness of the proposed approach are bianisotropic coupling terms. The details for the derivation
presented in Section IV. Conclusions are drawn in Section V. of passive, lossless, and reciprocal omega-type bianisotropic
GSTCs can be found in the work by Ataloglou et al [37].
II. FAR -F IELD C ONSTRAINTS TO S URFACE PARAMETERS Because we consider the EMMS within the model to have
Zse , Ysm , AND Kem zero thickness, we can substitute the average fields in (1) with
The first step of the inverse design process can be formu- the (known) incident and scattered fields. The scattered electric
lated as an optimization problem that searches for the surface field due to an electric current density along a thin strip on
properties to satisfy local power conservation as well as the the y-axis is given by
far-field constraints. After formulated as such, the problem can
be solved deterministically. The details of this macroscopic Z w/2
~ xscat (~ ωµ0 (2)
optimization procedure have been presented [36]. The first step E ρ) = − J~x (y 0 )H0 (k|y − y 0 |)dy 0 (3)
4 −w/2
is to describe the EMMS as a fully homogenized structure
using an integral equation approach implemented using the
method of moments (MoM). where w is the width of the thin strip and ρ ~ = y ŷ is a position
vector of the field point along the y-axis. Due to the fact that
M~ y is y-directed, it will not contribute to the scattered electric
A. EMMS Model from the Method of Moments field in the x-direction. We can similarly derive the y-directed
Starting with the generalized sheet transition conditions scattered magnetic field resulting from the magnetic current
(GSTCs) for a bianisotropic EMMS, density as
1 ~ ~ t,2 ) = Z̄¯se · J~s − K̄
¯ · (n̂ × M ~ s)
(Et,1 + E em (1a)
2 ~ scat (~ 1 ∂2
1 ~ H ρ) = − (k 2 + 2 )
(Ht,1 + H~ t,2 ) = Ȳ¯sm · M ¯ · (n̂ × J~ ).
~ s + K̄me s (1b) 4ωµ0 ∂y
2 Z w/2 (4)
These equations relate the tangential electric and magnetic ~ y (y 0 )H (2) (k|y − y 0 |)dy 0 .
M 0
~ t,1 , H
~ t,1 , to those on −w/2
fields on one side of the EMMS, E
~ ~
the other side, Et,2 , Ht,2 , through the surface parameters and
the surface currents. The surface parameters for this selected Similar to the scattered electric field, J~x is x-directed and will
representation are the surface electric impedance (Z̄¯se ), surface not contribute to the scattered magnetic field in the y-direction.
magnetic admittance (Ȳ¯sm ), and magneto-electric and electro- Using the MoM with pulse basis functions and point match-
magnetic coupling coefficients (K̄ ¯ , K̄ ¯ ). The surface cur-
em me ing along the EMMS we can create linear system of equations
rents are the electric surface current density (J~s ) and magnetic [36],
surface current density (M ~ s ).

Einc = −Escat + j[Xse ]Ie − [Kem ]Im (5a)

inc scat m e
H = −H + j[Bsm ]I + [Kem ]I (5b)
Escat = −[Ze ]Ie (5c)
scat m
H = −[Zm ]I , (5d)

Figure 3. Configuration of the homogenized EMMS model. where Ie and Im are the basis coefficients for the electric and
magnetic surface current, respectively. The elements of [Ze ]
and [Zm ] are
In this paper, we will use 2D TE-polarized examples. A
diagram illustrating this configuration is shown in Figure 3, Z v∆y
where the EMMS is periodic in the x-direction and varies ωµ0 (2)
zeuv = H0 (k|yu − y 0 |)dy 0 (6a)
nonuniformly in the y-direction. The EMMS will also be 4 (v−1)∆y
omega-type bianisotopic. As a result, the surface parameters Z v∆y
k2 (2)
in (1) reduce to scalars as a function of position. In addition, uv
zm = H2 (k|yu − y 0 |)
we will only consider purely passive, lossless, and reciprocal 8ωµ0 (v−1)∆y (6b)
(2) 0 0
EMMSs. This is because they are much easier to realize in + H0 (k|yu − y |)dy ,
practice. As a result,
Re(Zse ) = 0 (2a) where u and v are row and column indices, respectively.
Re(Ysm ) = 0 (2b) In order to impose far-field constraints, we also need an
equation to transform the discretized current densities and
Kem = Kme (2c)
corresponding current coefficients Ie and Im to the far-field.
Im(Kem ) = 0. (2d) We first examine the continuous form of the far-zone electric
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 4

field at a distance ρ resulting from electric and magnetic where f0 and fi are some objective and constraint functions
surface currents on the metasurface, for a certain design goal. The optimization variables in this
−jkρ
r case are the surface electric and magnetic current coefficients
~ fe f (θ) = − e √ ωµ0
E
2 jπ
e 4 Ie and Im along with the surface electric reactance [Xse ],
ρ 4 πk
Z w/2 (7a) magnetic susceptance [Bsm ], and magneto-electric coupling
~ 0 jky 0 sin θ 0 [Kem ]. Because we only optimize for passive and lossless
Jx (y )e dy
−w/2 surface parameters, our solution is passive and lossless by
−jkρ
r construction.
~ m (θ) = − e √ ωε0 η0
E
2 jπ
e 4 This optimization problem is non-convex due to the biaffine
ff
ρ 4 πk
Z w/2 (7b) [Xse ]Ie − [Kem ]Im and [Bsm ]Im + [Kem ]Ie terms. Biaffine
~ 0 jky 0 sin θ 0 terms are created when two variables are multiplied together.
My (y )e cos θdy .
−w/2 Unfortunately, non-convex problems are notoriously hard to
where θ is illustrated in Figure 3. In general, macroscopic solve in polynomial time. Therefore, we cannot simply opti-
optimization is more concerned with relative magnitudes of the mize for the surface currents to satisfy the objectives because
far-zone pattern rather than the absolute value at a certain dis- that would likely yield active and/or lossy EMMSs. Optimizing
tance. As a result, we focus on the far-zone fields irrespective the surface currents and parameters together yields a passive
of distance and frequency E ~ 0,f f . Now, using our discretized and lossless EMMS by construction.
e m
surface currents I and I , we can rewrite (7) as a linear To solve this non-convex problem, we will relax it us-
system of equations over M angular far-zone samples and N ing ADMM [38]. ADMM is an algorithm that alternatingly
spatial samples along the surface as minimizes the augmented Lagrangian of (10) with respect to
s ! surface currents and parameters with each iteration. This is
e e −jkρ λ well suited to solving biaffine problems. Now that we can
Ef f =E0,f f e (8a)
ρ solve EMMS optimization problems in the form of (10), we
can devise far-field objective and constraint functions [36].
Ee0,f f =[Ge ]Ie (8b)
s ! 1) Main Beam Level: To try to force a beam to achieve
−jkρ λ a certain level (or get as close as possible), we can add a
Em m
f f =E0,f f e (8c) `2 -norm minimization term to the objective. For example, to
ρ
m m
force a beam at θ0 to a level M Blevel the objective function
Em
0,f f =[G ]I , (8d) can take the form
where [Ge ] ∈ CM ×N and [Gm ] ∈ CM ×N are matrices fM B (θ0 ) =k[Ge ](θ0 )Ie + [Gm ](θ0 )Im
expanding the discretized electric and magnetic surface cur- (11)
+ Einc 2
f f (θ0 ) − M B level k2 ,
rent coefficients over their pulse basis functions and then
transforming them to samples of far-zone electric field. These where Einc
f f (θ0 ) is the incident field across the extent of the
matrices have elements surface in the far-field.
r
ωµ0 2 j π jkyv sin(θu ) 2) Maximum Sidelobe Level: Enforcing a maximum per-
e,uv
g =− e 4e ∆y (9a) missible sidelobe over a certain set of angles can be done
4 r πkλ
with the inequality constraint
ω0 η0 2 j π jkyv sin(θu )
g m,uv = − e 4e cos(θu )∆y, (9b)
4 πkλ |Ge (SL)Ie + Gm (SL)Im + Einc
f f (SL)|≤ τ + slackSL , (12)
where u and v are row and column indices, respectively. These
where we force the modulus of the total field over the sidelobe
elements were solved for by expanding J~x and M ~ y in (7) over
region SL to be less than a certain sidelobe level τ . The
their pulse basis functions and integrating with the midpoint
modulus is permissible here because it remains convex when
rule.
used outside of an `2 -norm. We have included a slack variable
slackSL so the following term must be added to the objective:
B. Optimization Formulation
We can now frame the EMMS macroscopic design as an fSL = kslackSL k22 (13)
optimization problem This is required in order to make the inequality active.
minimize f0 (Ie , Im , [Xse ], [Bsm ], [Kem ]) (10a) 3) Surface Current Smoothness: In order to allow for more
Ie ,Im ,[Xse ],[Bsm ],[Kem ]
feasible EMMS designs, we add a surface current coefficient
Einc = [Ze ]Ie + [Xse ]Ie Ie and Im smoothness constraint. This can be done with the
subject to (10b)
− [Kem ]Im inequality constraints
Hinc = [Zm ]Im + [Bsm ]Im |[D]Ie |≤ Dmax
e
+ slackDe (14)
(10c)
+ [Kem ]Ie |[D]Im |≤ Dmax
m
+ slackDm , (15)
fi (Ie , Im , [Xse ], [Bsm ], [Kem ]) ≤ 0
where [D] is the discrete second derivative matrix. Similar
i = 1, ..., n, to (13), we have a function fD in the objective function to
(10d) minimize the `2 -norm of the slack variables.
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 5

4) Complete Formulation: We can assemble a complete single split ring [Figure 4 (c)], complete ring [Figure 4 (d)],
formulation if we fill in the objective function and inequality dog bone [Figure 4 (e)]. The primitives shown in Figure 4
constraints of (10) with those described above. This can be (a)-(e) are used for the scatterers at the air-dielectric interface,
written as i.e. the top and bottom layers. Often alternating inductive
and capacitive behavior is required in bianisotropic EMMSs
αM B fM B (M B) + αN U fN U (N U ) utilizing an odd number of layers. Therefore, we extend the
minimize range of scatterer geometries in the middle layer to include the
Ie ,Im ,[Xse ],[Bsm ],[Kem ], + αSL fSL + αD fD
slackSL ,slackDe ,slackDm complementary Jerusalem cross and complementary rectangu-
(16a) lar patch in Figure 4 (f)-(g) as well. The values and ranges of
Einc = [Ze ]Ie the different dimensions of each primitive are listed in Table
subject to e m (16b) I. Based on the range of each parameter, the numbesr of each
+ [Xse ]I − [Kem ]I primitive are also listed in Table I.
β(Hinc = [Zm ]Im
(16c) 𝑦 𝑦 𝑦
+ [Bsm ]Im + [Kem ]Ie ) 𝑤𝐽𝐶,𝑗𝑥 𝑙 𝑤𝑆𝑆𝑅
𝑃,𝑥
𝑙𝐽𝐶,𝑗𝑥
|[Ge ](SL)Ie + [Gm ](SL)Im 𝑙𝐽𝐶,𝑥
𝑟𝑆𝑆𝑅
𝑥 𝑥 𝑥
+ Einc 𝑔𝑆𝑆𝑅
f f (SL)|≤ τ + slackSL
𝑤𝐽𝐶,𝑥 𝑙𝑃,𝑦
(16d)
e e
|[D]I |≤ Dmax + slackDe (16e)
(a) (b) (c)
m m
|[D]I |≤ Dmax + slackDm , (16f) 𝑦 𝑦
𝑤𝐶𝑅 𝑙𝐷𝐵,𝑗

where M B ∈ [0◦ , 360◦ ] represents the set of main beam 𝑟𝐶𝑅 𝑙𝐷𝐵
𝑥 𝑥
angles, N U ∈ [0◦ , 360◦ ] represents the set of null angles, 𝑤𝐷𝐵
SL ∈ [0◦ , 360◦ ] is the set of angles comprising the side- 𝑤𝐷𝐵
lobe region, and αM B , αN U LL , αSL , αD are predetermined
weights for different terms in the objective function. Due (d) (e)
to different magnitudes of the electric and magnetic current 𝑦 𝑦
MoM equations, a scaling term is needed. The β term in 𝑤𝐽𝐶,𝑗𝑦
𝑤𝐽𝐶,𝑦 𝑙𝑃,𝑥
(25c) is used to scale the magnetic current MoM equation.
This is because when the augmented Lagrangian is formed, 𝑥 𝑥
the equality constraints compete for minimization. Here, β is 𝑙𝐽𝐶,𝑦 𝑙𝑃,𝑦
𝑙𝐽𝐶,𝑗𝑦
experimentally determined to be 1000. The weights are used if
different portions of the optimization should be stressed more. (f) (g)
For example, if the optimizer is having trouble satisfying the Metal Dielectric (𝜀𝑟 = 2.2)
sidelobe level constraint, αSL could be increased relative to
the other weights. This is done on an experimental basis. Figure 4. Primitives used for training: (a) Jerusalem cross (JC), (b)
rectangular patch (RP), (c) single split ring (SSR), (d) complete
ring (CR), (e) dog bone (DB), (f) complementary Jerusalem cross
III. F ROM OPTIMIZED S URFACE P ROPERTIES TO THE (compJC), and (g) complementary rectangular patch (compRP).
EMMS P HYSICAL S TRUCTURE
The primitives are simulated from 1 GHz to 19 GHz with
To realize the optimized surface parameters that are obtained unit cell period of 5.3 mm. Periodic boundary conditions on
in Section II with actual physical EMMS, we consider the x- and y-sides are enforced. We numerically determine the
solution space of 3-layer bianisotropic unit cells. This solution generalized scattering matrix (GSM) containing not only the
space is composed of diverse scatterers with various dimen- fundamental but also 10 other higher order modes, so that
sions that are separated with different substrate thicknesses and we can employ the cascading technique [39] to determine
permittivity(ies). To obtain the optimum structure efficiently, the 3-layer unit cell response [29]. This number of higher-
we exploit and explore this solution space using a global order modes provides sufficiently accurate results to capture
optimizer integrated with machine-learning surrogate models. the interlayer coupling for the specified frequency range and
The surrogate models are trained on a limited set of training unit cell period. Each scatterer is translated to meshes using
data to predict the unit cells’ properties accurately and expedite Rao-Wilton-Glisson (RWG) basis functions and fed to our in-
the microscopic optimization process. house periodic MoM-based simulation tool.
The GSM calculation for each primitive is accelerated
A. Generation of the Training Data using a model-based parameter estimation (MBPE) method
The training data set is composed of 3-layer unit cells with for the periodic MoM [40]. This interpolation method pro-
a specific scatterer on each layer. Based on experience to vides accurate computation of the scatterers’ GSMs over a
curate the training data set, we choose certain scatterer shapes, broad frequency band by performing the calculation only at
called primitives, here. The shape of the primitives include the a few random frequencies in the band. Patch-like scatterers
Jerusalem cross [Figure 4 (a)], rectangular patch [Figure 4 (b)], can be simulated using MBPE by solving MoM with the
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 6

TABLE I. Dimensions of the Primitives in Figure 4. varying dimensions of the scatterers, as shown in Figure 4,
Primitive Parameter Value (mm) Num. of Primitives are represented with three values between 0 and 1, where 0
lJC,x/y [2.2 : 0.2 : 4.0] and 1 represent the minimum and maximum of each parameter
JC lJC,cx/cy lJC,x/y − 1.6 mm 100 (JCs) &
& compJC wJC,x/y 0.4 100 (compJCs) in Table I, respectively. This is an important adjustment to
wJC,cx/cy 0.45 emphasize on the impact of each physical parameter in the
lP,x/y [2.0 : 0.2 : 5.0] 256 (RPs) & scatterer design. The order of the normalized dimensions are
RP
& comRP 256 (compRPs) based on their appearance in Table I for each scatterer. It
rSSR [1.4 : 0.2 : 2.6] is worth noting that the dependent and fixed dimensions of
SSR gSSR [0.3 : 0.2 : 1.5] 202
wSSR [0.2 : 0.2 : 1.2]
the scatterers, e.g. lJC,cx/cy and wJC,cx/cy for the Jerusalem
rCR [1.4 : 0.2 : 2.6]
cross, are not included in this representation. For all the
CR 49 shapes except the single split ring (SSR) shown in Figure
wCR [0.1 : 0.2 : 1.3]
lDB [2.2 : 0.2 : 4.8] 4 (c), the third value is fixed to 0 since only two physical
DB lDB,c [1.0 : 0.2 : 4.8] 280 parameters are to be tuned. This is done to keep the length
wDB 0.2
of the representation for all the unit cell combinations fixed.
Therefore, the top and bottom layer scatteres are described
by 8 variables each and the middle layer is described by 10
electric-field integral equation (EFIE) [40]. To handle general
variables. Adding the substrate thickness representation, this
cases including complementary scatterers, an extension on
results in a total of 30 variables to describe each sample
the existed MBPE method has been developed for magnetic-
3-layer unit cell in the training data. The mapping from
field integral equation (MFIE)-based MoM. Therefore, we are
the physical unit cell space to the 30-dimensional solution
capable to efficiently generate a large data set required for the
space is shown through an example in Figure 5. This type of
ML analysis including all types of periodic structures.
representation greatly helps the surrogate models to understand
3-layer EMMS training samples are generated by randomly
both the categorical and continuous aspects of the input data,
selecting three different GSMs (to implement bianisotropy)
i.e. physical structures of the EMMSs. We elaborate on this
and cascading them with different dielectric thicknesses in-
further after introducing the surrogate models in Section III-C.
cluding 0.254, 0.508, 0.787, and 1.575 mm. About 70000 sam-
ples out of about 2 billion possible combinations are generated
in this way. This electromagnetic-aware approach to generate
the training data set is more efficient than the traditional TABLE II. The code and categorical variables assigned to
approach of simulating the 3-layer unit cells. The training data each scatter shape and each substrate thickness.
generation was parallelized using the multiprocessing package Type Code Categorical Variables
from Python. The training set is denoted by Xt and the training Scatterer Type
subsets for each substrate thickness are denoted by Xt,h here. If on top/bottom layer: (0, 0, 0, 0, 1)
CR 0
If on middle layer: (0, 0, 0, 0, 0, 0, 1)
B. The EMMS Unit Cell Representation JC 1
If on top/bottom layer: (0, 0, 0, 1, 0)
To explore the design space, new 3-layer combinations of If on middle layer: (0, 0, 0, 0, 0, 1, 0)
the primitives shown in Figure 4 need to be evaluated beyond If on top/bottom layer: (0, 0, 1, 0, 0)
RP 2
those generated in the 70000 samples described in Section If on middle layer: (0, 0, 0, 0, 1, 0, 0)
III-A. For that, we exploit the samples in the training data If on top/bottom layer: (0, 1, 0, 0, 0)
SSR 3
set to train deep learning neural networks to provide fast If on middle layer: (0, 0, 0, 1, 0, 0, 0)
predictions of the scattering parameters of the candidate under If on top/bottom layer: (1, 0, 0, 0, 0)
DB 4
test. If on middle layer: (0, 0, 1, 0, 0, 0, 0)
The solution space contains different scatterer primitives, If on top/bottom layer: N/A
compJC 5
standard substrate thicknesses, and continuous scatterer di- If on middle layer: (0, 1, 0, 0, 0, 0, 0)
mensions. Therefore, we describe each EMMS unit cell with If on top/bottom layer: N/A
compRP 6
a combination of categorical and continuous variables. The If on middle layer: (1, 0, 0, 0, 0, 0, 0)
categorical parts are used to describe which scatterer primitive Substrate thickness
for each layer and what substrate thickness are used in the 0.254 mm 0 (0, 0, 0, 1)
unit cell while the continuous variables describe the relevant 0.508 mm 1 (0, 0, 1, 0)
dimensions of the scatterers for each layer. 0.787 mm 2 (0, 1, 0, 0)
Table II shows the code and the categorical variables given 1.575 mm 3 (1, 0, 0, 0)
to each scatterer shape and the substrate thickness. It is worth
noting that five types of the scatterer primitives shown in
Figure 4 (a)-(e) can be used for the top and the bottom C. Surrogate Models
layers, while all the shapes in Figure 4 (a)-(g) can construct We use the 30-variable representation, xinput , for the phys-
the middle layer of a unit cell. Therefore, five digit and ical structure of each unit cell along with the normalized
seven digit variables are used to describe the top/bottom- frequency point as the inputs to deep neural networks (DNNs),
layer and the middle-layer scatterers, respectively. Only the acting as surrogate models, to predict the magnitudes and
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 7

perparameters including the number of the hidden layers and

𝑤𝐶𝑅
= 0.3 mm the number of the neurons in each layer. More importantly,
CR
one can consider the inputs to the DNN the key factor in its
𝑟𝐶𝑅
= 2.0 mm
success.
5.3 mm Figure 6 (a) shows the mag-DNN. One might consider
CompRP using only xinput to predict a vector of scattering parameters,
𝑙𝑃,𝑥 = where elements of the vector represent the predicted scattering
4.4 mm
ℎ = 0.787 mm parameter at different frequency points in the band. However,
our tests reveal that it is simpler and more successful to train an
average-size DNN to predict the scattering parameters at each
frequency based on the frequency point input to the DNN. This
𝑤𝑆𝑆𝑅 =
1.2 mm SSR choice greatly helps the DNN to make accurate predictions
especially at the resonant frequency(ies) of the unit cell.
𝑔𝑆𝑆𝑅 𝑟𝑆𝑆𝑅
= 0.9 mm = 2.0 mm We implemented the DNNs using TensorFlow-backend
Keras libraries [41] in Python. The trained mag-DNN has six
hidden layers with 100, 500, 1000, 1000, 500, 100 neurons
SSR compRP CR ℎ with ReLU activation functions, from the most shallow to
the deepest layers, respectively. The output layer has four
𝑥𝑖𝑛𝑝𝑢𝑡 = (0,1,0,0,0, 0.5,0.5,1, 1,0,0,0,0,0,0,0.8,0,0, 0,0,0,0,1,0.5,0.2, 0, 0,1,0,0)
neurons for the magnitude of each TE-scattering parameter
𝑟𝑆𝑆𝑅,𝑛 𝑤𝑆𝑆𝑅,𝑛 𝑙𝑃,𝑥,𝑛 𝑟𝐶𝑅,𝑛 with sigmoid activation function due to the bounded values
𝑔𝑆𝑆𝑅,𝑛 𝑙𝑃,𝑦,𝑛 𝑤𝐶𝑅,𝑛
between 0 and 1. The DNN is trained with a batch size of
Figure 5. An example of how the physical structure of each unit cell 2048 using the ADAM optimizer [42] and backpropagation
is translated to a 30-variable input for the surrogate models.
method [43] on 85% of the training set and tested on the
remaining 15%. The loss function used to train the mag-DNN
is based on the mean squared error (MSE) between the actual
phases of the dispersive scattering parameters in an acceler- and predicted values as described by
ated. To find the optimized EMMS design, a loss function that
quantifies the difference between the desired surface properties Lmag−DN N (17a)
and the surface properties offered by the example under test X
needs to be evaluated. Therefore, for each unit cell that is = {MSE(|S11 |)+MSE(|S12 |)+MSE(|S21 |)+MSE(|S22 |)},
being examined, the surface properties need to be obtained
in the optimization loop. Traditionally, full-wave solvers are where
used to perform this part. However, due to the large amount TE TE
MSE(|Sij |) = |||Sij,actual |−|Sij,pred |||22 ; i, j ∈ 1, 2. (17b)
of time and resources they need, this step is the bottleneck
of the optimization through a large and high-dimensional
solution space such as the one of a bianisotropic EMMS.
To accelerate this part, we replace the full-wave solver with
surrogate models. Unlike the magnitudes of the scattering parameters, the
Provided with enough training examples, DNNs have scattering parameter phases are not limited to values between
proved to be efficient tools to learn the underlying patterns in 0 and 1. In fact, different unit cells have a different phase
large data sets and use this patterns to provide fast predictions range in the frequency band considered. Moreover, there are
for the features of the new inputs that they have never 180◦ - and 360◦ -phase jump(s) based on the order of unit
encountered before. These deep neural networks can then be cell’s responses. These features make the accurate predictions
integrated into optimization cycles as surrogate models to for the scattering parameter phases more challenging than the
accelerate the search for the optimized features. ones for the magnitudes because DNNs are better at predicting
Here, we choose to predict the scattering parameters of the bounded values within a specific range. Therefore, we consider
EMMS unit cell instead of Zse , Ysm , and Kem . This is because the following points to implement a DNN to predict each TE-
the magnitudes and phases of the scattering parameters are scattering parameter’s phase at different frequencies within the
intrinsically bounded and less prone to drastic changes at the band.
scatterer’s resonant frequencies. Without the loss of generality, 1) Adding the magnitude of the scattering parameter to
we focus on developing the surrogate models for TE-polarized the input of the phase-DNN: It is known that when
scattering parameters. The DNNs used to predict the magni- the magnitude of the scattering parameter is zero, its
tudes and phases of the scattering parameters are designated phase experiences a jump. Therefore, to predict the
mag-DNN and phase-DNN, respectively. The discussions here frequency of these phase jumps, we use the magnitude
can be applied to predict TM-scattering parameters in the information as an auxiliary input to the DNN besides the
same models of the TE-scattering parameters or in separate normalized frequency point and the 30-variable physical
models. However, building the right DNN that can accurately representation of the EMMS. This greatly improves the
predict the required features necessitates setting different hy- accuracy of the predictions made by the phase-DNN.
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 8

DNN is trained in a similar way to the mag-DNN with the

𝑓𝑟𝑒𝑞 following loss function to be minimized:

𝑥1
X
TE
𝑎1 Lphase−DN N = {MSE(]Sij,min )
𝑓1 𝑇𝐸
|𝑆11 | (19)
TE TE
+ MSE(]Sij,range ) + MSE(]Sij,n )}.
𝑥2 𝑎2 𝑓2
𝑇𝐸 TE
|𝑆12 | Firstly, one phase-DNN is trained to predict phases of S12 of
different unit cells at different frequency points. Then, transfer
𝑇𝐸 learning method [44] is employed to train similar neural
𝑎𝑀−1 𝑓𝐿−1 |𝑆21 | TE
𝑥26 networks to S12 ’s phase-DNN for the predictions of the rest
TE
of the scattering parameter phases, e.g. phase of S11 . In this
𝑎𝑀 𝑓𝐿 𝑇𝐸
|𝑆22 | method, instead of beginning the training of the new DNN
𝑥30
from random weights, the DNN is initialized with the weights
(a) of the pretrained S12 TE
’s phase-DNN. This approach requires
𝑓𝑟𝑒𝑞 much fewer number of epochs to minimize the prediction error,
which considerably reduces the training time.
𝑥1 Figures 7 (a) and (b) show the MSE distribution over the
𝑎1 𝑓1
training data set. It can be seen that most of the samples are
𝑇𝐸 predicted with an error of around 5.0×10−4 and 6.5×10−3 rad
𝑥2 𝑎2 𝑓2 ∡𝑆𝑖𝑗, 𝑚𝑖𝑛
by the mag-DNN and phase-DNN, respectively. The average
𝑇𝐸 error over the training set for mag-DNN and phase-DNN are
∡𝑆𝑖𝑗, 𝑟𝑎𝑛𝑔𝑒
8.8 × 10−4 and 2.4 × 10−3 rad, respectively.
𝑎𝑀−1 𝑓𝐿−1
𝑥26 𝑇𝐸
∡𝑆𝑖𝑗, 𝑛
𝑎𝑀 𝑓𝐿
𝑥30 𝑖, 𝑗 ∈ {1,2}

𝑇𝐸
|𝑆𝑖𝑗 |
(b)
Figure 6. Deep neural networks of the surrogate models to predict
(a) the magnitude of the TE-scattering parameters (mag-DNN), and
TE
(b) the phase of the S12 (phase-DNN).

0.5 1.5 2.5 3.5 4.5

TE TE TE
2) Predicting ]Sij,min , ]Sij,range , and ]Sij,n instead of MSE (× 10−3 )
TE
]Sij , where i, j ∈ {1, 2} and (a)
TE TE
]Sij,min = min(unwrap{]Sij }), (18a)
TE TE TE
]Sij,range = max(unwrap{]Sij } − ]S12,min ),
(18b)
TE TE
TE
unwrap{]Sij } − ]Sij,min
]Sij,n = TE
. (18c)
]Sij,range
TE TE TE TE
]Sij = ]Sij,n × ]Sij,range + ]Sij,min (18d)
TE
Instead of predicting ]Sij that can range from
TE TE
]Sij,min and ]Sij,range at each frequency, phase-DNN 0.65 1.8 3.0 4.2 5.4
TE
predicts ]Sij,min TE
and ]Sij,range that are constant for MSE (× 10−3 )
TE
a certain unit cell over the frequency points and ]Sij,n (b)
that is bounded between 0 and 1 for all unit cells. Figure 7. MSE distribution over all the training samples for (a) the
mag-DNN in Figure 6 (a), and (b) the phase-DNN in Figure 6 (b)
Both of these adjustments greatly help in understanding the TE
trained to predict ]S12 .
relation between the inputs and outputs by phase-DNN.
phase-DNN has six hidden layers with 100, 500, 2000,
2000, 500, 100 neurons with ReLU activation functions, from D. Optimization Over the Solution Space
the most shallow to the deepest layers, respectively. The Practical implementation of the EMMS dictates that a
TE TE
output layer has three neurons for ]Sij,min , ]Sij,range , and certain substrate thickness must be chosen for all the unit
TE
]Sij,n . One of these output neurons that is intended to provide cells that realize the desired surface parameters. To find which
TE
prediction of ]Sij,n has a sigmoid activation function and the substrate thickness must be used, the following fast approach
other two do not employ any activation functions. The phase- is implemented. For each EMMS unit cell, the minimum error
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 9

TE
between the desired scattering parameters and the training mostly impacted by the magnitude and phase of the S12 , we
samples of each substrate thickness, Xt,h , are calculated as modify the PSO loss function at 9.4 GHz to
L(h) (20)
i=N
X LP SO,i = ||S11,desired,i
= min(||S11,desired,i − S11,Xt,h ||22 + ||S12,desired,i − 2
S12,Xt,h ||2 ),
− S11,pred ||22 + |||S12,desired,i |−|S12,pred |||22 .
i=1

where N is the total number of unit cells in the EMMS. L(h) (23)
is calculated for each of the h ∈ {0.254, 0.508, 0.787, 1.525}
mm and the thickness that corresponds with the minimum To analyze the optimized EMMS for each case, we use
value of L(h) is chosen as the optimized sub-solution space scripts that read the file containing the primitives and dimen-
to be explored. sions of the EMMS scatterers and generate their model in
After finding which sub-domain provides the maximum Ansys HFSS. Once the 1D array of unit cells is generated,
number of matches to the desired scattering properties, the perfect boundary conditions (PEC) are assigned at x = −2.65
optimization through this space is performed as follows. To mm and x = 2.65 mm to replicate a periodic array in that
obtain the physical structure of the i-th unit cell of the direction and achieve a 2D simulation. This setup is shown in
EMMS, particle swarm optimization (PSO) is used. The i-th Figure 8.
desired scattering parameters converted from the optimized i-
th Zse , Ysm , and Kem obtained from Section II are defined as
the targets in the loss function of the PSO. This loss function
is defined as LP SO,i = |[S]desired,i − [S]pred |, where [S]pred
𝑧
is efficiently computed using the surrogate models developed
in Section III-C. PEC walls
The swarm of the PSO includes P particles and it is run Optimized
for I iterations. The positions of each particle are denoted 𝜃 EMMS
by xn , where n ∈ [1, P ] and x denotes the EMMS unit cell 𝑦
representation under test that has 30 variables. xn is updated
at the (m + 1)th iteration based on [45]

xn (m + 1) = xn (m) + vn (m + 1), (21)

𝑥
where

xn (m + 1) = w × vn,i (m) + c1 × (pn − xn ) + c2 × (pg − xn ).

(22)
pn is the particle’s historically best position and pg is the Feed
swarm’s best position regardless of which particle had found it. (𝑥𝑓 , 𝑦𝑓 , 𝑧𝑓 )
c1 and c2 are the cognitive and social parameters, respectively.
They control the particle’s behavior given two choices: (1) to
follow its personal best or (2) follow the swarm’s global best
position. Overall, this determines if the swarm is explorative
or exploitative in nature. In addition, a parameter w controls
the inertia of the swarm’s movement. The performance of the Figure 8. Simulation setup of the optimized EMMS.
PSO [45] is controlled by the choices of P , I, c1 , c2 , and
w. Here, these parameters are empirically selected to obtain
the best results. The PSO is implemented using the PySwarms
library in Python [46].
A. Example 1: 45◦ -Refraction Optimization
IV. D ESIGN E XAMPLES
In this example, the 15λ0 EMMS is excited with two line
In this section, two examples are provided with the far- sources that are located at xf = 0 mm, yf = ±D/8, and
field constraints including main beam(s) levels, sidelobe level, zf = −D/4. The far-field constraints are defined to transform
and null position applied on the transmitting side of the the excitation field to a beam directed at θ = 45◦ with sidelobe
EMMS. The design of both surfaces is carried out at 9.4 levels 12 dB below relative to the main beam for angles θ ∈
GHz where the unit cell’s period, i.e. 5.3 mm, is λ0 /6, where [0, 38]◦ and 10 dB below the main beam for the remainder of
λ0 is the free space wavelength. The EMMS has a size of the angles. Moreover, two null regions in the front, θ = 0◦ ,
D = 15λ0 composed of 90 unit cells. The EMMS unit cells and back, θ ∈ {179, 180, 181}◦ , of the EMMS are defined.
are located at z = 0 mm and distributed in the interval of The details of this optimization problem are listed in Table III
y ∈ [−7.5λ0 , 7.5λ0 ]. Since the field on the transmitted side is and (24).
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 10

TABLE III. 45◦ -Refraction optimization parameters

Parameter Value
Incident √ 1 e−jkr1 + √2r1
e−jkr2
2r1p 2
Field (E inc ) r1 = p(y − D/8)2 + (D/4)2
r2 = (y + D/8)2 + (D/4)2
Angular Sampling Points (M ) 361
Spatial Sampling Points (N ) 90
Max Iterations 295
Initial ([Xse ]0 , [Bsm ]0 , [Kem ]0 ) = 0,
Conditions ρ = 10, (µ0Ze , µ0Zm ) = 0
Main Lobe Angle (M B) θ = 45◦
Main Lobe Level (M Blevel ) 16.8 + j16.8 Solution from macro step
Sidelobe {7.4272, 0◦ ≤ θ ≤ 38◦ }, Solution from micro step
Level (τ ) {5.9679, 51◦ ≤ θ ≤ 180◦ },
and Angles (SL) {7.4272, −180◦ ≤ θ ≤ −90◦ },
{7.4272, −90◦ ≤ θ ≤ 0◦ }
TE TE
Null Angles (N U ) θ = {0◦ , 179◦ , 180◦ , 181◦ } Figure 10. |S12 | and ]S12 obtained from the macroscopic and
microscopic optimizations for the first example.

Unit cell #27

10000fM B (M B)+
minimize (24a)
Ie ,Im ,[Xse ],[Bsm ],[Kem ], 10000fN U (N U ) + fD
slackDe ,slackDm

subject to Einc = [Ze ]Ie + [Xse ]Ie − (24b)

m
[Kem ]I (24c) Unit cell
1000(H inc m
= [Zm ]I + #34
(24d)
[Bsm ]I + [Kem ]Ie )
m

|[Ge ](SL)Ie + [Gm ](SL)Im

(24e)
+ Einc
f f (SL)|≤ τ
|[D]Ie |≤ 0.1 + slackDe (24f) Figure 11. A section of the optimized EMMS for example 1.
m
|[D]I |≤ 25 + slackDm (24g)

Figure 12. Far-field (FF) patterns obtained from the homogenized

Figure 9. (a) Desired Zse , Ysm , Kem for transforming normal plane model (dashed red line) and the actual model (blue solid line)
wave from one side of the EMMS to 45◦ -directed pencil beam. simulated in Ansys HFSS for the specified main lobe, SLL, and null
The desired Zse , Ysm , and Kem for this transformation are constraints (magneta) at 9.4 GHz.
shown in Figure 9. These macroscopic properties are converted
to S-parameters to obtain S11,desired and S12,desired . The PSO
loss function based on these parameters is defined in (23). To middle layer but they are considered in the design space of
minimize this loss function, the sub-domain with h∗ = 1.525 the next example. Figure 10 shows the optimized S12 from the
mm is found as the sub-optimized domain to be explored as macroscopic step and S12 provided by the optimized physical
described in Section III-D. For this example as an initial test, unit cells in this domain. It can be seen that most of the desired
we have simplified the problem and limited the solution space properties are met with a small error.
to the unit cells with non-complementary primitives on their A subsection of the optimized EMMS is shown in Figure
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 11

11, where it is shown that the adjacent unit cells have com-
pletely different primitives. The optimized physical EMMS
is simulated and the far-field radiation pattern is shown in
Figure 12, where it is compared with the far-field obtained
from the optimized Zse , Ysm , and Kem based on the approach
explained in Section II. The comparison between the opti-
mized far-field and the imposed constraints shows that the
macroscopic step is successful in outputting the optimized
surface parameters. While the constraints are mostly met by
the synthesized EMMS, the null region in the simulated far-
field is not realized and the sidelobe levels are slightly higher
than the desired level. This is due the fact that the mutual
coupling between the unit cells has significant impact on both
of these properties. This physical mutual coupling between the
actual unit cell conductors is not fully captured by the inter-
cell mutual coupling predicted by the homogenized model
discussed in Section II. Despite this difference, the proposed Figure 13. Desired Zse , Ysm , Kem for transforming cylindrical wave
approach is successful in matching physical unit cells to a from the line source located at F = D/4 from one side of the EMMS
to the multi-beam far-field on the other side.
wide variety of surface parameters and generating the desired
radiation pattern. Figure 13 shows the Zse , Ysm , and Kem that are the
B. Example 2: Multi-Beam Optimization optimized solutions to the problem described in (25). Figure
TE TE
In the second example, we design an EMMS that transforms 14 shows the desired |S12 | and ](S12 ) that are converted
the fields from a line source located at xf = yf = 0 mm and from the desired Zse , Ysm , and Kem . Based on these desired
zf = −D/4 into an asymmetric multi-beam far-field radiation scattering parameters, the domain with the substrate thickness
pattern. The details of the optimization problem and the far- h∗ = 1.525 mm was chosen as the sub-solution space to
field constraints are listed in Table IV and (25). be explored. Unlike the previous example, in this example
to meet the desired surface properties, the domains with
both complementary and non-complementary primitives on the
TABLE IV. Multi-Beam optimization parameters middle layer are considered. The |S12TE
| and ](S12TE
) provided
Parameter Value by the optimized physical unit cells are shown in Figure 14,
Incident √1 e−jkr
pr where they are in excellent agreement with the optimized
Field (E inc ) r= y 2 + (D/4)2
Angular Sampling Points (M ) 361
parameters from the macroscopic step. Similar to the previous
Spatial Sampling Points (N ) 90 example, the optimized EMMS here is also composed of
Max Iterations 295 different primitives of scatterers. A subsection of the array
Initial ([Xse ]0 , [Bsm ]0 , [Kem ]0 ) = 0, is shown in Figure 15, where it can be seen that the unit cells
Conditions ρ = 10, (µ0Ze , µ0Zm ) = 0
Main Lobe Angles (M B) θ = {30◦ , −20◦ }
74 − 75 and 79 − 80 have complementary rectangular patch
Main Lobe {10.65 + j10.65, θ = 30◦ }, (compRP) scatterers on their middle layers.
Level (M Blevel ) {12.56 + j12.56, θ = −20◦ } Figure 16 shows the far-field pattern calculated from the
Sidelobe {3.85, 0◦ ≤ θ ≤ 11◦ },
Level (τ ) {3.95, 49◦ ≤ θ ≤ 90◦ },
desired Zse , Ysm , and Kem shown in Figure 13, the applied
and Angles (SL) {2.55, 90◦ ≤ θ ≤ 180◦ }, constraints, and the simulated far-field pattern from the phys-
{2.55, −180◦ ≤ θ ≤ −90◦ }, ical optimized EMMS. Despite the fact that the optimized
{3.85, −90◦ ≤ θ ≤ −36◦ },
{3.95, −4◦ ≤ θ ≤ 0◦ }
Null Angles (N U ) θ = {0◦ , 180◦ }

10000fM B (M B)+
minimize (25a)
Ie ,Im ,[Xse ],[Bsm ],[Kem ], 100fN U (N U ) + fD
slackDe ,slackDm

Einc = [Ze ]Ie + [Xse ]Ie

subject to (25b)
− [Kem ]Im
1000(Hinc = [Zm ]Im
(25c)
+ [Bsm ]Im + [Kem ]Ie )
Solution from macro step
|[Ge ](SL)Ie + [Gm ](SL)Im Solution from micro step
(25d)
+ Einc
f f (SL)|≤ τ
e
|[D]I |≤ 0.1 + slackDe (25e) TE TE
Figure 14. |S12 | and ]S12 obtained from the macroscopic and
m
|[D]I |≤ 25 + slackDm , (25f) microscopic optimizations for the second example.
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 12

where a combination of machine learning and optimization

Unit cell #73 techniques is employed to first obtain the optimized surface
parameters and then the physical EMMS. This problem has
many degrees of freedom including scatterer primitives and
dimensions and can have more than one solution. Therefore,
instead of using ad hoc and heuristic methods, we implement
the alternating direction method of multipliers (ADMM)-based
convex optimization to deterministically obtain the EMMS sur-
Unit cell #81
face parameters derived from the method of moments (MoM)
to satisfy constraints on the main lobe level(s), sidelobe levels,
and null positions on the radiated far-field. Then, the particle
swarm optimization integrated with deep neural networks used
as surrogate models is employed to explore and exploit the
Figure 15. A section of the optimized EMMS for example 2. physical solution space for the optimized EMMS.
Compared to the image-based representations of unit cells
[29], the proposed method to represent 3-layer bianisotropic
unit cells is more successful in capturing both the categorical
and continuous nature of the scatterers. Hence, it can interpo-
late the scattering properties of individual layers and interlayer
coupling in a new unit cell based on the samples in the
training data with reduced error, which removes the need for
further full-wave simulations. The surrogate models expedite
the exploration through the solution space so much that the
designer can choose large swarms or number of iterations for
the PSO and find the global optimum.
Two passive and lossless nonuniform bianisotropic meta-
surfaces are designed to transform the excited fields from
line source(s) to complex far-fields that successfully satisfy
the constraints on main lobe(s), sidelobes, and null regions.
The constituent unit cells of the optimized EMMSs can have
completely different primitives from their neighbors, which is
necessary to match the desired optimized surface parameters
Figure 16. Far-field (FF) patterns obtained from the homogenized
model (dashed red line) and the actual model (blue solid line) but it also results in discrepancies of the results due to the
simulated in Ansys HFSS for the specified main lobe, SLL, and null unpredictable levels of additional mutual coupling introduced
constraints (magneta) at 9.4 GHz. beyond that captured in the homogenized model. In the exam-
ples we have shown, the mutual coupling has mostly impacted
the null realization and sidelobe levels of the far-field radiated
constituent unit cells have different primitives and structures, it from the physical EMMS with acceptable error. Nevertheless,
can be seen that the optimized EMMS produces a far-field that the presented results demonstrate the effectiveness of the
satisfies the constraints everywhere except the maximum side- proposed approach in the inverse design of different EMMSs
lobe level for θ ∈ [−43, −35]◦ . It is worth mentioning that this for various types of application.
is due to the fact that the behavior of the unit cells in the non- Despite the demonstrated success of the method, several
periodic optimized EMMS deviates from their properties with practical considerations must be mentioned. First and fore-
the periodic boundary condition, due to the additional mutual most, the macroscopic optimizer described in Section II must
coupling between cells that are not modeled by the inter-cell be supplied with a roughly feasible problem to optimize for
coupling in the homogenized model. Nonetheless, through this if the designer is to meet their objectives. There are physical
example, we have shown that the proposed approach can be limitations to the directivity, sidelobe level, number of beams,
successfully employed for the inverse design of a nonuniform etc. that an EMMS can realize. Although the macroscopic
bianistropic metasurface that transforms the excitation from a optimizer can be used as a rough guide for the feasibility of a
line source to a complex far-field pattern comprising two main problem (i.e. if it converges or not) it is not rigorous. The
beams with different radiation intensities. macroscopic optimizer might also return optimized surface
parameters that may not be achievable from a 3-layer unit cell
V. C ONCLUSION given the constraints on substrate thicknesses and permittivi-
In this paper, the inverse design problem of a multilayer ties. Designers might also need to include scatterers with more
nonuniform metasurface based on the high-level far-field con- complex primitives and a wider range of scattering properties
straints is solved in a systematic and end-to-end approach in the solution space to match the desired surface properties
for the first time. We divide this problem into two inverse with a better accuracy. Furthermore, if the surface parameters
problems: a macroscopic problem and a microscopic problem, change abruptly and frequently across the surface, it means
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 13

that the physical structure of the EMMS also needs to follow [14] G. Xu, S. V. Hum, and G. Eleftheriades, “Augmented Huygens’ meta-
the properties, thereby causing more uncertainty due to the surfaces employing baffles for precise control of wave transformations,”
IEEE Trans. Antennas Propag., vol. 67, no. 11, pp. 6935-6946, Jun. 2019.
mutual coupling. To circumvent this issue, limiting constraints [15] R. Gomez-Bombarelli et al., “Automatic chemical design using a
on the spatial derivatives of the surface parameters can be datadriven continuous representation of molecules,” American Chemical
added to the macroscopic optimization problem, which will Society Central Science, vol. 4, no. 2, pp. 268-276, Jan. 2018.
[16] D. R. Prado, J. A. López-Fernández, G. Barquero, M. Arrebola , and
be considered in the future. F. Las-Heras, “Fast and accurate modeling of dual-polarized reflectarray
There remain some areas for future research to augment unit cells using support vector machines,” IEEE Trans. Antennas Propag.,
the utility of this method. Firstly, a method to automatically vol. 66, no. 3, pp. 1258-1270, Mar. 2018.
determine the optimization weights of the macroscopic opti- [17] T. Qiu et al, “Deep learning: a rapid and efficient route to automatic
metasurface design,”Adv. Sci., no. 6, pp. 1900128 (1-12), 2019.
mizer would reduce the amount of tuning required. Secondly, [18] V. Richard, R. Loison, R. Gillard, H. Legay, and M. Romier, “Loss
a minimum level field constraint would also be very useful. analysis of a reflectarray cell using ANNs with accurate magnitude
This would enable the specification of minimum directivity prediction,” in Proc. 11th Eur. Conf. Antennas Propag. (EuCAP), Paris,
France, Mar. 2017, pp. 2402–2405.
levels for beams or more sophisticated patterns (such as [19] D. Kampouridou and A. Feresidis, “Machine learning-driven design
isoflux patterns). Thirdly, moving this design scheme into optimization for a multi-layer metasurface antenna,” 14th Eur. Conf.
three dimensions would provide more degrees of freedom Antennas Propag. (EuCAP), Copenhagen, Denmark, Mar. 2020.
[20] D. Caputo, A. Pirisi, M. Mussetta, A. Freni, P. Pirinoli, and R. Zich,
to optimize for objectives at different elevations. Finally, a “Neural network characterization of microstrip patches for reflectarray
fully generative machine-learning model [29] combined with optimization,” in Proc. 3rd Eur. Conf. Antennas Propag. (EuCAP), Berlin,
generalized reliable surrogate models can be integrated into Germany, Mar. 2009, pp. 2520-2522.
the microscopic optimization step to efficiently explore the [21] P. Robustillo, J. Zapata, J. A. Encinar, and J. Rubio, “ANN charac-
terization of multi-layer reflectarray elements for contoured-beam space
physical solution space of 3-layer bianisotropic unit cells to antennas in the Ku-band,” IEEE Trans. Antennas Propag., vol. 60, no. 7,
a greater degree and achieve desired surface parameters with pp. 3205-3214, Jul. 2012.
reduced error. [22] M. Salucci et al., “Efficient prediction of the EM response of reflectarray
antenna elements by an advanced statistical learning method,” IEEE
Trans. Antennas Propag., vol. 66, no. 8, pp. 3995-4007, Aug. 2018.
R EFERENCES [23] Z. Liu, D. Zhu, S. P. Rodrigues, K-T. Lee, and W. Cai, “A generative
model for inverse design of metasurfaces,” Nano Lett., vol. 18, no. 10,
[1] O. Quevedo-Teruel et al, “Roadmap on metasurfaces,” Journal of Optics, pp. 6570-6576, Sep. 2018.
vol. 21, no. 7, pp. 073002 (44pp), Aug. 2019. [24] X. Shi1, T. Qiu, J. Wang, X. Zhao, and S. Qu, “Metasurface inverse de-
[2] K. Achouri, M. A. Salem MA, C. Caloz, “General metasurface synthesis sign using machine learning approaches,” Journal of Physics D: Applied
based on susceptibility tensors,” IEEE Trans Antennas Propag, vol. 63, Physics, vol. 53, no.27, pp. 275105 (7pp), 2020.
no. 7, pp. 2977-2991, Jul. 2015. [25] J. Jiang, D. Sell and J, A. Fan, “High efficiency metasurface design
[3] A. Epstein and G. V. Eleftheriades, “Arbitrary power-conserving field based on deep generative models,” Advanced Photonics Congress (IPR,
transformations with passive lossless omega-type bianisotropic metasur- Networks, NOMA, PVLED, SPPCom), OSA 2019.
faces,” IEEE Trans. Antennas Propag., vol. 64, no. 9, pp. 3880–3895, [26] S. An et al, “Multifunctional metasurface design with a generative
2016. adversarial network,” arXiv:1908.04851, Aug. 2019.
[4] M. Chen and G. Eleftheriades, “Omega-bianisotropic wire-loop huygens’ [27] J. Jiang and J. A. Fan, “Simulator-based training of generative models
metasurface for reflectionless wide-angle refraction,” IEEE Trans. An- for the inverse design of metasurfaces,” Nanophotonics, vol. 75, no. 5,
tennas Propag., vol. 68, no. 3, pp. 1477-1490, Mar. 2020. pp. 1059–1069, Nov. 2019.
[5] M. Selvanayagam and G. V. Eleftheriades, “Discontinuous electromag- [28] W. Ma, F. Cheng, Y. Xu, Q. Wen, and Y. Liu, “Probabilistic represen-
netic fields using orthogonal electric and magnetic currents for wavefront tation and inverse design of metamaterials based on a deep generative
manipulation,” Opt. Express, vol. 21, no. 12, pp. 14409–14429, jun 2013. model with semi-supervised learning strategy,” Advanced Materials, vol.
[6] G. Oliveri et al., “Synthesis of multi-layer WAIM coatings for planar 31, no. 35, pp. 1901111 (9 pp.), 2019.
phased arrays within the system-by-design framework,” IEEE Trans [29] P. Naseri and S. V. Hum, “A generative machine learning-based approach
Antennas Propag, vol. 63, no. 6, pp. 2482-2496, Jun. 2015. for inverse design of multilayer metasurfaces,” IEEE Trans. Antennas
[7] P. Naseri et al., “Efficient evaluation of gradient transmit-arrays through Propag., Feb. 2021 (Early Access).
an equivalent dispersive dielectric description,” IEEE Trans Antennas
[30] G. Gosal, “The use of inverse neural networks in the fast design of
Propag, vol. 67, no. 9, pp. 5997-6007, May 2019.
printed lens antennas,” [Link]. Thesis, Electrical & Computer Engi-
[8] J. Budhu and A. Grbic, “Perfectly reflecting metasurface reflectarrays:
neering Department, University of Ottawa, Ottawa, ON, Canada, 2015,
mutual coupling modelling between unique elements through homoge-
[Link]
nization,” IEEE Trans. Antennas Propag., vol. 69, no. 1, pp. 122–134,
[31] G. Gosal, E. Aljamali, D. McNamara, and M. Yagoub, “Transmitarray
2020.
antenna design using forward and inverse neural network modeling,”
[9] T. Brown, Y. Vahabzadeh, C. Caloz, and P. Mojabi, “Electromagnetic
IEEE Ant. Wireless Prop. Lett., vol. 15, no. 1, pp. 1483-1486, 2016.
inversion with local power conservation for metasurface design,” IEEE
Antennas Wirel. Propag. Lett., vol. 19, no. 8, pp. 1291 – 1295, 2020. [32] G. Oliveri et al., “System-by-design multi-scale synthesis of task-
[10] M. Salucci, A. Gelmini, G. Oliveri, N. Anselmi, and A. Massa, “Synthe- oriented reflectarrays,” IEEE Trans. Antennas Propag., vol. 68, no. 4,
sis of shaped beam reflectarrays with constrained geometry by exploiting pp. 2867-2882, Apr. 2020.
nonradiating surface currents,” IEEE Trans. Antennas Propag., vol. 66, [33] C. Yeung et al., “Multiplexed supercell metasurface design and opti-
no. 11, pp. 5805–5817, 2018. mization with tandem residual networks,” Nanophotonics, vol. 10, no. 3,
[11] H.-D. Lang, S. V. Hum, and C. D. Sarris, “Optimization of reactively pp. 1133-1143, Jan. 2021.
loaded reflectarrays via semidefinite relaxation,” in 2018 IEEE Int. Symp. [34] J. Jiang and J. Fan, “Global optimization of dielectric metasurfaces using
Antennas Propag. Usn. Natl. Radio Sci. Meet., Boston, MA, 2018, pp. a physics-driven neural network,” arXiv: 1906.04157v2, Jul. 2019.
1593–1594. [35] J. Budhu, E. Michielssen, and A. Grbic, “The design of dual band
[12] M. Bodehou, C. Craeye, E. Martini, and I. Huynen, “A quasi-direct stacked metasurfaces using integral equations,” arXiv:2103.03676, Feb.
method for the surface impedance design of modulated metasurface 2021.
antennas,” IEEE Trans. Antennas Propag., vol. 67, no. 1, pp. 24–36, Jan [36] S. Pearson and S. V. Hum, “Optimization of scalar and bianisotropic
2019. electromagnetic metasurface parameters satisfying far-field criteria,”
[13] D. R. Prado, “Advanced techniques for the analysis and synthesis of arXiv:2011.09016, Nov. 2020.
reflectarray antennas with applications in near and far fields”, Ph.D. [37] V. G. Ataloglou, M. Chen, M. Kim, and G. V. Eleftheriades, “Microwave
Thesis, Technologies of Information and Communications, University of Huygens’ Metasurfaces: Fundamentals and Applications,” IEEE J. Mi-
Oviedo, Oct. 2016. crowaves, vol. 1, no. 1, pp. 374–388, 2021.
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2021 14

[38] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed Op-

timization and Statistical Learning via the Alternating Direction Method
of Multipliers,” Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1–122,
2011.
[39] C. Wan and J. A. Encinar, “Efficient computation of generalized scatter-
ing matrix for analyzing multilayered periodic structures,” IEEE Trans.
Antennas Propag., vol. 43, no. 11, pp. 1233–1242, Nov. 1995.
[40] Z. Wang and S. V. Hum, “A broadband model-based parameter estima-
tion method for analyzing multi-layer periodic structures,” IEEE Trans.
Antennas Propag., Feb. 2021 (Early Access).
[41] F. Chollet et al., Keras, 2015, Available at:
[Link]
[42] D. P. Kingma and J. L. Ba, “Adam: a method for stochastic opti-
mization,” International Conference on Learning Representations (ICLR),
pages 1–13, San Diego, CA, USA, May 2015.
[43] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning internal
representations by backpropagating errors,” Nature, vol. 323, no. 6088,
pp. 533-536, 1986.
[44] K. Weiss, T.M. Khoshgoftaar, and D. Wang, “A survey of transfer
learning”, Journal of Big Data, vol. 3, no. 9, 2016.
[45] Kennedy and R.C. Eberhart, “Particle swarm optimization,” in Proc. of
the IEEE International Joint Conference on Neural Networks, pp. 1942-
1948, Nov. 1995.
[46] L. J. Miranda, “PySwarms: a research toolkit for Particle Swarm
Optimization in Python,” Journal of Open Source Software, 3(21), 433,
2018.