Optimisation of flat dielectric lenses using an interior point method

Optimisation of flat dielectric lenses using an

interior point method

Department of Mathematics, Linköping University Jonatan Ek

LiTH-MAT-EX–2021/02–SE

Credits: 30 hp Level: A

Supervisor: Yurii Malitskyi,

Department of Mathematics, Linköping University Examiner: Fredrik Berntsson,

Department of Mathematics, Linköping University Linköping: June 2021

Abstract

This thesis aims to study how flat dielectric lenses can be designed. The usage of flat lenses is steadily increasing as they are smaller and less bulky than tra-ditional convex lenses. Instead of a lens with a curved surface the permittivity in the lens is varied to achieve the same effect. Two different computational methods were investigated when approaching this problem: physical and geo-metrical optics. In physical optics the incoming radio waves are treated as waves in contrast to geometrical optics where it is considered as rays. Both methods are used as approximations of Maxwell’s equations.

The variation of permittivity in the lens was formulated as an optimisation problem where the lens’ focusing abilities were maximised. The optimisation was implemented with an interior point method. Both arbitrary permittivity distributions as well as predetermined distributions were examined in this work. All optimised lens models were then simulated in a full wave commercial simu-lation software to verify and compare the two.

The simulations showed that both approaches gave promising results as they focused the electromagnetic wave in a satisfying way. However the physical op-tics approach was more prominent as the focused radio waves had a much higher magnitude than the approach based on geometrical optics. The conclusion was therefore that physical optics is the preferred approach.

Keywords:

Lens, Permittivity, Optimisation, Interior Point Method, Wave Propaga-tion

Acknowledgements

This master thesis marks the end of my education in Applied Physics and Electrical Engineering, Master of Science degree in Applied Mathematics at Linköping University.

This work was carried out at Saab Dynamics in Linköping. Thanks to everyone who helped me in some way during this thesis work. Especially my supervisors Philip Bergander and Sören Poulsen deserve to be mentioned. Thanks for fruit-ful discussions and guidance throughout this work. Your help was invaluable. I also want to thank Fredrik Berntsson who took the responsibility of being the examiner for this project and Yurii Malitskyi who acted as supervisor. Lastly I want to direct a big thanks to my family and friends. You supported and motivated me along this journey and for that I am truly grateful.

Contents

1 Introduction 1 1.1 Background . . . 1 1.2 Purpose . . . 2 1.3 Goal . . . 2 1.4 Approach . . . 2 2 Theory 3 2.1 Maxwell’s equations . . . 3 2.2 Geometrical optics . . . 4 2.2.1 Snell’s law . . . 5 2.3 Physical optics . . . 6 2.4 Poynting vector . . . 9 2.5 Reciprocity . . . 9 2.6 Related works . . . 10 3 Method 11 3.1 Interior point algorithm . . . 11

3.2 Geometrical optics implementation . . . 17

3.2.1 Arbitrary distribution . . . 18

3.2.2 Radial distribution . . . 19

3.3 Physical optics implementation . . . 20

3.4 Simulation and verification . . . 21

4 Results and discussion 25 4.1 Implementation results . . . 26 4.1.1 Geometrical optics . . . 26 4.1.2 Physical optics . . . 28 4.2 Simulation results . . . 29 4.3 Comparison . . . 34 Ek, 2021. vii

Chapter 1

Introduction

Lenses have a wide area of use. Cellphone cameras, radar systems and eye-glasses are just a few examples. Lenses are traditionally curved which enables them to redirect, focus and scatter radio waves. However, the curvature makes the lenses bulky which is not desirable when many of the applications push for smaller and smaller devices. The curvature can also be hard to construct when the structure is complicated and the dimensions are small. Because of this, a new type of lens was introduced [13] which instead has a flat surface. For it to still be able to steer and focus radio waves the same way as before it now relies on permittivity changes within the lens. How the permittivity should vary within the lens depends on the desired effect of the lens.

This thesis will study flat dielectric lenses and focus on how the permittiv-ity should vary in the lens. Two different approaches will be tested: a physical and a geometrical optics approach. To find the best permittivity profile for the lens it will be formulated as an optimisation problem which is solved with an interior point method. In the optimisation of the lenses we try to optimise their ability to focus radio waves. Focus will in this work mean to convert a spherical wave to a plane wave or conversely.

1.1

Background

This project was proposed by Saab Dynamics in Linköping. Saab produces radomes which is a dielectric structure and is therefore knowledgeable in elec-tromagnetic calculations and modelling. They were interested in how flat di-electric lenses could be used in measuring equipment. The interest resulted in

this master thesis which studies how such a lens can be implemented in differ-ent ways and how well it performs. The focus of this thesis is not on dielectric lenses specific for radar usage, however the lenses will be designed for waves in the radio spectrum.

1.2

Purpose

The implementation of a flat dielectric lens is central for this project. The purpose is to see how well such a lens can perform. Since it is a dielectric lens some sort of modelling of the electromagnetic waves is needed. One part of this thesis is therefore to see how different ways of modelling affect the end result. To compare the result of the different ways of modelling, optimisation of the lenses will be done. The purpose of the optimisation is to see the full potential of the lenses. This project is theoretical which means the implementation will consist of modelling, visualisation and simulation.

1.3

Goal

The goal for the project is to model a flat dielectric lens that can focus radio waves in a satisfying way. To quantify this measure a comparison will be made with lenses from other works. The lens should be able to focus the incoming radio waves to a desired point. Optimisation of the lens will be done so that we get maximum efficiency of a lens with a certain size. Simulations will be used for analysis.

1.4

Approach

This thesis will be constructed as following. First the theoretical basis of this thesis is introduced which provides the reader with all necessary background information to understand the rest of the work. Then the chosen implementation of computational and optimisation methods will be shown before results are presented. Results will be discussed and analysed before the report is concluded.

Chapter 2

Theory

This chapter introduces the theory that is necessary to understand the com-putational and optimisation methods in the report. The main focus is on the propagation and modelling of electromagnetic waves. It is divided in multiple sections to distinguish between the different areas. Implementations and results later on in the report are based on the contents of this chapter.

2.1

Maxwell’s equations

Maxwell’s equations form the foundation of classical electromagnetic theory. They relate the electric field and the magnetic field to charge and current [4]. The four equations, in time harmonic form, are

∇ × H = J + jω0rE, ∇ × E = −jωµ0µrH, ∇ · B = 0 and ∇ · D = ρ, (2.1)

where H, B are magnetic fields and E, D electric fields. J is the electric current density and ρ is the electric charge density. We also have relations between the fields in the form of

D = r0E

and B = µrµ0H,

(2.2)

where µ0 and 0 are constants, and r and µr are relative permittivity and

rel-ative permeability respectively. Equation 2.2 only holds for isotropic materials which is described later on in this section. Relative permittivity or the dielectric constant is a value that affects how the electric field behaves in a medium [4]. Relative permeability is similar but for magnetic fields. In vacuum the relative permittivity is one while other materials have increasing values. The relative permittivity is the most important aspect of this work since this is the parameter which will be optimised in the material to affect how the electromagnetic wave propagates in the lens. Instead of having curved dielectric lenses, the values for the relative permittivity will vary throughout the lens, converting the incoming wave to a desired transmitted wave. Higher relative permittivity leads to the wave propagating slower.

Relative permeability and relative permittivity can take a wide range of val-ues. In some materials the properties can differ in different directions of the medium. Materials like these are called anisotropic and are more complex to work with. Materials that have uniformity between the different directions are called isotropic non-magnetic materials. In this work only isotropic materials will be considered. Hence, relative permeability is also not relevant to this work, henceforth it will always be assumed to be 1, the same as in vacuum. For com-pactness relative permittivity will be referred to as permittivity.

The connection between refractive index and permittivity is

n =√rµr. (2.3)

Since the relative permeability is assumed to be 1 in this work the relation is straightforward and easy to use.

The solutions to Maxwell’s equations gives a complete description of the re-lation between electromagnetic fields and currents and charges. They can also be used to describe how the fields propagate in different media. However, the complicated structure of nested differential equations makes this a hard prob-lem to solve analytically. Instead approximations are frequently used. Since we are interested in dielectric lenses for radio waves, approximations for high frequencies are studied.

2.2

Geometrical optics

One way of modelling the electromagnetic wave propagation is with optics in a geometrical sense [11]. Here we assume that the electromagnetic field propagates

2.2. Geometrical optics 5

as light rays. This assumption holds well as long as the dimensions of the physical objects involved are much greater than the wavelength [6]. With this assumption, crucial laws that describe the propagation of light rays can be postulated.

2.2.1

Snell’s law

Snell’s law, or the law of refraction, explains how the direction of light rays change when entering a new medium. The law is

n1sin θ1= n2sin θ2 (2.4)

where n1, n2are the two refractive indices, θ1is the angle of incidence and θ2 is

the angle of refraction [11]. When a ray moves from one medium into another it will change direction, refract, assuming the new medium has a different refrac-tive index, see Figure 2.1. This is commonly used in lenses to steer light rays in desired ways. The angle of incidence and angle of refraction are measured with respect to the normal of the surface. Depending on if n2 or n1 is larger, the ray

will either bend towards or from the normal (in Figure 2.1 n2 > n1). Not all

of the light ray will refract, some of it will reflect instead. For those rays the angle of reflection will always be the same as the angle of incidence according to Snell’s law for reflection [11].

θ2

n1

n2

θ1

n1sin θ1= n2sin θ2

Figure 2.1: Snell’s law

Since some of the light ray will be reflected instead of transmitted, the energy of the transmitted ray will decrease. We can denote this with R for reflectance and T for transmittance. The following equations are under the assumption that absorption is zero and that the permeability µ1 = µ2 = 1. They can be

the assumption of no absorption), according to the conservation of energy. The equations for R and T are given by Fresnel’s equations which takes into account the polarisation of light. This leads to two equations, one for each polarisation (s and p), for both R and T. Fresnel’s equations are [11]

Rs=  n1cos θi− n2cos θt n1cos θi+ n2cos θt 2 , Rp=  n2cos θi− n1cos θt n2cos θi+ n1cos θt 2 . (2.5) Here n1 and n2 are the refractive indices, θi is the angle of incidence and θtis

the angle of transmission. Using Snell’s law and trigonometric identities we can write them as Rs=  sin(θi− θt) sin(θi+ θt) 2 , Rp=  tan(θi− θt) tan(θi+ θt) 2 . (2.6)

These equations are not defined when the angle of incidence and the correspond-ing angle of transmittance are zero. In that special case we can rewrite Equation 2.5 into a single equation

Rs= Rp=

 n2− n1

n1+ n2

2

. (2.7)

These equations are useful when designing optical equipment. In some appli-cations there are desired values of how big the reflectance or transmittance is allowed to be which is easy to calculate with Fresnel’s equations.

2.3

Physical optics

Physical optics (PO) is one way of approaching the modelling of electromag-netism. It is similar to geometrical optics (GO) since it’s also a high-frequency approximation method. However it is a less rough approximation because it takes into account wave effects. It can be seen as a middle-ground between geometrical optics and full wave electromagnetism. Radio waves can still be vi-sualised as rays even though it now has wave properties that should hold. Since it is also a high-frequency approximation we assume wavelengths are smaller than the dimensions of the physical objects. Due to this assumption all waves can locally be seen as plane waves at the boundary of a media. The area of physical optics is large and not everything is of importance to this work. This section will therefore only present the relevant theory and results, for full details see for example [1] and [7]. In this report physical optics will mostly be seen as a framework which enables methods to be feasible.

2.3. Physical optics 7

We are interested in calculating how electric and magnetic fields propagate through a lens. In physical optics this can be seen as a scattering problem which is calculated through integral equations. These equations are characterised by having the unknown quantity in the integrand. The integral representation for the scattered electric field can [1], in a simplified way, look like this

ES = 1 C∇ × ∇ × Z Z S JS(r)e−ikrdS (2.8)

where JS is the electric current density over a surface S that is being radiated

by a source and C is a constant. However, JSis not known and is difficult to

cal-culate, even with numerical methods. Instead the physical optics approximation is used

JS = 2ˆn × Hi, (2.9)

where ˆn is the normal to the surface S and Hi _{is the incoming magnetic field}

generated by the source. This holds for the part of S that is illuminated, other-wise JS is zero. This approximation allows us to calculate the scattered electric

field and further how the field changes throughout the lens.

Another way of calculating the propagation of fields within physical optics is with a propagator. In [15], Rikte, Kristensson and Andersson demonstrate how to derive the expression for the propagator straight from Maxwell’s equations for planar structures. The propagator is an operator that takes an incoming field and a medium and returns the outgoing field. In the article they assume that the medium is bianisotropic. Since we are only dealing with isotropic materials, which can be seen as a special case, only parts of their work is relevant here. The concept of a dyadic [12] is also used in the propagator. However, in the implementation the dyadics are handled as matrices. A dyadic D is constituted by two vectors a = (a1, a2)T and b = (b1, b2)T as

D = ab = abT =a1b1 a1b2 a2b1 a2b2



. (2.10)

The propagator takes the incoming fields as input and returns the fields at the end of some media. This creates the relationship

 Exy(z) η0J · Hxy(z)  = P(z, z1) ·  Exy(z1) η0J · Hxy(z1)  (2.11)

where P is the propagator in a media that stretches from z1 to z. J is a

two-dimensional rotation dyadic (in this case rotation of π/2 in the xy-plane) and η0=pµ0/0 is the impedance in vacuum. The propagator is defined as

P (z, z1) = eik0(z−z1)M = I4cos(k0(z − z1)λ) +

i

λM sin(k0(z − z1)λ). (2.12) Here λ2 _{= }

r− k2t/k02 where k0 and kt are the vacuum wavenumber and the

tangential wavenumber. The connection between the two is kt= k0sin(φi)with

φi being the angle of incidence. The fundamental dyadic M can be found from

M = 0 −I2+ 1 rk02 ktkt −rI2−_k12 0 J · ktkt· J 0 ! (2.13) where kt is the tangential wave vector and can be written as kt= ktˆek. Note

that it is the dyadic product between the tangential wave vector and itself. The vector ˆek and its counterpart ˆe⊥ can be seen in Figure 2.2.

y ˆ e⊥ kt ˆ ek kt φi _x

Figure 2.2: Overview of vector system

In the figure φi is as previously the angle of incidence. We can formulate the

connection between the cartesian coordinates and the local coordinate system as ( ˆ e_k= cos φix + sin φˆ iy,ˆ ˆ e⊥= − sin φix + cos φˆ iy.ˆ (2.14) With this notation introduced and d = z − z1we can reformulate the final form

2.4. Poynting vector 9 of the propagator as P (z, z1) =  (ˆekˆek+ ˆe⊥ˆe⊥) cos(k0dλ) −i(ˆekˆekλr + ˆe⊥eˆ⊥ 1 λ) sin(k0dλ)

−i(ˆekeˆkλr + ˆe⊥ˆe⊥λ) sin(k0dλ) (ˆekeˆk+ ˆe⊥ˆe⊥) cos(k0dλ)

 . (2.15) In Equation 2.15 it is the dyadic product between the local vectors. The prop-agator is now defined for propagating the fields through one layer of media. The parameters needed are the angle of incidence, the relative permittivity in the media, the frequency and the thickness of the media. One important detail to note is that the propagator assumes the media has a homogeneous relative permittivity. This is not an issue since the propagator can be repeatedly called every time the relative permittivity changes. The propagator also assumes that the homogeneous surface is infinitely wide. Since the surfaces in our modelling is finite the propagator will be seen as an approximation.

2.4

Poynting vector

The Poynting vector represents the energy flux density and is a way of seeing what direction and magnitude the energy in a electromagnetic field has [1]. It is closely connected with the electric and the magnetic field as the definition is

S = 1

2Re(E × H

∗_{). (2.16)}

When visualised, the poynting vector creates a good overview of the intensity of electromagnetic fields, often in a heat map-format. Hence, it is useful to visualise outputs to see how intense the fields are in different points of space. The poynting vector will especially be useful when studying the lens’ focusing abilities.

2.5

Reciprocity

Reciprocity is a concept that explains a connection between an electromagnetic system and its inverse. Electromagnetic system translates in this case to a radiating source, some media and a receiver. Reciprocity states that the source and the receiver can swap places, creating an inverted system, without changing the observed wave [14]. Therefore one can design a media, in this case a lens, that converts waves from one type to another and then inverse the system to do the opposite. This way one can easily change between spherical waves and plane waves as input depending on what the purpose is, without having to redo the design of the lens.

2.6

Related works

The concept of replacing convex dielectric lenses with flat equivalents has been studied in many other works and is not unique to this project. In this section we will therefore briefly look how some of these works has approached the problem and what they concluded.

The authors Isakov et al. produced an article [9] which studied how a flat dielectric lens with a radial permittivity distribution could be used in combi-nation with an antenna. The goal was to direct the signal of the antenna into a narrow beam. With transformation optics, similar to geometrical and phys-ical optics, they derived a formula for the permittivity of the lens. They then 3D-printed the lens and tested it in practice. The 3D-printed lens significantly reduced the size of the system with no loss of ability to direct the signal. In [10] a physical optics approach, similar to the one previously described, was used. The authors optimised the flat lens such that the transmission was max-imised while still producing a plane wave as output. The resulting lens had five layers, each with a declining permittivity as the radius grew. Symmetry was naturally present in every layer but also between the layers, i.e.layer 1 and 5 was the same. The results from the study indicated that this approach is preferred over transformation optics and ray optics. This was concluded based on the transmission of the lens and also the simplicity of the method.

A simple model of a flat dielectric lens was presented by Zhang in [17]. Also here the lens was intended to be used in combination with an antenna to direct the radiation. The simplicity of Zhang’s model makes for a straightforward way to calculate the permittivity. The focus in this article was not to maximise transmission but to obtain good directivity. The results were good in the sense that the side lobes were significantly lower than the main lobe.

Lastly we will look at a conference paper [8] written by He and Eleftheriades. They present a lens that has a gradient in its permittivity profile. The focus of the paper was to produce a flat lens with low reflections. The solution they came up with was to coat the lens with multiple anti-reflection layers. These anti-reflection layers were designed such that the permittivity had a smooth in-crease or dein-crease rather than having a large discrete jump. With these layers the reflections were severely lowered and the transmittance was only slightly decreased for the lens.

Chapter 3

Method

Chapter 3 describes the methods used in this study. Firstly, a proper introduc-tion to the optimisaintroduc-tion algorithm is made. Some derivaintroduc-tions are made and the important steps are described. Secondly the implementation of the geometrical optics approach is described. Here the integration of optimisation and wave modelling is motivated. Lastly the physical optics implementation is explained. Similarly to the second section but instead of geometrical modelling a physical optics model is introduced.

3.1

Interior point algorithm

The optimisation algorithm that is used in this work is the interior point algo-rithm [5]. It is efficient for solving nonlinear optimisation problems. Interior point methods are an alternative to active set methods (which follows the edge of the feasible region) when it comes to handling inequality constraints. The basis of the algorithm is sequential quadratic programming which is strong when the constraints are non-linear. This section will go trough the fundamentals of the algorithm as well as the underlying theory. Full derivations and details can be seen in [2] and [3].

We start by formulating the optimisation problem as min x f (x) subject to ( h(x) = 0, g(x) ≤ 0, (3.1) Ek, 2021. 11

where f : Rn _{→ R, h : R}n _{→ R}` _{and g : R}n _{→ R}m _{are smooth functions.}

It is assumed that the first and second derivatives of the functions exist. We introduce slack variables s to have equality constraints

min x f (x) subject to      h(x) = 0, g(x) + s = 0, s ≥ 0. (3.2)

Instead of using an active set approach, the interior point method is chosen. The problem then turns into a barrier problem where a penalty parameter µ > 0 appears: min x f (x) − µ m X i=1 ln si subject to      h(x) = 0, g(x) + s = 0, s > 0. (3.3)

Trivially, as µ approaches 0, the solutions will converge towards the solution of the original problem. The algorithm aims to approximately solve problem (3.2) for fixed values of µ. To solve the problem it is needed to fulfil the optimality conditions. We begin by introducing some notation to be able to write the conditions in a simple way. Firstly, the Lagrangian of (3.3) is

L(x, s, λh, λg) = f (x) − µ m X i=1 ln si+ λThh(x) + λ T g(g(x) + s), (3.4)

where λh∈ R`, λg∈ Rmare Lagrange multipliers of the equality and inequality

constraints, respectively. An optimal solution (x, s) will satisfy ∇xL(x, s, λh, λg) = ∇f (x) + Ah(x)λh+ Ag(x)λg= 0

and

∇sL(x, s, λh, λg) = −µS−1e + λg= 0,

(3.5) where Ah(x) = (∇h1(x), . . . , ∇h`(x)), Ag(x) = (∇g1(x), . . . , ∇gm(x)) are

con-straint gradient matrices and where e = (1, . . . , 1)T_{, S = diag(s}

3.1. Interior point algorithm 13 simplicity we introduce z = (x, s)T, ϕ(z) = f (x) − µ m X i=1 ln si, c(z) = (h(x), g(x) + s)T, (3.6) which enables us to rewrite Problem 3.3 as

min

z ϕ(z)

subject to c(z) = 0. (3.7)

To proceed with solving the problem the sequential quadratic programming method is used. Sequential quadratic programming is an iterative method with step size d = (dx, ds)T. The quadratic problem is formulated

min ∇ϕ(z)T_{d +}1 2d T_{W d} subject to ( ˆA(z)Td + c(z) = 0, kdk2≤ ∆, (3.8)

where ˆA(z)T _{is the Jacobian of c(z) and is defined as}

ˆ A(z)T =Ah(x) T ₀ Ag(x)T I  , (3.9)

and W is the Hessian of the Lagrangian, see (3.4), and is therefore given by W = ∇2zzL(x, s, λh, λg) = ∇2 xxL(x, s, λh, λg) 0 0 µS−2  , (3.10)

and ∆ represents the radius for a trust region. In the trust region, the model and constraints are considered good approximations of the original problem. Sometimes problems can occur with the constraints in (3.8) because of the step d. Steps might lie outside the trust region ball in order to satisfy the constraints. One way of assuring this does not happen is with the usage of a vertical step, v. The step v is taken within the trust region and will never leave the region. The step is taken such that it tries to satisfy the constraints. A parameter ξ ∈ (0, 1) is introduced and v can be formulated as

min k ˆA(z)Tv + c(z)k2

subject to kvk2≤ ξ∆.

As we can see from the objective function, v is a measure of how well the constraints will be fulfilled. With v defined we can modify (3.8) into

min ∇ϕ(z)T_{d +}1

2d

T_{W d}

subject to ( ˆA(z)Td = ˆA(z)Tv, kdk2≤ ∆.

(3.12) This modification allows us to guarantee feasibility. A trivial solution can always be d = v which automatically will be inside the trust region, while satisfying the constraints to some degree. To ensure that the steps are taken in good directions, a merit function is introduced. The merit function should decrease when a new step is taken, otherwise the step is rejected and a new one is required. The merit function is defined as

φ(z, ν) = ϕ(z) + νkc(z)k2 (3.13)

where ν > 0 is another penalty parameter. It is chosen according to how much each step should decrease the merit function.

The previous calculation of step d was done with sequential quadratic pro-gramming. However, one can also do this with Newton’s method applied to the optimality conditions. The approach in this report is therefore tightly connected with traditional methods and the results have strong similarities. To visualise this, the calculation of step d is now done with Newton’s method. The starting equations are the Karush-Kuhn-Tucker conditions

    ∇f (x) + Ah(x)λh+ Ag(x)λg −µS−1_{e + λ} g h(x) g(x) + s     = 0. (3.14)

Newton’s method is then applied to the KKT-conditions to obtain a iteration where the step size d, at iteration k, satisfies

    ∇2 xxLk 0 Ah(xk) Ag(xk) 0 µS_k−2 0 I AT h(xk) 0 0 0 AT g(xk) I 0 0         dx ds dλh dλg     = −     ∇f (xk) + Ah(xk)λh+ Ag(xk)λg λg− µS−1e h(xk) g(xk) + sk     . (3.15)

3.1. Interior point algorithm 15

Rewriting this with λ+

h = λh+ dλh and λ + g = λg+ dλg gives     ∇2 xxLk 0 Ah(xk) Ag(xk) 0 µS_k−2 0 I AT h(xk) 0 0 0 AT g(xk) I 0 0         dx ds λ+_h λ+ g     = −     ∇f (xk) −µS−1_e h(xk) g(xk) + sk     . (3.16)

The solution of this system of equations gives a value for d that is also obtained from the sequential quadratic programming approach. Under the assumption that the problem in (3.8) is strongly convex, it has been shown that the solu-tion aligns with the one from (3.16). The convexity is equivalent to W being positive definite on the null space of ˆA(z)T_{. The connection with traditional}

methods motivates the relevance of this method. Byrd, Gilbert and Nocedal [2] emphasises the similarities in local behaviour and the differences in global behaviour between the alternatives. Since the similarities rely on the problem being strongly convex this method differs from other interior point methods. So far we have written the steps in the main algorithm for solving the opti-mization problem. However, some of the parameters that are used are unknown and need to be calculated every iteration. The Lagrange multipliers are such parameters and are found with a least square approach. The equation follows as λk = λh λg  = λLS(xk, sk, µ) = ( ˆAkTAˆk)−1AˆTk −∇f (xk) µe  , (3.17) where ˆ Ak = Ah(xk) Ag(xk) 0 Sk  . (3.18)

Also the trust region needs to be determined at each iteration. The trust re-gion should ensure that the approximations of the model and the constraints hold and that the slack variables are positive. This is to guarantee feasibility. The first criteria needs a simple bound such as k(dx, ds)k ≤ ∆. This not quite

enough, since the shape of the trust region is of importance as well. It is not desirable that the slack variables approach zero too soon. This is accomplished by punishing dswhen near boundaries. The result of this is k(dx, Sk−1ds)k2≤ ∆.

Here S−1

For the second criteria we use an equation that origins from the fraction to the boundary rule in the form of

sk+ ds≥ (1 − τ )sk =⇒ ds≥ −τ sk, (3.19)

where τ ∈ (0, 1) is a predetermined parameter. This equation allows us to formulate the trust region as

ds≥ −τ sk and k(dx, S_k−1ds)k2≤ ∆. (3.20)

The barrier parameter µ will approach to zero with each iteration but exactly how it will reduce is not obvious. Both [2] and [3] use advanced methods to calculate µ for each iteration. These methods introduces further parameters just for this purpose. However, other articles [16] have shown that a much easier method performs as well as the more complicated one. If the last barrier problem (3.3) was solved sufficiently enough in less than three iterations, we let

µk+1= µk/100. (3.21)

If it takes three iterations or more we instead let

µk+1= µk/5. (3.22)

To conclude the section about the algorithm a summary is given. It is divided into steps that the algorithm proceeds through:

1. Set the values of the parameters ξ ∈ (0, 1), τ ∈ (0, 1) and ν > 0. Set initial values of µ > 0, ∆ > 0, x and s.

2. If problem (3.2) is solved good enough with regards to a stopping tolerance, stop.

3. If problem (3.3) is solved good enough with regards to a stopping tolerance, go to 6.

4. Calculate λ, W , ˆA(z)T. Also compute the vertical step v from (3.11). 5. Calculate the step d by solving (3.12).

6. If the step d reduces the merit function φ sufficiently, update x and s-values and go to 3. Otherwise, decrease the trust region and go to 4. 7. Decrease µ and go to 2.

3.2. Geometrical optics implementation 17

3.2

Geometrical optics implementation

This section describes how the geometrical way of modelling the wave propaga-tion was implemented. It also accounts for how the optimisapropaga-tion was integrated with the implementation, such as what the objective function was and so on. The problem formulation included a radiating point source and a lens. The model was implemented in two dimensions but can easily be generalised to three dimensions. This generalisation will be explained later on. Since this is a geometrical implementation the source was radiating rays. Ideally, the number of rays would be infinite as that would be the closest to a real source. This was in practice not realisable because of computational times. Instead a finite number N was chosen for the number of rays. We assumed that the rays radi-ated from the point source in a spherical pattern, resembling a spherical wave. The goal for the lens was then to convert this spherical wave into a plane wave. Figure 3.1 shows how the intended set-up for the implementation looked like.

Figure 3.1: Setup for the geometric implementation. The dot represents the point source and the rectangle the lens.

Since the goal for the lens was to return a plane wave, the outgoing rays should be as parallel as possible, pointing forward. Another way of seeing this was that the rays should leave the lens with a small angle between themselves and the normal of the lens. The objective function was constructed with this in mind. It takes a permittivity profile for the lens as input and returns a vector of the angles the rays have (relative to the normal) when leaving the lens. This vector was then minimised with the interior point algorithm. To calculate the departure angles the function needed to propagate each ray through the lens according to

the given permittivity profile. The propagation was based on repeated usage of Snell’s law and ray tracing. Some consideration was also made to ensure that not all of the energy in the rays were reflected. Fresnel’s equations were used for these calculations.

Two different ideas were implemented to try to achieve an output in form of a plane wave. The first alternative discretised the lens both vertical and hori-zontal. Every cell could then take an arbitrary number as its permittivity. This version will therefore be called arbitrary distribution. The second alternative used a known distribution and varied the parameters of the distribution. This version will be called radial distribution. This is the more common approach to this problem, and has been well studied.

3.2.1

Arbitrary distribution

In the arbitrary distribution the lens is divided into a number of discrete cells where each one has an individual permittivity. Such a lens can be seen in Figure 3.2. Lens 11 m1 mn 12 21 1n 2n ... ... ... . . . . . . . . . 22 . . . m2

Figure 3.2: Overview of an arbitrary distribution where the lens is divided into mrows and n columns.

3.2. Geometrical optics implementation 19

Here, the parameters 11to mnformed an input matrix

 =      11 12 . . . 1n 21 22 . . . 1n ... ... ... ... m1 m2 . . . mn      . (3.23)

The values for the permittivity were not completely arbitrary since there are some constraints that needed to be upheld. A lower bound was introduced to make sure the chosen values actually were feasible to achieve in practice, i.e. larger or equal to one. An upper bound was also introduced because of the same reason. In mathematical form the constraints looked like

lb ≤ ij≤ ub i = 1, . . . , m, j = 1, . . . , n (3.24)

where lb and ub stands for lower and upper bound respectively. Because of the arbitrary approach, no constraints existed between different epsilon values.

3.2.2

Radial distribution

In the radial distribution we assume that the permittivity profile in the vertical direction follows a known distribution with a number of parameters. In the literature there exists many different distributions but they all share a gradient that decreases radially. Some initial tests showed that the difference between the different distributions was small when taking into account different values on the parameters. Mathematically the distribution is defined as

f (r) = A

cosh2(ωr) (3.25)

where A and ω are the parameters for the center value and the steepness of the slope, respectively. Two parameters were enough to create a wide range of distributions for the permittivity. The distribution of choice can be seen in Figure 3.3 where the parameters are set to different values.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Radius 2 3 4 5 6 7 8 9 10 Permittivity A = 10, = 1.1 A = 9, = 1.3 A = 8, = 0.8 A = 7, = 0.5 A = 6, = 1

Figure 3.3: Various distributions where the parameters are varied. Like in the arbitrary case the number of columns can be varied, meaning that the lens can consist of multiple layers with different distributions. The param-eters are also here restricted with constraints. To avoid the permittivity being less than one everywhere in the lens while also being reasonably small similar constraints as in (3.24) were formulated.

3.3

Physical optics implementation

The following section presents how the modelling based on physical optics was implemented. It explains how the wave is propagated through the lens and what the objective of the optimisation is.

In the physical optics implementation there are both similarities and differ-ences compared to the geometrical implementation. The problem formulation depicted in Figure 3.1 is the same except for one thing. Instead of the source radiating rays, it is now radiating a wave towards the lens. The lens is then used to convert the incoming spherical wave to a plane wave. Also here the implementation is done in two dimensions but the same generalisation for three dimension holds. Since the problem now will consist of waves a frequency is needed to be set. The frequency was chosen to be 10 GHz, relatively high such that the approximations hold.

3.4. Simulation and verification 21

Beginning from the source, a spherical wave will be generated. Similar to the geometrical case, the wave will also be discretised due to computational times. The discretisation can be seen as looking locally at the wave at N different positions. The desired output of the lens is a plane wave and a property of a plane wave that is simple to deal with is phase. In a plane wave, the phases in different positions of the wave front should be the same. The objective function in this case will therefore be that all N measurement points should have phases as close to each other as possible. Initially the phase for all measurement points will be the same. However, since they are leaving the source in a spherical pat-tern they have different distances to the lens. This leads to a difference in phase for the measurement points already at the boundary of the lens. The lens will then try to minimise the differences by delaying the phase for some points. The delay comes from the permittivity profile of the lens.

The crucial part of this implementation is naturally how the wave and fields propagate through the lens. As the name implies the foundation for this im-plementation was physical optics. Physical optics and its assumptions let us use the propagator for calculating the fields throughout the lens. To be able to use the propagator we first needed to find the permittivity values for each mea-surement point. The initial point of contact with the lens and the permittivity profile was known. An assumption was then made that each measurement point travelled at that same height through the whole lens. In a similar work [9] the same assumption was used with good results.

In the physical optics implementation, only the approach with a radial dis-tribution was studied. This was done because of initial results from geometrical optics showed that a radial distribution was superior to an arbitrary. Like in the geometrical case the number of layers with a distribution could vary. The same constraints on the parameters in the distributions as in the geometrical case are also applied here.

3.4

Simulation and verification

The previous sections contain information about the different models and how they are implemented. In order to compare these to each other in a just way, a way of quantifying them was needed. The trivial choice was to simulate the models in a full-wave modelling tool. This let us not only compare the lens models between themselves but also see how valid the approximations were. The simulations were done in the software CST Studio Suite. CST is a

compu-tational electromagnetics tool that can solve Maxwell’s equations and visualise the fields movement. There are endless alternatives and choices in the software but the results that were most interesting in this work were power flow and the radiating pattern. Power flow is the equivalence to the poynting vector and will give a measure of how good focus the lens can achieve. The radiating pattern will show the directivity of the energy. Both are of significant importance to a lens and are therefore good ways to assess it. One can also choose which solver CST will use to solve the problem. Both the integral solver and the time domain solver were tested but because of time complexity the time domain solver was preferred. The time domain solver is based on Finite-difference time-domain method which is a numerical analysis technique.

CST solves and visualises the problems in three dimensions but our optimi-sation and modelling was done on a two-dimensional problem. The step from two dimensions to three is however easy. Since we have symmetry around one axis, in this case the x-axis, we can see the two-dimensional lens as a vertical cross section in the middle of a cylinder. In Figure 3.4 the situation is illustrated. The red dotted rectangle is the modelled lens and the cylinder is the correspond-ing three-dimensional lens. The middle value of the lens will correspond to the radially innermost part of the cylinder. The cylinder will then obtain the rest of its values by rotating the cross section plane, making it rotationally symmetric.

Figure 3.4: Illustration of the 3D-generalisation.

In the setup the cylindrical lens was in between the source and the receiver. It is more straightforward to have a source emitting plane waves rather than spherical waves in CST so such a source was chosen. Due to reciprocity this was not an issue. The source and the receiver simply changed places and the cylinder was flipped. The flip is needed since there is no symmetry around the y-axis of the lens. To be able to measure and see the power flow, some sort of receiver is needed. A simple vacuum plane was introduced just for visual purposes. The plane begins immediately after the lens such that the power flow can be calculated everywhere on the plane. The setup can be seen in Figure 3.5

3.4. Simulation and verification 23

where the plane wave approaches the lens from the left with an incidence angle of zero. The lens is in this case multilayered and can be seen to the left in the picture. The rectangle is the vacuum plane which will be used as a receiver.

Chapter 4

Results and discussion

This chapter presents the computational results acquired in this thesis work. It also presents a comparison between the acquired models and two models of the works described in Section 2.6. The purpose of this section is to illustrate the results in a perspicuous way such that the analysis can be easily followed. Some discussion is also made around the results and how to interpret it.

During the optimisation it was noticed that the algorithm sometimes fixated on a local minimum and hence did not find a global solution. To avoid this, the goal function was evaluated with a large amount of random values. The one that gave the best result was then chosen as the initial value in the optimisation algorithm. This gave us a much better chance of finding the global minimum, although not guaranteed.

The parameters in the optimisation algorithm were chosen in line with [3]. τ ∈ (0, 1) and ξ ∈ (0, 1) were set to 0.995 and 0.8 respectively. The penalty parameters µ and ν were set to 0.1 and 1. Lastly the value for the trust re-gion radius ∆ was also set to 1. The optimisation was performed in MATLAB. The lens also had some parameters to decide, the width and the radius. In the geometrical implementation these values did not matter since the calculations were independent of what unit the distances were. However the physical optics implementation as well as the simulation program CST took these values into consideration so they were needed to be set. For the approximations to hold the size of the structures should be larger than the wavelength. Since we chose a frequency of 10 GHz the wavelength was close to 3 cm, calculated with the speed of light c ≈ 3 · 108_{. The values for the width and radius were then set to}

5 cm and 15 cm, respectively.

This section will be divided into three parts: one where the results from the implementation is shown, one for the results of the simulations and one where two lens models from Isakov [9] and Zhang [17] are presented for comparison.

4.1

Implementation results

4.1.1

Geometrical optics

In the arbitrary distribution the number of rows and columns had to be cho-sen beforehand. A few different combinations were tested initially with mixed results. In this section two of the combinations will be displayed. In Figure 4.1 the results of these distributions can be seen. The blue lines are rays that are being traced through the lens. In this implementation the lens had a ratio between the radius and the width that was not the same as the decided relation in the previous section. When this model later was simulated its radius and width was aligned with the predetermined dimensions.

(a) Arbitrary permittivity dis-tribution with 16 rows and 12 columns.

(b) Arbitrary permittivity dis-tribution with 80 rows and 40 columns.

Figure 4.1: Two arbitrary distributions with 20 rays

As we can see in the figure the results are lacklustre as the rays are far from leaving the lens with a normal incidence. Notably is that the result is signifi-cantly better in the case with a larger discretisation than the other one. This is intuitive since it increases our chances of steering the rays in the correct

direc-4.1. Implementation results 27

tions. Another important detail to note is the large amount of cells that are the same colour, i.e. cells that have the same permittivity. The rays will not pass through every cell and because of that some of them will be redundant. For the cells in the beginning of the lens this is due to the chosen contact points for the rays. They were chosen such that no ray left the lens vertically. The unused cells will then keep their initial values causing many cells to not be optimised. This is an effect of the number of rays and, as previously stated, ideally we would want to have an infinite amount of them. Some efforts were made with interpolation to make all cells have a somewhat optimal value but the results were barely improved. These poor results indicated that the arbitrary approach was not the best approach and focus shifted instead over to the radial approach. In Figure 4.2 two radial implementations can be seen. One has 1 layer with a radial distribution and the other one has 8 layers. The lens width is always the same, only the thickness of each layer changes. The rays are now almost perpendicular to the lens while leaving it. This is a clear improvement compar-ing to the arbitrary approach. In this case the permittivity is a little bit higher but not by a large margin. An important detail to note is that the lens with 8 layers does not perform better than the one with 1 layer. Actually, in the optimisation, all 8 layers had close to identical values on the parameters in the distribution. One reason for this is that a transition between two layers with a large discrepancy in permittivity will lead to a large change in angle according to Snell’s law. This large change can then be troublesome for the second layer to diminish while still maintaining the radial distribution. This property is likely a consequence of the geometrical approach.

(a) Lens with 1 radial layer. (b) Lens with 8 radial layers

This distribution is not as dependent on the number of rays since all permittivity values still will have optimal values according to the radial distribution. It is also computationally faster since the number of optimisation variables decreases with the introduction of a known distribution. However, a downside is that we limit ourselves within the boundaries of the distribution. This can in some cases be an issue if the desired permittivity profile is complex. The distribution would then have a hard time to follow and replicate such a profile.

4.1.2

Physical optics

Lastly, the results for the physical optics implementation is presented. Since this implementation was not based on rays like the others the visualisation of the result is completely different here. The optimisation was in this case based on the phase of different measurement points of the wave front. To compare this implementation with the others we will instead look at the results from the simulations.

Like in the geometrical implementation of the radial distribution, we will look at one lens with 1 layer and one lens with 8 layers. Figure 4.3 shows how the wave front looks in the two cases. The number of measurement points was in these cases chosen to be 20.

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 y-axis of lens (m) 0 50 100 150 200 250 300 350 Phase (degrees)

(a) Phase of the wave front after the lens with 1 layer.

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 y-axis of lens (m) 0 50 100 150 200 250 300 350 Phase (degrees)

(b) Phase of the wave front after the lens with 8 layers.

Figure 4.3: Phase plots

The ideal result is for the phase to be equal at all measurement points of the wave. With this in mind it is clear that the lens with 8 layers performs

signifi-4.2. Simulation results 29

cantly better than the one with 1 layer. Unlike in the geometrical approach the layers’ parameters now take different values instead of the same. Both results are good with small deviations in phase between the measurement points. The permittivity values in the single layer begins at 10.5 and drops to 4. In the lens with 8 layers the permittivity values varies between 20 and 2. The values are comparable to the ones in the geometrical cases. The parameter values for models with 1 layer and 8 layers is presented in Table 4.1 and Table 4.2.

Table 4.1: Parameters for the lenses with 1 radial distribution

Model A ω

GO 20.2 0.92

PO 10.6 1.05

Table 4.2: Parameters for the lenses with 8 radial distribution

Model A1 A2 A3 A4 A5 A6 A7 A8 ω1 ω2 ω3 ω4 ω5 ω6 ω7 ω8

GO 19.4 19.4 19.4 19.4 19.4 19.4 19.4 19.4 1.10 0.94 0.94 0.95 0.96 0.97 0.97 0.98 PO 6.0 20.0 13.8 7.7 6.3 20.0 20.0 6.0 1.10 0.93 0.77 1.10 1.10 0.85 0.91 1.10

4.2

Simulation results

In this section the simulation results from CST will be presented. To be able to interpret the results we need a reference. The reference in this case will be a simulation where the material in the lens was replaced with vacuum. The power flow and radiating pattern of the vacuum lens is displayed in Figure 4.4. Figure 4.4a shows the Poynting vector over the lens and the vacuum plane. Since the lens is made out of vacuum the power flow is close to constant everywhere as it should. For our lenses we want to have high values at the focal point and close to zero otherwise. Figure 4.4b shows the normalised radiating pattern. The most interesting thing in these plots is how big the reflection lobes are relative to the main lobe. The main lobe is at 0 degrees which is straight forward, and the reflection lobes are at −180 and 180 degrees which is straight backwards. Since the lens is made out of vacuum the reflections are low in this case. The size of the side lobes is an indication of the focus abilities of the lens. However, this graph is based on the fields a large distance from the lens. We are more interested in

what happens close to the lens where it should focus, which will be clear when we look at the power flow. The reflection lobes are still interesting even if they are taken some distance away from the lens. Note that the magnitude on the y-axis is normalised and is in dB.

(a)

(b)

Figure 4.4: Simulations of lens made of vacuum.

To get a good overview and clearly see the differences between the models all radiating pattern plots will be presented first. After this all power flow plots will be presented in the same way. Figure 4.5 contains the images of the radiating patterns of all models. As we can see the structure of all graphs are similar but they differ in the magnitude of the side lobes and the reflection lobe. The

4.2. Simulation results 31

size of the side lobes indicate how well the lens can focus the fields and the reflection lobes tell us how much of the fields that reflect instead of transmit. Both arbitrary distributions have a distinguished main lobe and a reflection around 8 dB lower than the main lobe. This translates to roughly 40% of the amplitude of the main lobe.

(a) Small arbitrary permittivity distribution.

(b) Large arbitrary permittivity distribution.

(c) Geometrical optics radial dis-tribution with 1 layer.

(d) Geometrical optics radial dis-tribution with 8 layers.

(e) Physical optics radial distribu-tion with 1 layer.

(f) Physical optics radial distribu-tion with 8 layers.

Figure 4.5: Radiating patterns for all models

The geometrical optics models have even higher reflection values at 6 dB lower, resulting in a reflection amplitude of 50% of the main lobe. Their side lobes

are low compared to the main lobe. The physical optics models have reflection lobes of 8 dB lower which is the same as for the arbitrary distributions. They do not have an equally distinguished main lobe as the other models. Overall the models show high values of reflection which trivially is not desirable. We transmit less power and the reflections can cause disturbances if they begin to reflect multiple times. One example of this is if the reflected wave once again reflects but now on the source of the wave, causing another wave front towards the lens. This is troublesome for the receiver since those reflections will make the signal noisy. The receiver can in those cases struggle to interpret the signal. The power flow of all the models can be seen in Figure 4.6 in the same or-der as in previous figure. Compared to the vacuum lens we now see a clear spot where the lenses try to focus the fields. Most of the plots show good focusing abilities of the lens but the magnitude of the fields differ. The outlier is the arbitrary distribution with a larger discretisation which barely can focus the fields. This is due to the large amount of cells that were unoptimised and hence had the same values. This makes the power flow have the same characteristics as the vacuum lens. The small arbitrary distribution had an overall low mag-nitude in its power flow. The focus is lacking since the power flow does not converge. With that both arbitrary distributions perform badly which aligns with the results from the implementation.

The geometrically optimised gradient models are naturally similar since the parameters almost were identical. They show good focus in the focal point and have a magnitude that is four times greater than the one in the arbitrary ones. The one with 8 layers perform slightly better due to the small changes in the parameters. Lastly the physical optics optimised gradient models both have resembling structures but the scale is important to note. The model with one layer is slightly worse than the two geometrically optimised gradient models in terms of magnitude but have a similar focus. The one with 8 layers however have a magnitude that is around 7 times greater than the geometrical instance with 8 layers. In regards to magnitude it is outstanding compared to the other models. Its focus is also good which makes for a strong overall impression.

4.2. Simulation results 33

(a) Small arbitrary permittivity distribution.

(b) Large arbitrary permittivity distribution.

(c) Geometrical optics radial dis-tribution with 1 layer.

(d) Geometrical optics radial dis-tribution with 8 layers.

(e) Physical optics radial distribu-tion with 1 layer.

(f) Physical optics radial distribu-tion with 8 layers.

4.3

Comparison

To get a perception of how good the acquired models were, a comparison was made with other similar works. Those lens models then acted as references and could hint if our models were performing well or not. One lens from Isakov’s pa-per [9] was chosen and one from Zhang [17]. In respective papa-per the dimensions of the lenses were not the same as in this work. In order to get a fair comparison both models were then transformed into the right dimensions. Small deviations in focal length are present. They were then simulated in the same manner as previously. In Figure 4.7 the radiating patterns can be seen for both models. Most notably are the reflection values where the Zhang-lens is better than the Isakov-lens and performs similar to the 8-layered lens that was modelled with physical optics.

(a) Isakov (b) Zhang

Figure 4.7: Radiating patterns for the Isakov-lens and the Zhang-lens

No major conclusions could be drawn from the radiating pattern so the power flow was instead studied, see Figure 4.8. In the Isakov-lens we can see the shortened focal length where the focus has shifted slightly to the left. However, the magnitude of the focus is high and is only second to the 8-layered lens that was modelled with physical optics. The Zhang-lens shows convergence in power flow but the magnitude is low compared to the best models. The strong impression made by the 8-layered lens that was modelled with physical optics is strengthened after this comparison. It has a better power flow and a better radiating pattern than both of the two models used in the comparison.

4.3. Comparison 35

(a) Isakov (b) Zhang

Chapter 5

Conclusions and outlook

In this thesis we investigated how a flat dielectric lens can be used to focus ra-dio waves instead of a traditional convex lens. The dielectric lens was optimised with two different computational methods: geometrical optics and physical op-tics. The lenses were optimised with regards to their focusing abilities. This section aims to compile the thesis and present the conclusions that were drawn. From the results it is clear that the lenses are capable of transmitting and focusing electromagnetic waves at a certain distance from the lens. The two arbitrary distributions had the worst power flow which is seen as the most im-portant result in this work. They also had similar values on the reflection as the others and is therefore seen as the worst performing model. The geometri-cal approach with radial distributions had higher magnitude in its power flow making it a stronger choice. The physical approach with radial distributions did even better with the highest magnitude and a good focus. The radial distribu-tions are also advantageous since the computational times in the optimisation is significantly lowered. Also in comparison to the Isakov-lens and the Zhang-lens our lenses modelled with physical optics showed good results. The downside with all implemented models is that the reflections are quite large. One way to handle the reflections could be that the receiver only is active in a short time window, timing the original signal and avoiding the delayed reflections.

The optimisation algorithm of choice was an interior point method. The reason-ing behind this was that the interior point methods are considered to be strong and efficient when handling non-linear problems. Much of the methods strength lies in how it handles non-linear constraints in a simple way. The problem for-mulation in this work did however not have any constraints except for upper

and lower bounds. The lack of constraints could let less advanced algorithms be used with similar results but faster computational times. One small downside with the interior point method was that there was no guarantee of finding a global minimum when optimising. A guarantee of finding the global optima every time is extremely hard and is not something that is easily achieved. In-troducing some heuristic would be a potential way of improvement.

A natural next step could be to actually produce one of the designed lenses and see how it performs in practice. This would indicate whether the designed lens could be useful in some application. In future similar works one could explore a wider range of distributions. The distribution that performed best in this thesis might not be ideal which would be interesting to immerse in. One approach could be to employ a high-degree polynomial and then optimise the coefficients in the polynomial. To compare different optimisation methods against each other would also be an interesting task. Especially when the num-ber of variables increase it is important to have a solver that can find a good solution in a reasonable time space. A future work could also study other com-putational methods to model the electromagnetic wave such as transformation optics. Lastly, it may be of interest to implement some sort of multi-objective optimisation where we not only look at maximising focusing abilities but also minimising the reflection.

Bibliography

[1] A.C. Balanis. Advanced engineering electromagnetics. John Wiley & sons, 1989.

[2] R.H. Byrd, J.C. Gilbert, and J. Nocedal. A trust region method based on interior point techniques for nonlinear programming. Mathematical pro-gramming, 89(1):149–185, 2000.

[3] R.H. Byrd, M.E. Hribar, and J. Nocedal. An interior point algorithm for large-scale nonlinear programming. SIAM J. Optim., 9(4):877–900, 1999. [4] D.K. Cheng. Field and Wave electromagnetics. Pearson Educated Limited,

2014.

[5] I.I. Dikin. Iterative solution of problems of linear and quadratic program-ming. Dokl. Akad. Nauk SSSR, 174(4):747–748, 1967.

[6] J. Gutiérrez-Meana, J.A. Martínez-Lorenzo, and F. Las-Heras. Electromag-netic wave propagation in complex matter: High Frequency Techniques: the Physical Optics Approximation and the Modified Equivalent Current Ap-proximation (MECA). InTech, 2011.

[7] R.F. Harrington. Time-harmonic electromagnetic fields. McGraw-Hill Book Company, 1961.

[8] Y. He and G.V. Eleftheriades. A highly-efficient flat graded-index dielectric lens for millimeter-wave application. In 2017 IEEE International Sympo-sium on Antennas and Propagation USNC/URSI National Radio Science Meeting, pages 2655–2656, 2017.

[9] D. Isakov, C.J. Stevens, F. Castles, and P.S. Grant. 3d-printed high dielec-tric contrast gradient index flat lens for a directive antenna with reduced dimensions. Adv. Mater. Technol., 1(16), 2016.

[10] S. Jain, M. Abdel-Mageed, and R. Mittra. Flat-lens design using field transformation and its comparison with those based on transformation op-tics and ray opop-tics. IEEE Antennas and Wireless Propagation Letters, 12:777–780, 2013.

[11] G. Jönsson. Våglära och optik. Teach Support, 2016.

[12] I.V. Lindell. Methods for electromagnetic field analysis. Oxford University Press, 1992.

[13] R.K. Luneburg. Mathematical theory of optics. Providence R.I., 1944. [14] R.A. Paulus. The lorentz reciprocity theorem as a consistency test for

propagation models. In Multiple Mechanism Propagation Paths (MMPPs): Their Characterisation and Influence on System Design. Papers presented at the Electromagnetic Wave Propagation Panel Symposium, held in Rot-terdam, The Netherlands, 1992.

[15] S. Rikte, G. Kristensson, and M. Andersson. Propagation in bianisotropic media - reflection and transmission, volume TEAT-7067 of Technical Report LUTEDX. 1998. Published version: IEEE Proc. - Microwaves, Antennas and Propagation, 148(1), 29-36, 2001.

[16] R. Waltz, J. Morales, J. Nocedal, and D. Orban. An interior algorithm for nonlinear optimization that combines line search and trust region steps. Mathematical programming, 107(3):391–408, 2006.

[17] S. Zhang. 3d printed dielectric fresnel lens. In 2016 10th European Con-ference on Antennas and Propagation (EuCAP), pages 1–3. IEEE, 2016.

Linköping University Electronic Press

Copyright

The publishers will keep this document online on the Internet – or its possible replacement – from the date of publication barring exceptional circumstances.

The online availability of the document implies permanent permission for anyone to read, to download, or to print out single copies for his/her own use and to use it unchanged for non-commercial research and educational purpose. Subsequent transfers of copyright cannot revoke this permission. All other uses of the document are conditional upon the consent of the copyright owner. The publisher has taken technical and administrative measures to assure authentic-ity, security and accessibility.

According to intellectual property law the author has the right to be men-tioned when his/her work is accessed as described above and to be protected against infringement.

For additional information about the Linköping University Electronic Press and its procedures for publication and for assurance of document integrity, please refer to its www home page: http://www.ep.liu.se/.

Upphovsrätt

Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare – från publiceringsdatum under förutsättning att inga extraordinära omständig-heter uppstår.

Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner, skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för ickekommersiell forskning och för undervisning. Överföring av upphovsrätten vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av dokumentet kräver upphovsmannens medgivande. För att garantera äktheten, säkerheten och tillgängligheten finns lösningar av teknisk och administrativ art. Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i den omfattning som god sed kräver vid användning av dokumentet på ovan beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan form eller i sådant sammanhang som är kränkande för upphovsmannens litterära eller konstnärliga anseende eller egenart.

För ytterligare information om Linköping University Electronic Press se för-lagets hemsida http://www.ep.liu.se/.