close

Вход

Забыли?

вход по аккаунту

?

An improved adjoint method for design of microwave devices with three-dimensional finite elements

код для вставкиСкачать
INFORMATION TO USERS
This manuscript has been reproduced from the microfilm m aster. UMI films
the text directly from the original or copy submitted. Thus, som e thesis and
dissertation copies are in typewriter face, while others may be from any type of
computer printer.
The quality o f th is reproduction is d ep en d en t up o n th e quality o f th e
copy subm itted. Broken or indistinct print, colored or poor quality illustrations
and photographs, print bleedthrough, substandard margins, and improper
alignment can adversely affect reproduction.
In the unlikely event that the author did not send UMI a complete manuscript
and there are missing pages, these will be noted.
Also, if unauthorized
copyright material had to be removed, a note will indicate the deletion.
Oversize materials (e.g., maps, drawings, charts) are reproduced by
sectioning the original, beginning at the upper left-hand com er and continuing
from left to right in equal sections with small overlaps.
Photographs included in the original manuscript have been reproduced
xerographically in this copy.
Higher quality 6? x 9" black and white
photographic prints are available for any photographs or illustrations appearing
in this copy for an additional charge. Contact UMI directly to order.
ProQ uest Information and Learning
300 North Z eeb Road. Ann Arbor. Ml 48106-1346 USA
800-521-0600
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
An Improved Adjoint Method
for Design of Microwave Devices
with 3D Finite Elements
Hatem Akel
i
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
National Libraiy
of Canada
Bibliothdque nationale
du Canada
Acquisitions and
Bibliographic Services
Acquisitions et
services bibliographiques
395 WaKngton Street
Ottawa ON K1A0N4
Canada
395, rue Weffngton
Ottawa ON K1A0N4
Canada
YourOm W W iW m i a
O u r N o b w r U tn n e m
The author has granted a nonн
exclusive licence allowing die
National Libraiy of Canada to
reproduce, loan, distribute or sell
copies of this thesis in microform,
paper or electronic formats.
L?auteur a accorde une licence non
exclusive pennettant a la
Bibliotheque nationale du Canada de
reproduire, preter, distribuer ou
vendre des copies de cette these sous
la forme de microfiche/film, de
reproduction sur papier ou sur format
electronique.
The author retains ownership of the
copyright in this thesis. Neither the
thesis nor substantial extracts from it
may be printed or otherwise
reproduced without the author?s
permission.
L?auteur conserve la propriete du
droit d?auteur qui protege cette these.
Ni la these ni des extraits substantiels
de celle-ci ne doivenl etre imprimes
ou autrement reproduits sans son
autorisation.
0-612-64494-4
Canada
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Abstract
A numerical method is described for determining the derivatives o f the scattering
parameters o f an arbitrary 3D microwave device with respect to changes in model
geometry. The new technique is based on the adjoint variable method as applied to finite
element analysis, but does not require the separate calculation o f an adjoint solution. The
new technique is applied to the optimization o f the design o f microwave components,
using a derivative-based, quasi-Newton method.
The derivatives o f the scattering parameters o f three devices are computed and
compared to analytical or finite difference values. The agreement is very good. The three
devices are: a length o f short-circuited waveguide, an E-plane miter-bend, and a
waveguide transformer.
The short-circuited waveguide length is then optimized to achieve a specified phase of
the reflection coefficient. The miter-bend and the transformer are optimized for minimum
return loss. All optimizations reached the (local) optimum solution after only 4-12
computations o f the field.
ii
with permission of the copyright owner. Further reproduction prohibited without permission.
R6sum6
Une methode numerique pour calculer les derivees de la matrice de dispersion d?un
circuit micro-onde 3D en fonction des parametres geometriques est presentee.
Similairement a la methode des elements finis, la nouvelle technique fait appel a la
methode de la variable adjointe, mais ne requiert pas le calcul de la solution adjointe. La
nouvelle technique est utilisee conjointement avec une methode quasi-newton basee sur
les derivees pour concevoir des composantes micro-ondes optimales.
Les derivees des parametres de dispersion pour trois composantes sont ainsi calculees et
comparees, selon la disponibilite des resultats, a des valeurs analytiques ou a des resultats
de differences finies. Une tres bonne correlation entre les differents resultats a ete
observee. Les trois composantes s o n t: un guide d ?onde termine par un court-circuit, un
coude a onglet dans le plan E et un transformateur guide d ?onde.
Le guide d?onde court-circuite a ete optimise en fonction de la phase du coefficient de
reflexion. Les deux autres composantes ont, quand a elles, ete optimisees pour minimiser
les pertes par reflexion. Toutes les simulations ont atteint la solution optimale (locale)
apres seulement 4 a 12 iterations.
with permission of the copyright owner. Further reproduction prohibited without permission.
Table of Contents
Abstract .................................................................................................................................... ii
Resume ...................................................................................................................................... iii
Table o f Contents .................................................................................................................... iv
List o f Figures .......................................................................................................................... vi
CHAPTER 1 Introduction
..................................................................................................
l
CHAPTER 2 Computing the Scattering Parameters .......................................................... 6
2.1 Introduction ....................................................................................................................... 6
2.2 Maxwell?s Equation ......................................................................................................... 6
2.3 Boundary Conditions ....................................................................................................... 8
2.3.1 Dirichlet Boundary Condition for E ....................................................................... 8
2.3.2 Neumann Boundary Condition for E ..................................................................... 10
2.3.3 Absorbing Boundary Conditions ........................................................................... 10
2.3.4 Port Boundary Conditions ....................................................................................... 10
2.4 Weak Form o f the Problem ............................................................................................. 14
2.5 Finite Element Discretization ......................................................................................... 15
2.6 Scattering Matrix .............................................................................................................. 18
2.7 Propagation Modes ..........................................................................................................
19
2.6.1 Rectangular Waveguide .......................................................................................... 19
2.6.2 Circular Waveguide ................................................................................................ 21
CHAPTER 3 The Gradient o f Scattering Parameters .........................................................
3.1 Introduction ......................................................................................................................
3.2 The Adjoint Method .........................................................................................................
3.3 The Nature o f the Adjoint Method ................................................................................
3.4 Evaluating the Derivative o f [ K ] ....................................................................................
23
23
23
25
27
CHAPTER 4 Optimization o f Microwave Devices ............................................................ 29
4.1 Introduction ........................................................................................................................ 29
4.2 Cost Function o f Microwave Devices ........................................................................... 29
4.3 Conditions for a Local Minimum .................................................................................. 30
4.4 Multivariate Gradient Methods o f Minimization ......................................................... 31
4.5 Quasi-Newton Methods .................................................................................................. 33
4.6 Line Search ....................................................................................................................... 35
4.6.1 Establishing Bounds for the Search ....................................................................... 36
4.6.2 Locating a Minimum ............................................................................................... 37
CHAPTER 5 Computer Program Implementation .............................................................
5.1 Introduction .......................................................................................................................
5.2 Program Structures ..........................................................................................................
iv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
40
40
5.2.1 Quasi-Newton Optimizer .......................................................................................
5.2.2 Line Search ..............................................................................................................
5.2.3 Mesh Generator .......................................................................................................
5.2.4 Nodal-to-Edge Converter ......................................................................................
5.2.5 3D Finite Element Solver .......................................................................................
5.2.6 Remodeller ..............................................................................................................
5.3 Data Files .......................................................................................................................
5.3.1 Geo.dat .....................................................................................................................
5.3.2 2D.dat .......................................................................................................................
5.3.3 Input.dat ..................................................................................................................
5.3.4 Ports.dat ...................................................................................................................
5.3.5 Mats.dat ....................................................................................................................
5.3.6 Param.dat .................................................................................................................
5.3.7 Tracing.log ..............................................................................................................
5.3.8 Tets.dat .....................................................................................................................
5.3.9 Edges.dat ..................................................................................................................
5.3.10 Results.dat .............................................................................................................
40
41
42
42
42
43
43
43
43
44
44
44
44
44
45
45
45
CHAPTER 6 Results ............................................................................................................
6.1 Introduction .....................................................................................................................
6.2 Short-circuited Rectangular Waveguide .....................................................................
6.2.1 Verifying the Gradient o f Z S U ............................................................................
6.2.2 Optimization ...........................................................................................................
6.3 E-plane Miter Bend Waveguide ...................................................................................
6.3.1 Verifying the Gradient ..........................................................................................
6.3.2 Optimization ...........................................................................................................
6.4 Single Stage Waveguide Impedance Transformer .......................................................
6.4.1 Verifying the Gradient ...........................................................................................
6.4.2 Optimization ...........................................................................................................
6.5 Summary .........................................................................................................................
46
46
46
48
49
51
53
55
58
60
61
63
CHAPTER 7 Conclusion ......................................................................................................
7.1 Summary ..........................................................................................................................
7.2 Original Contributions ....................................................................................................
7.3 Suggestions for Further Work .......................................................................................
66
66
67
67
APPENDIX A: Matrix Assembly for 3D Edge Elements ...............................................
69
APPENDIX B: Quasi-Minimal Residual Method ............................................................
72
APPENDIX C: Derivations ..................................................................................................
74
References ...............................................................................................................................
83
V
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table of Figures
Figure 2.1
Figure 2.2
Figure 2.3
Figure 2.4
Figure 2.5
Figure 2.6
Figure 5.1
: Closed and Open Suspended Strip Line ......................................................... 9
: An N-port Microwave D evice............................................................................ 11
: n is a Unit vector Outward from the Device .................................................. 13
: 2D Elements ......................................................................................................
16
: 3D Elements ....................................................................................................... 17
: (a) Rectangular Waveguide (b) Circular W aveguide ................................... 20
: Representative Flow Chart o f the Relation
Between the Six Main Modules ......................................................................
41
Figure 6 .1 : Rectangular Waveguide Short-Circuited at a Distance L from the Port ...
47
Figure 6.2 : Error in % o f the Phase o f S n and its Gradient ............................................
49
Figure 6.3 : Error in the phase o f S// Versus Solution Number ....................................... 50
Figure 6.4 : Cost Function in the Vicinity o f the Optimum Solution .............................. 51
Figure 6.5 : E-plane Miter Bend ........................................................................................... 52
Figure 6.6 : Displacement o f the Mesh Nodes Located on the Surface
o f die Inclined W all .........................................................................................
54
Figure 6.7 : % Error in the Gradient o f |S//| o f the miter-bend example ......................... 55
Figure 6.8 : Miter-Bend at Two Different Positions .......................................................... 56
Figure 6.9 : Miter-Bend Position as a Function o f the Solution N um ber.......................... 56
Figure 6.10 : Cost (dB) Versus Solution N um ber................................................................ 57
Figure 6.11 : Return Loss o f the Optimum Miter-Bend .................................................... 58
Figure 6.12 : Single Stage Rectangular Transform er........................................................... 59
Figure 6.13 : % Error in the Gradient o f IS//I ..................................................................... 61
Figure 6.14: The Cost (dB) as a Function o f the Solution Number ................................. 62
Figure 6.15: Return Loss o f the Initial and the Optimized Geometry ............................ 63
vi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 1
INTRODUCTION
M any modem communications and radar systems employ microwave devices, e.g.
antennas, filters, polarizers, phase shifters, circulators, and diplexers. O ver the last
decade, the application o f high frequency devices has grown very rapidly, and this has
generated a tremendous demand for accurate and easy-to-use analysis and design
software.
The performance o f any microwave device can be described in terms o f a scattering
matrix. The scattering matrix is an 7VxN matrix relating the input and output waves at the
N ports o f the device. Although the scattering parameters o f some microwave devices can
be obtained as closed form expressions, which can be used to design and optimize them,
in general this is not possible. Some components can be modeled using a lumped-element
circuit. However many o f these circuits turn out to be poor approximations, especially for
wide-band applications. The finite element method (FEM) is one o f many field-based
numerical techniques which can analyze an arbitrarily-shaped device and obtain very
good results compared to measurements [1]. However, a drawback o f FEM even with
advanced computers and improved algorithms, is still its execution time. This problem is
worse when FEM is used in an optimization routine where many solutions are required.
Consequently, it is important to find ways to reduce the number o f solutions needed to
optimize a microwave device.
This dissertation describes a method o f calculating the derivative (or sensitivity) o f the
scattering parameters o f a microwave device with respect to changes in its geometry. The
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
calculations are obtained directly from the finite element solution o f a boundary value
problem for the electric field. These derivatives are then used in an optimization routine,
which aims to improve the design o f the device.
Lee and Etoh [2] demonstrated in their paper the adjoint method where, for an Af-port
problem, 2N solutions are required to calculate unlimited numbers o f (first) derivatives o f
all the scattering parameters. Garcia and Webb [3] presented formulae to calculate the
derivatives o f the admittance parameters o f 2D microwave models using only N
solutions. The derivatives o f the scattering parameters are derived from the derivatives o f
the admittance parameters. In this dissertation, new formulae are presented for the
derivatives o f the scattering parameters directly, and for 3D problems, using only N
solutions [4]. These derivatives are then used in an optimization scheme.
Optimization techniques can be classified into two main categories: stochastic, and
deterministic. Stochastic techniques search the space o f the unknowns randomly, seeking
a global minimum o f the cost (objective) function. Deterministic methods generally
follow a trajectory in search space, but often only find the local minimum that is closest
to the starting point. There are two subclasses o f deterministic methods: direct, and
gradient, depending on whether they exploit the derivative o f the cost function. The
reader can refer to [5] for an exhaustive review o f optimization techniques.
Stochastic methods are best used with simple problems which have analytical
solutions, or can be solved in very short o f time, because they require very large numbers
o f cost function evaluations. There are many techniques that can be classified as
stochastic, e.g. random search, evolutionary strategy, and genetic algorithms [6].
Rechenberg [7] was the first to use the evolution technique in 1973. Arndt et al.[8-18] in
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
their mode-matching papers used the evolution technique to optimize all kinds o f
microwave components, e.g. filters, polarizers, diplexers, couplers, and orthomode
transducers. This was a successful approach because in mode matching the gradients
cannot be extracted easily, but on the other hand the analysis time is very short. Zhang
and Weiland [19] combined the evolution method with deterministic gradient techniques
to optimize waveguide-to-microstrip transitions. Holland [20] was the first to use genetic
algorithms in 1992, and since then many papers have appeared on using genetic
algorithms for low and high frequency electromagnetic applications. John and Jansen
[21] used the genetic algorithm to design and optimize M(MIC) components, Jones and
Joines [22] used the genetic technique to design Yagi-Uda antenna, and Altshuler and
Linden [23] applied it to a loaded monopole.
Direct techniques are used when no gradient information is available, or when the
target function is discontinuous. Some direct methods rely on pattern finding approaches
such as the simplex method, Hooke and Jeeves method, and Rosenbrock pattern search.
Other direct methods seek conjugate directions, such as Powell?s method. Okoshi et al.
[24] used Powell?s method to design and optimize a hybrid ring directional coupler.
Miyoshi and Shinhama [25] used Powell?s method to design planar circulators, and Aloss
and Guglielmi [26] used it also to design and optimize microwave filters.
The direct methods are considered inefficient techniques, due to the fact that they do
not make use o f gradient (derivative) information. Examples o f gradient methods are:
steepest descent, Newton Raphson, quickprop, quasi-Newton, and gradient associated
conjugate direction (GACD). The earliest literature regarding the use o f gradients goes
back to Hadamard [27] in 1910. Lee et al. [28] used the steepest descent method to design
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ridged waveguides, and Rakanos [29] used the quickprop technique for low frequency
applications. Garcia and W ebb [3] [6] used quasi-Newton techniques to optimize
rectangular waveguide H-bends and circulators. Suzuki and Hosono [30] used the quasiNewton method to optimize the cross section o f a waveguide giving the m inim um
conductor loss for the dominant mode, and Lee and Itoh [2] used the quasi-Newton
methods to optimize waveguide-to-microstrip transitions. Zhang [31] [19] proposed the
GACD method which he claims behaves better than the quasi-Newton methods, for some
cases.
Chapter 2 introduces the basic electromagnetic formulae and boundary conditions
encountered in microwave problems. The weak form is given as an alternative to solving
the vector wave equation in terms o f the electric field. The finite element modeling is
then presented, and the formula used to calculate the scattering parameters is introduced.
In Chapter 3, the adjoint method is presented, and its main formulae are derived. Then
it is shown how the gradients o f the scattering parameters can be obtained with no
separate adjoint solution. Finally, we present the details o f calculating the gradients for
tetrahedral edge elements.
Chapter 4 presents a general cost function and its gradient with respect to a
geometrical parameter. Then, a brief description o f the mathematical background o f some
gradient methods o f optimization are reviewed. The chosen gradient method (BFGS) is
presented in detail.
The computer program implementation is outlined in Chapter 5. The different data
files used in the program are discussed briefly.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The results are given in Chapter 6. In the first example the accuracy o f the calculations
o f the gradient is verified, by comparison w ith the analytical solutions. In the second
example, a one-dimensional optimization at a single frequency is presented. In the last
example, a multidimensional optimization over a range o f frequencies is considered.
The dissertation is concluded in chapter 7. Some recommendations for future work in
the microwave dom ain are given.
5
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 2
Computing the Scattering Param eters
2.1 In tro d u ctio n
The main objective o f this chapter is to summarize the application o f the finite element
method to finding the network parameters o f microwave devices. Finite element methods,
or other numerical analysis techniques, are used because most real problems do not have
analytical or closed form solutions. The numerical analysis o f any microwave device can
be formulated into a set o f linear equations. This set o f linear equations can then be
solved directly or iteratively to get the field distribution, and consequently, the scattering
parameters. This chapter presents the electromagnetic formulae, and how they are
translated into a set o f linear equations, which provide the means to calculate the field
distribution. We also present the necessary formulae used to compute the scattering
parameters from the field solution.
In section 2.2 we give Maxwell?s equations, and in section 2.3 the most commonly
used boundary conditions are described. Using Maxwell?s Equations and the boundary
conditions the weak form o f the mathematical problem is constructed in section 2.4. In
section 2.5 the finite element modeling is described. In section 2.6 the scattering matrix
and the formulae used to calculate it are presented, and in section 2.7 we discuss some
aspects of the theory o f modes in guided structures.
2.2 M axwell?s E q u a tio n s
Electromagnetic phenomena are governed by a set o f differential equations known as
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Maxwell?s Equations. For a wave travelling in a free space or in a guided device, with
linear, isotropic, sourceless, and homogeneous or inhomogeneous media, the electrical
(E) and magnetic (H) field must satisfy the following four relationships [32]:
? ?
VX E =
5H
a
-L L ------
(2.1a)
V x H = г - + ( tE
dt
(2'lb)
V -e ?E = 0
(2.1c)
V -/* H = 0
(2. Id)
where <r, fx and e' are the conductivity, permeability and permitivity o f the media of
concern, respectively. When a time harmonic electromagnetic wave is assumed, complex
phasors can be introduced, i.e.,
E(x,y,z, t) = Re{E(x,y,z)eJ'QJ' }
(2.2)
where co is the angular frequency o f the wave.
Maxwell?s equations can then be rewritten as:
V x E = -jcojuH
(2.3a)
V x H = jcoe E
(2.3b)
V eE = 0
(2.3c)
VjiH =0
(2.3d)
where e is the complex permitivity:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
With some mathematical manipulation, one can combine the four formulae o f (2.3) to
obtain one formula with E or H as the unknown. The new formula is called the Curl-Curl
Equation:
V x ? V x E - J t :s E = 0
(2'5)
M r
where k is the free space wave number:
k =c
(2 6 )
c is the speed o f light in free space. /jr and e T are the relative permeability and
permitivity o f the media. A similar formula for the H-field could be derived. However, in
this thesis, only the E-field formula (2.5) is considered.
2 .3 B o u n d a ry C o n d itio n s
The finite element method can only analyze finite volume models. Consequently,
whether we are solving a guided structure, or an antenna problem, the device must be
enclosed by a boundary (Figure 2.1).
Four types o f boundary condition are commonly encountered in guided microwave and
antenna problems: Dirichlet, Neumann, Absorbing, and Port.
2.3.1 Dirichlet Boundary Condition for E
This condition constrains the tangential component o f the E-field over the surface:
E x n = E0
(2.7)
where n is the unit vector normal to the surface. If Eo is equal to zero, we obtain the
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
homogeneous Dirichlet boundary conditions. This case is imposed on surfaces which are
touching perfectly conducting materials. The Dirichlet Boundary Conditions can also be
applied on some symmetrical planes where the tangential components o f the E-field, due
to symmetry, must be equal to zero.
Ca)
F ictitio u s
" -\B ou n d ary
Figure 2.1: (a) Closed and (b) Open Suspended Strip Line. In (b) a fictitious boundary was introduced to
truncate the infinite Space.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.3.2 Neumann Boundary Condition for E
This condition constrains the tangential component o f the curl o f E over the surface:
Vx Ex n = 0
(2.8)
H xn = 0
(2.9)
which also implies that:
i.e., the tangential part o f the magnetic field vanishes. This condition is satisfied on some
symmetrical surfaces; i.e. surfaces that split a microwave device into two parts, that are
geometrically and magnetically symmetric.
2.3.3 Absorbing Boundary Conditions
Absorbing boundary conditions are used with open boundary problem, e.g. antenna
problems, to truncate the infinite space. Many versions o f the Absorbing Boundary
Condition can be found in the literature [l][33-40]. The theory discussed in this thesis is
applicable to both closed and open problems. However, in this research, only closed
structures are investigated, so no absorbing boundary condition is implemented.
2.3.4 Port Boundary Conditions
Consider an N-port microwave device (Figure 2.2), where each port may or may not
support a propagating mode. Let the transverse part o f the incident wave o f mode m at
port i be:
10
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 2.2: An N-port microwave device.
me = a mem
m m
^inc
(2?l0a)
H me
(2-lOb)
m m
where a? is the wave amplitude o f the incident mode, and e?,and h OTrepresent the
transverse electric and magnetic field distribution across the port for mode m,
respectively. em and h m are related as follows:
11
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where n is a unit normal vector outward from the device, i.e., opposite to the direction of
propagation (Figure 2.3). Z mis the wave impedance:
J 'm
z
(2-12)
m
where y mis the propagation constant o f mode m. e mand h mare normalized such that
when
is equal to 1, mode m will be then carrying a unit power, i.e. e mand h msatisfy
the relationship:
fem
I
nt x h?,
nt n d S = -1
(2.13a)
s,
Also, from mode orthogonality,
j e mxh? n dS = 0
s,
m *n
(2.13b)
The outward wave o f mode m at port i can also be written in terms o f e mand h m:
E
me^ m
^ out = bm
H out,
**
where
= - b ╗(0h
╗*"?╗*╗
(2.14a) '
v
(2.14b)
'
7
is the wave amplitude o f the outward mode m.
Hence, the total E- and H-fields due to mode m are:
E = Vm*e?nt
(2.15a)
H = I^m h nt
(2.15b)
Where
12
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 2.3: n is a unit vector outward from the device
Vо
Vо
m = aо
m +m6о
(2.16a)
Iо
m
(2.16b)
'0
m -A m
and Iо are called the generalized voltage and current of mode m at port i,
respectively. From mode orthogonality, the voltage can be calculated from the total Efield using the following projection [41]:
13
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.17)
V * (E )= -/E x h ..n d S
S,
where St is the cross section area o f port /.
Let port i be excited by unit pow er mode m only, i.e.
is equal to 1, with all other
ports perfectly matched, i.e. no modes are coming into the device through these ports. Let
the electric field solution within the device under these conditions be denoted
E о. Then,
<
the port boundary condition for port j can be shown to be equal to:
H ? = 2 S , h ? - X v f ( E ; ') h
(2.18)
n
where V f (eо ) is the voltage o f mode n at port j, due to exciting port i with unit power
mode m, and 5,y is the Kronecker delta.
The wave equation (2.5) and the above boundary conditions fully represent the field in
free space and inside any microwave-guided device. Solving them together analytically or
numerically gives the field distribution (E- or H-field).
2 .4 W e a k F o rm o f th e P ro b le m
Given the boundary value problem defined by the differential equation (2.5), and
subject to the boundary conditions described by formulae (2.7-8, and-18), the solution
E ^ to the boundary value problem m ay be obtained by solving:
л( E л 1w ) * 2 P л ( w )
(2.19a)
for all weight functions w, where a is the bilinear form
(2.19b)
+ Z г * ? ( g )* 7 (f)
V I
14
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where
tj0
is the intrinsic impedance o f the free space and it is approximately equal to 377
a.
Formula (2.19a) is called the weak form or weighted-residual equation. The same
result can be obtained if a variational formulation is used.
2 .5 F in ite E le m e n t D isc retiza tio n
In the finite element method, the region o f concern, is discretized into a set o f small
subregions. In 2D, the subregions could be triangles, rectangles, or 2D curvilinear shaped
elements, as shown in Figure 2.4. In 3D, the subregions could be tetrahedral, hexahedral
or 3D curvilinear shaped elements, Figure 2.5.
In the 3D finite element method, tetrahedrals are the most commonly used elements,
due to their capability to approximate any complicated 3D structure and the possibility o f
automatic generation o f tetrahedral meshes. In the 3D solver developed for this thesis,
tetrahedral elements were used.
By discretizing a 3D problem into a set o f M tetrahedra, the volume integral o f
formula (2.19b) can then be decomposed into M integrals:
6V r
л ( f . g ) - t - r 2- KVx g-A*;' -V x f - k 1g.er .f}dO.
I .J krio a
+ Z 2 > / < g) r / ( f )
*
1
where Vm is the volume o f element m, and
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.20)
C a)
CbJ
C c)
Figure 2.4: 2D Elements, (a) Triangle, (b) Rectangle, and (c) Curvilinear Element
(2.21)
C,r is the r th simplex coordinate o f a 3D tetrahedral element.
The field inside any element can be written as [42]:
16
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(b)
(c )
Figure 2.5: 3D Finite Elements, (a) Tetrahedral, (b) Hexahedral, and (c) 3D Curvilinear Element
E ( x , y , z ) = г г rN r(x ,y ,z)
(2.22)
N r is the 3D basis function associated with the r th unknown. The unknowns in a 3D
FEM nodal-solver are the components o f the E-field at each node, while in an edge-solver
they are the E-field along each edge in the mesh. In this thesis an edge solver was
developed.
17
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Let g in formula (2.20) be replaced by formula (2.22) for all elements, and replace f
with all possible Nu, i.e., u= 1.. U, where U is the total number o f unknowns. This
substitution generates a set o f U equations in the form:
[at]{ o
= { * :░ )
<2-33)
where {E} is a vector o f E-field values along all the edges, and
6V
At
r
j{ V x N r ./t;lV x N , - f c X - f i r N j d n
(2.г4)
m=l J t t l o a
+
Z E K*(IW
1 I
( N .>
<2-25)
The [K\ matrix is generated from the wave equation and all the boundary condition?
applicable to the problem o f concern. However, the {R^░} vector is constructed from port
boundary conditions only.
2 .6 S c a tte rin g M atrix
The performance o f any microwave device can be described by a scattering matrix [.S].
[5] relates the incident modes at all ports to the reflected ones:
{*} =[$]{<*}
(2-26)
where {a} and {b} are column vectors o f incident and reflected amplitudes,
respectively, o f all the propagating modes at all the ports. The diagonal terms in [5] are
the reflection coefficients, while the off-diagonal terms are the transmission coefficients.
Assuming again that the unit power mode m is incident at port /, while the other ports
were matched, then according to formula (2.26), the voltage o f an outward wave o f mode
18
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
n at port j is equal to the scattering parameter
, i.e.
V f(E 2 )= S r
(2.27)
The voltage o f the outward wave o f mode m at port i is equal to
vm
*(Eо)
2m? + i
\ nt / = Sm
(2-28)
The general formula for the scattering parameters is, then :
s r = V о (E j)-4 r
(2.29)
This formula is used to obtain [5] from the computed solution Eо .
2.7 P ro p a g a tio n M o d e s
The field distribution o f a mode o f an arbitrarily shaped waveguide can be obtained
using 2D finite element eigen-solver. Given the geometrical parameters and media
characteristics, the solver computes the field distribution and the propagation constants o f
many propagating and evanescent modes. These can then be used by a 3D finite element
solver to impose port conditions and to calculate the scattering parameters.
However, for many uniform guided structures there are analytical formulae for both
the field distribution and the propagation constants [43][44]:
2.6.1 R ectangular W aveguide
Figure 2.6a shows a rectangular waveguide with its port placed on the x-y plane.
Assume that a TE or TM wave is propagating inside the waveguide in the positive z
direction. Then, the transverse components o f the mnth TE mode are:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Cb)
Figure 2.6: (a) Rectangular waveguide (b) Circular waveguide
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
nn
e = - Z lm'n2
X
cos
(2.30a)
( mn \ .
x sm
\ a )
tttft
4 a tk L n
b
4 l + S0m4l + S0
(2.30b)
. f mn
(n n \
sin
x cos ? y
mn
v a
U
J
1/2
e.. = Z mn
4 abkc.mn ░
4 {jrS0m 4 l + Si
and the transverse components o f the mnth TM mode are:
___ 2
e = Z -1/2 Yn
jm
=
f ib k ,:
mn
( mn
^ . ( nn \
cos
x sin ? v
a
I л J lb yJ
(2.3 lb)
2
nn . f m n ^
f nn
,-------sin l------ xx | cos
< ? y
e v = Z~l
mnn' 2 ?
j m [^ JV ^ b k f Z b
{ a
C *M t
*
^ b
(2.31a)
'
where
Ym n
=
/c?c,mn =.
J \lk l
m n\
c r )
(2.32)
k c.n
(2.33)
( u ttY
+W
2.6.2 Circular Waveguide
Figure 2.6b shows another example o f a regular shape waveguide, that is the circular
waveguide. Again we assume that the port is placed on the r-<|) plane, and the wave is
propagating in the positive z direction.
Then, the transverse components o f the mnih TE mode are:
1/2
er = Z mn
j
2
+
n
)
J m( p ? r /a )
r 4 Pmn - m ' J * (Pmn ) ,
f-sin(mp)]
[cos(/np) J
21
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.34a)
e =Z
U2
^ mn
I
2
P mn
+
a
cos(/i<p)]
{
J m ( P mn r / a )
(2.34b)
- r n 2J m( p mn)y sin (m p)J
and the transverse components o f the mnthTM mode are:
^ - 1 / 2 y mn
_
r
^ mn
ip
P mn
I
2^
a
J n {pmnr/a) Vcos(/z(p)]
+ S 0m ) P mn
J m
(2.35a)
(P mn ) J l Sin(/I<p) J
(2.35b)
= Z mn
m1' 2
X ( l + to m )
P m n J m ( P mn
) J |c0s(/l<p)
where a is the radius o f the waveguide, p mn and p mn are the n-th roots o f the /w-th Bessel
function and its derivative, respectively, and
Ymn
= j \ j k2
(2.36)
~ k l?
IP m n 1 a
T E
IP m n 1 a
T M
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.37)
CHAPTER 3
The Gradient of Scattering Param eters
3.1 In tro d u c tio n
The scattering parameters of microwave devices change as their geometries are
modified. To improve the performance o f a device, one needs to modify some
geometrical parameters in a manner that leads to an optimum. To modify the geometry
properly, gradient information is very useful, i.e. information about the rate o f change o f
the scattering parameters with respect to a set o f geometrical parameters. The goal o f this
chapter is to evaluate these gradients, using the finite element method.
3.2 T h e A d jo in t M eth o d
From formulae (2.29) and (2.17), the scattering parameter Sf? can be expressed in a
general form as [2] :
(3.1)
where { E
} is the computed electric field within the device when port i is excited with
mode m. {Eг }} is itself a function o f g. Differentiating (3.1) with respect to g:
(3.2)
where
esU
i)
nmK
j.
is the column vector,
(7X0
23
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
U is the total number o f unknowns. Differentiating both sides o f formula (2.23):
d{E?}
dg
d[K]
dg
0) _ d{R ?}
dg
(3-3)
[/CJ ? ------ + ? ? \ L m ) ------ ------
Substituting (3.3) in (3.2) we obtain:
d s 'jr _ 3 s a r . <
:
*
?
a ,
Formula (3.4) involves the calculation o f the inverse o f the symmetric [AT] matrix. To
avoid such an operation, we introduce the adjoint variable vector {A.}, which is defined
as:
f ,
. dSm n
(x )= [ /:] '
3{гл(
(3.5)
Substituting in (3.4),
dg
= a s% r_ + {A(r a
dg
dg
(
})
<3-6>
Formula (3.5) can be rewritten as
[* ]W =
a s!? *
0-7)
3 ( 0
Formula (3.6) indicates that the adjoint method requires two 3D solutions in order to
calculate the derivative. First, formula (2.23) is solved to obtain {e (^ ]}, and then formula
(3.7) is solved to calculate the adjoint variable vector. However, once the adjoint solution
is found, derivatives with respect to any number o f different g ?s can be found.
Assume that g is an internal dimension which does not effect the location and the
shape o f any port in the device o f concern. Then the derivative o f
} with respect to g
can be dropped. Also since the scattering parameters are evaluated on the ports only, the
24
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
first partial derivative in formula (3.6) is equal to zero. Consequently, formula (3.6) can
be simplified to
=
(3'8)
dg
dg
Notice that formula (3.8) still requires two solutions to calculate the derivative.
3.3 T h e N a tu re o f th e A djoint V ariab le
In Section (2.6) we explained how
could be extracted from the computed
electric field using [4]:
S rtm
Jam = Vо
ft \(e 2
nt /) -
(3.9)
nnt
Differentiating formula (3.9) with respect to one o f the unknowns, say the field value
on the лth edge:
a g
3K u
_ d V f( E * )
(3-10)
d& L
Substituting (2.17) into the right hand side o f (3.10),
dSm
r)Ff0
m.u
d e
^ p(0
J m
m.u S ,
(3.11)
"
J
The field inside any element can be described using formula (2.22). Substituting (2.22)
into (3.11) yields,
dSm>
f
- ^ - = - / N ? x h ? .n r f S y
SE_
(3.12)
But from (2.17) the right hand side is simply ^nu )(N u) which, from (2.25), is equal to
?
If
2 n.u
╗
i
p-
25
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
So now the adjoint problem in formula (3.7) is simply equal to,
(3l4)
Comparing this to (2.23) we see that
{ ;.} = ! { г ? }
(3 1 5 )
i.e. it is in fact half the value o f the field generated from exciting port j with mode n. This
useful conclusion can be used now in formula (3.8) to obtain,
g
"
dg
(3l6)
2
"
dg
Here we are assuming that g is an internal parameter that does not affect the location
or the shape o f any port.
For an N-port model, to obtain the entire scattering matrix requires N field solutions,
given that one dominant mode is supported at each port. With the usual adjoint formula
(3.8), to getthe gradient requires an additional N adjointsolutions. However,
the new
formula (3.16)gives the gradients at no additional cost.Moreover, if just the return loss
o f the dominant mode o f one port is to optimized, then only one solution is required to get
Si i and its gradient.
3 .4 E v a lu a tin g th e D eriv ativ e o f [K]
In order to evaluate (3.16), we need the derivative o f [K\ with respect to g. Using the
chain ru le :
d\_K] = y y d[K ]drl
dG
V t dr I dG
(3.17)
where the first summation is made over all vertices which move when parameter g is
26
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
changed, and the second summation is over the three coordinates x, y, and z o f each
vertex, rj is the coordinate I o f vertex k.
From formula (2.24), each entry o f
is equal to:
A
dK n
drl
(3-18)
d (6V)
. . i.
/ ( v x N ^ . ^ - 'V x N . - * X N r .N ? )d n
A A n on
-
A'
+
1 J ( v x N r .A i;'V xN u - * .2er N r j O л л
A n 0n
(6F)
JI
jkoHo n I
x N. + V x N r.ju;
A
5(V x N b)
ar/
</Q
M
A'
17
, n ,m .
,
,
.
, .
dV a V x N
5Nr
Formula (3.19) involves three important derivative terms; ? - , ------ ?<-, ? f- . It can
dr[
dr[
drl
be shown that,
(3.19)
dV = v d гk
drlk
dr?
Further, using
(3.20)
dr
d r1
we get
= 2V gt x
drl
aw .
U r'
f
(3.21)
Sa
?
- ^ r ^ C ,b
3 r#
S*
ac╗
where N r and N uare,
27
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.22)
Nr - c .v c ,- e * v c .
(3.23)
(3.24)
Introducing the following definitions:
(3.25)
C a, =
dr1
(3.26)
d+
(3-27)
= 5 > ? ,c ?
/= I
(3.28)
e * * = 4 (V f. x V f 4) .( V ^ x V f ? ) = 4
- <*?</?)
formula (3.18) can now be rewritten in simple form as,
SAT
dr/
I
6V
A rja
(3.29)
'abed
?
г r ( ^ o c C^bd
I ad ^ b e
^b e ^ ad + ^ b d ^ ac )
6V
+
r*r
( C bl e kacd
C al e kbcd + C dl e kcab
C cl e kdab)
?6 V k a er(cal Ibc d u ?ca[ I m d kc ?cbl Iac d kd + cbl I ad d kc
+ Cct Ida d u, ?ccl I db d ka ?cdl I ca dkb + cdl Icb d^ )]
Using (3.17) and (3.29), it is straightforward and inexpensive to obtain the derivatives
o f [K\ for any geometrical parameter g.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
Optimization of Microwave Devices
4.1 In tro d u c tio n
The main purpose o f optimizing microwave devices is to meet design objectives such
as return loss, insertion loss, and coupling. There are many techniques for optimizing
microwave devices, but, since we have gradient information available at very little cost,
only optimization techniques that make use o f gradient information are o f interest.
In section 4.2 we present a general formula for a cost functions used in microwave
design, and the corresponding gradient formula. In section 4.3, we describe the conditions
for a point to be a local minimum, then in section 4.4 we discuss, in general, the theory o f
the gradient methods. Different kinds o f optimization techniques are discussed in section
4.5, while in section 4.6 we present the Quasi-Newton methods, which are the methods
used in this research. In section 4.7 the Line Search method is presented in details.
4 .2 C o s t F u n c tio n f o r M icrow ave D e v ic e s
A general cost function o f an N-port microwave device is [6]:
(4.1)
where a ',
and P ~ are weight functions, S - { f t ) are the desired scattering parameters,
and the frequency range has been discretized into Ah frequencies. The a ~ and p~
weights control the emphasis on the magnitude and phase repectively o f the scattering
parameters. From this general formula, one can also deduce the formula o f the gradient o f
29
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the cost function with respect to any geometrical parameter, say g. The derivative o f the
cost function is a function o f the scattering parameters and the derivative o f the amplitude
and/or the derivative o f the phase o f the scattering parameters. The derivative o f the
amplitude is equal to:
r f |W < ) |
dg
R efc, ( / , )}R e(8,?; ( / ' *1 +
dg
( / , )}lro(35╗ ( / ' }
dg
(4.2)
? W /)I
while the derivative o f the phase is:
(4.3)
<US,
dg
In terms o f these, the general formula for the derivative o f the cost function is:
f-? ? ?
(J , )
dg
f c ,( .f, )| - | s , ( / , )|)+
d Z S jX fi) (
*
(4'4>
\
4 .3 C o n d itio n s fo r a L o cal M inim um
In general, the cost function C can be written as a function o f a vector {g} o f
geometrical parameters. Let {g}={g*} be a local minimum, and C be continuously
differentiable at {g*}. Then the following condition must be satisfied [45-48]:
V C \g* = 0
(4.5)
i.e. the slope o f C is zero at {g*}. However, satisfying (4.5) does not necessarily
mean that {g*} is a local minimum; it could be a local maximum. To obtain minima only,
an additional criterion is needed. The new criterion is derived from the second derivative
30
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
o f the cost function C.
Let C be twice continuously differentiable at {g*}. Then, one can expand C in the
neighborhood o f {g*} to obtain the formula:
c({g-}+ {a})= cfe*})+
where
)
(46>
is a very small perturbation vector, and the [H] matrix is the Hessian o f C:
(4.7)
>J
Sgi8gj
g
*
Higher order terms are neglected in formula (4.6). For {g*} to be minimum, the left-hand
side o f (4.6) should be greater than the cost function at {g*}. Consequently, the additional
required criterion for minima is:
(4.8)
{5}r [tf]{5 }> 0
V {S }* 0
Formula (4.8) can be satisfied if [//] is a positive definite Hessian matrix.
4 .4 M ultivariate G ra d ie n t M e th o d s o f M in im izatio n
In all gradient methods, a vector {g} o f geometrical parameters is iteratively updated
until the objective cost, C, becomes lower than a prespecified threshold. The A+1st update
is o f the following form:
(4.9)
where {p(kJ} is certain search direction, and
is a step size in the direction o f {p(k)}-
Gradient methods use the slope (gradient) o f C at {g?k)} to construct a suitable search
direction. Some gradient methods do not require derivative information, like Powell?s
31
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
method. However, we are only interested in those techniques that make use o f the
derivative information such as, steepest descent, Newton, quasi-Newton, and the
quickprop method.
The gradient methods differ from each other in the w ay the vector {p(k>} is updated.
The value v(k) can be constant, or it can be calculated using single-variable optimization
(line search) technique.
In the steepest descent method, developed by Cauchy, the vector {p(k>} is given by:
{/>"}=-vcrfe};
(4I0)
i.e. the direction o f search is along the line o f maximum rate o f change o f C. The steepest
descent method has a linear convergence rate for a quadratic function. This is due to the
fact that the second derivative is not used. Thus, the method is not usually recommended
for general applications.
On the other hand, the Newton Raphson method uses the second derivative
information to update the search direction vector {p(k>}- The formula used is:
{p'l))= [H m y 'c u g } )
<4 -n>
where [//] is the Hessian matrix o f equation (4.7) which includes the second derivative
information.
The main disadvantage o f the Newton-Raphson method is the cost o f calculating [//],
which rises dramatically as the number o f variables increases. In addition, singular
Hessians, indefinite matrices, and overshooting could be possible sources o f problems to
this technique.
For conjugate gradient methods, the search direction {p?}, satisfies the conjugacy
condition [48]:
32
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
m ^n
(4.12)
One can show that {p(k)} vector can be generated iteratively without constructing the
Hessian matrix [//], using formula:
(4.13)
{p(m,}= -W C (k>+
Conjugate gradient methods are highly efficient when the number o f variables is large,
and no information is available on the second derivatives, or it is very difficult to obtain
them. An improved convergence rate can be accomplished through exact line searches to
obtain the optimum step size v(k) in the direction o f {p(k)}.
In the quickprop method [29], the cost function is approximated by a parabola which
has arms open upwards. The assumption made is that the slope o f the cost with respect to
a variable is not affected by the change in the values o f all other variables in the model o f
concern. Then, one can easily show that:
(4.14)
Rekanos showed that the QP technique is 3-4 times faster than the steepest descent
method. No comparisons with Newton or quasi-Newton techniques are available yet.
4 .5 Q u a si-N e w to n M e th o d s
To avoid the calculations o f the exact Hessian matrix in the Newton-Raphson method,
the quasi-Newton methods use an iterative scheme to approximate it. The gradient and the
33
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
cost function are used in the calculations, which take place at each iteration.
The search direction {p(k)} is found by:
(4.15)
{pw } = -[A fa ) ]{vCa ) }
and the [M\ matrix is updated using the following formulae [45][47-51][52]:
(4.16a)
(4.16b)
W
J
W
w
\ r )
(4.16c)
2(2(5+ l - 3 i j )
Where
{5(i)}r { v c (i_I)} ?
C (t) - C (*~l)
{/*>}= {vc(t)}- {vc(?_1)}
and,
Formula (4.16) is known as the Broyden-Fletcher-Goldfrab-Shanno (BFGS) method.
There is another quasi-Newton method called the Davidon-Fletcher-Powell (DFP)
method. Powell [53] showed that BFGS performs better than DFP in minimizing
quadratic functions o f two variables, which suggests that BFGS should be used also for
34
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
general non-linear functions.
To improve the overall performance o f the algorithm, a good initial value for [M] is
necessary [54]. A convenient choice is to set the initial value to v(0)[/o], where [70] is the
identity matrix. v(0) is the initial step size.
When using a line search, the value o f the step size J k) converges to 1.
The algorithm o f the quasi-Newton module has the following steps for minimization:
1. Begin with [M]=[/a].
2. Using the initial geometry, calculate the cost and its derivative.
3. Call the line search module, and perform an approximate line search along
{/>(0)}= -{V C (0)}. The line search module returns: the new step v(0), the Cost and its
derivative at v(0)
4. Set [A/0)]=
v (0)[ / o].
5. Evaluate formulae (4.16).
6. Evaluate formula (4.15).
7. Update the geometrical parameters using formula (4.9).
8. Call Line Search which returns under three conditions:
8.1 If the cost or its derivative are better than the target values, return.
8.2 If the geometrical parameters were not changed, return.
8.3 Otherwise proceed to step 9.
9. Redo 5-8.
4 .6 L ine S e a rc h
To improve the convergence rate o f any gradient method, a line search is used to
update the value o f the step size v/*; at each iteration. The goal o f the line search is to find
35
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the v╗ which makes C({g(k)}+ \fk> {p(k)}) locally minimal, i.e.
(4.17)
Using the chain rule, this becomes
{p(i)}r { v c (*+l,} = 0
Formula (4.18) represents the condition for the termination o f an exact search.
(4.18)
The line search module is composed o f two consecutive steps: establishing bounds for the
search, and allocating a minimum along the {p(k)} vector.
4.6.1 Establishing Bounds for the Search
In order to find a minimum, the line search module needs to identify an interval
bounding at least one local minimum. Making use o f the derivative information,
Davidon?s polynomial extrapolation method is used to allocate the upper or lower bound
o f this interval.
Let a be a given lower or upper bound for v, depending on the sign o f C '. Then an
estimate for the minimum is:
(4.19)
b =a +2
where C? is the cost at a, and C, is the target cost (a user-defined value). C ' is the slope
o f the cost function with respect to v at point a.
Formula (4.19) can be applied iteratively : the program starts with a given a, Ca, and
C ', then estimates the upper bound o f the interval i f C ' is negative, or the lower bound
o f the interval if C ' is positive. If the sign o f the slope o f the cost function at the new
position b is the same as at a, the program changes the value o f a to b, and seeks another
value for b. The process stops when the slope changes its sign.
36
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The algorithm o f the interval step includes the following steps:
1. Start with a=0.
2. I f the slope o f the cost function with respect to v is negative then,
2.1. Evaluate b using formula (4.19).
2.2. Evaluate formula (4.9) using v=b.
2.3. Solve the problem and calculate the cost function and its derivative.
2.4. I f the cost function is below the target cost function Ct, return.
2.5. I f the derivative o f the slope is negative, then,
Set a~b, Ca=Cb, and C ' = C 'b , then go to 2.1
2.7 Else return.
3. If the slope o f the cost is positive then,
3.1. Set b=a, Cb=Ca., C'b = C'a
3.2. Evaluate a using formula (4.19), where a is replaced w ith b and b with a.
3.3. Evaluate formula (4.9) using v=a.
3.4. Solve the problem and calculate the cost function and its derivative.
3.5. If the cost is below the target cost function C,, return.
3.6. If the derivative o f the slope is positive then,
Set b=a, Cb=Ca, and C'b = C ', then go to 3.1
3.7. Else return.
4.6.2 Locating a Minimum
Once an interval bounding a minimum has been determined, polynomial interpolation
methods can be used to approximate the cost function by a polynomial, from which a
37
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
possible minimum is calculated. In this thesis, cubic interpolation is used. The method
exhibits quadratic convergence.
It can be shown that if the cost function and its derivatives are known at the bounds o f
the interval ([a,b]), the minimum o f a cubic polynomial fitting these values is given by:
(4.20a)
(4.20b)
(4.20c)
w
Formulae (4.20) are used in an iterative process to obtain the real m inim um o f the cost
function within the given interval. For each new ?cubic? m inim um , the cost function and
its derivative are evaluated. I f the derivative o f the cost function is not small enough, then
the calculated point is not a real minimum. However the new information is not wasted,
but used to narrow the interval from one o f the two sides (lower or upper) depending on
the slope o f the cost function at v. This process is repeated for the new interval, until the
slope o f the cost is below certain threshold, or the interval is less than a given intervalthreshold.
The algorithm for the allocation o f a minimum is:
1. Calculate v using formulae (4.20).
2. Evaluate formula (4.9).
3. If the differences between the new and old geometrical parameters are below the
geometrical tolerance, return
4. Adjust the geometry and solve the new model.
5. Calculate the new cost function and its derivative.
38
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
6. If the new cost function or its derivative are better than the target cost Ct return.
7. If the derivative is negative then set a=v, Cfl=Cv, and C ' = C ', and go to 1.
8. Else set b ?v, Cb~Cv, and C ?b = C ', and go to 1.
39
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 5
Computer Program Implementation
5.1 In tro d u c tio n
Using the theory discussed in the previous chapters, a Fortran-99 program was
developed.
The program consists o f two main parts: the optimizer and the 3D finite element
package. The optim izer is made o f two modules: the Quasi-Newton and the Line Search
modules. The 3D finite element solver is made o f four modules: the m esh generator, the
node-to-edge converter, the 3D finite element solver, and the remodeller. The program
interface is via input and output files. The program needs four input files which are :
Geo.dat, 2D.dat, Input.dat, Ports.dat, Mats.dat, and Param.dat. It generates three
temperory files: Tets.dat, Edges.dat, and Results.dat, and the output o f the program is one
file: Tracing.log.
In Section 5.2 we describe the six modules o f the program, then in section 5.3 we
present the different data files used during the run.
5.2 P ro g ra m S tr u c tu r e s
The flow chart shown in Figure 5.1 demonstrates the interconnection between the six
main modules o f the program.
5.2.1 Quasi-Newton Optimize
The method implemented here uses the BFGS technique described in the previous
chapter. The efficiency o f this optim izer depends on how fast and accuratly the cost
40
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
function and the gradient information are obtained.
The Quasi-Newton method is responsible for generating new geometric param eters.
The Quasi-Newton calculates the search direction p, and calls the Line Search module to
obtain the minimum along p.
Start ^
E valuate Cost Function
Call
Initialization
\
?
f
S
Remodeller
j
Return
Mesh
Generator
/?
V.
( a d
'
y
/
Quasi-Ne^on
Module
C all
Nodal-to-Edge
Converter j
Return
Line Search
Module
Call
FEM Solver
Return
Figures. 1 Representative flowchart o f the relation between the six main modules : quasi-Newton, line
search, remodeller, mesh generator, nodal-to-edge converter, and the 3D FEM solver.
5.2.2 Line Search
The Line Search module is composed o f two steps: in the first step Davidson?s quadratic
polynomial extrapolation is used to construct an interval bounding at least one local
41
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
minimum, and in the second one, the cubic polynomial interpolation is implemented to
minimize the cost function along a search direction p.
The Line Search may not perform the two steps all the time, as discussed in the previous
chapter. The user may also choose to discard the use o f the Line Search option.
5.2.3 Mesh Generator
The mesh generator was written by Marc P. Choufani (Master Project) [54]. The user
supplies two data files: a module file and a design parameter file. The program generates
a data file which contains the (x,y,z) coordinates o f all the points, elements material labels
and nodes, and the boundary conditions o f each o f the four faces o f all the elements. The
program uses the Delauny triangulation method to construct the 2D mesh, then uses the
extrusion method to make the final 3D-nodal mesh.
5.2.4 Node-to-Edge Converter
The 3D mesher generates a 3D mesh in which each node is assigned a global number.
However, the 3D finite element solver is a 3D-edge solver. The Node-to-Edge converter
module assigns a global number to each edge in the 3D mesh. The output file, Edges.dat,
contains information on all edges and their associated nodes (two nodes per edge), and all
elements and their associated edges (6 edges per element).
Some elements have the number zero assigned to some o f their edges, this indicates
that theses edges are located on perfectly conducted surfaces, i.e., homogeneous Dirichlet
Boundaries.
5.2.5 3D Finite Element Solver
This module performs four different operations: constructing the [K] and [R] matrices,
solving the linear set o f equations, calculating the scattering parameters, and finally
calculating the derivatives. Formula (3.1a ) is used to construct the [K] and [R] matrices.
42
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The no-look-ahead QMR algorithm is used to solve the set o f equations developed in
chapter 2. The preconditioning used in the QM R solver is the diagonal o f the [K] matrix.
Formula (3.12) is used to calculate the derivatives o f the scattering parameters. These
derivatives are then used to calculate the gradient o f the cost function. Appendix .1
contains a copy of the QMR algorithm.
5.2.6 Remodeller
The remodeller modifies the geometry in accordance with the latest values o f the
geometric parameters. It modifies two data files: the model file, Geo.dat, which contains
the geometrical information and the derivative file, Param.dat, which contains
information on the perturbed surfaces. The mesher uses the model file, while the 3Dfinite element solver uses the perturbed file.
5 .3 D a ta F ile s
The program communicates through a set o f data files. Each procedure o f the program
reads data from specific data files, and outputs the results to another set o f data files.
5.3.1 Geo.dat
Geo.dat is the model file and it contains information about the geometrical dimensions,
the materials, and boundary conditions o f the device under investigation. The Remodeller
modifies Geo.dat prior to calling the FEM solver.
5.3.2 2D.dat
2D.dat is the mesh control file. 2D.dat contains some controlling parameters to be used
by the 3D Mesher. One important parameter in 2D.dat is the largest edge length allowed
in the system.
43
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
5 3 3 Inputdat
The main program reads Input.dat once at the beginning o f the optimization process.
Input.dat contains information about: the number o f variables and their initials, maximum
and minimum dimensions, frequency points, QM R solver tolerances, Quasi and Line
Search tolerances, and the number o f geometrical precessions tolerated.
5 3 .4 Ports.dat
Ports.dat contains information about the location, the shape, and the exited modes at
all the ports in the system. This file cannot be changed by the Remodeller, since the
theory assumes that the dimensions, locations, and excited modes o f all ports are fixed.
5.3.5 Mats.dat
Mats.dat contains a library o f many materials w ith their permitivity and permeability
parameters. Each material is given a label, and the same label should be used in the model
file.
5.3.6 Param.dat
G2mesh.dat is the derivative file. To calculate the derivatives o f the scattering matrix
with respect to any geometrical parameter, the program needs to know which surfaces are
related to these geometrical parameters. G2mesh.dat contains information about all
variables, the surfaces related to them, and the direction o f perturbation of each o f these
surfaces. The program supports only rectangular surfaces at present.
5.3.7 Tracing.log:
Tracing.log contains information about all intermediate results obtained during the
optimization. For each trial, the program saves the new values o f the variables, the cost
function and the derivative o f the cost.
44
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Tets.dat is the output o f the mesher. Tets.dat contains information on the cartesian
coordinates of all the points, the nodes o f each element identified by their global number,
the material label o f each element, and the boundary conditions on each surface o f each
element.
5.3.9 Edges.dat
Edges.dat is the output o f the node-to-edge converter. Edges.dat contains information
on the nodes o f each edge and the edges o f each element identified by their global
numbers.
5.3.10 Results.dat
The 3D finite element solver solves for the field value along each edge in the model.
The last solution is saved by the system in Results.dat.
45
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 6
Results
6.1 In tro d u ctio n
The twofold purpose o f this chapter is (1) to verify the theory o f the gradient
developed in Chapter 3, and (2) to demonstrate the optimization technique described in
Chapter 4.
Three examples are presented. The first example is a short-circuited rectangular
waveguide, the second is an E-plane miter bend, and the third is a single stage
waveguide transformer.
In the first example, the gradient is verified by comparing numerical results with
analytical ones, hi the second and third examples, the gradients are verified by numerical
perturbation.
In the first two examples, the devices have one variable and are optimized at only one
frequency, while in the transformer case, three variables are optimized simultaneously
and at eight different frequencies.
6 .2 S h o rt C ircu ited R e c ta n g u la r W a v e g u id e
Consider the uniform rectangular waveguide in Figure 6.1 with inner dimensions: a
(in the x-direction) and b (in the y-direction). The waveguide is short-circuited at a
distance L from the port. The waveguide is operating in the fundamental mode TEio. The
dimensions are as follows:
46
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
b/X =0.1
Xg/X = 1.515
L is the design parameter. X is the wavelength o f the operating frequency.
Short?Circu ited
Side
b
Port
Figure 6.1: Rectangular waveguide short-circuited at a distance L from the port
The design criterion is to change L until the desired phase o f Su is obtained. The cost
function is defined as:
C=
z s ^ fj-zs,^
(6 . 1)
where the cost is evaluated at a single frequency point f a. The gradient o f the cost
47
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
function is:
(6.2)
The gradient o f the cost with respect to L is a function o f Si i and its gradient with
respect to L
6.2.1 Verifying the Gradient o f Z S?
Analytical formulae can be obtained for the gradient o f г u o f a rectangular
waveguide short-circuited at a distance L from the port:
S u = -cos(2kzL) ?jsvn(2kzL)
dL
= 2k,(sm(2kzL) ?j cos(2 k zL))
(6.3)
(6.4)
where kt is the phase constant o f the TEio mode. Substituting (6.3-4) in (6.2) to obtain
the exact formula for the gradient o f C. The values extracted from FEM solution using
formula (3.16) were compared with the exact values. The results are displayed in Figure
6.2 for LI A =0.0667, ..., 0.667. Good agreement between the finite element solutions
and the analytical ones can be concluded. The error in the phase o f S u is also displayed
in Figure 6.2 to demonstrate the accuracy o f the 3D finite element solver. Maximum
percentage error in the phase occurred at the location where the phase is almost zero, as
expected.
48
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
8
6
<0
c
a
?o 4
2.
Error in the
Phase %
Error in the
Gradient %
л
10 2
a.
h
o
0
JO
*A
-a-*"A .
Phase of
S11 (radians)
з ?2
UJ
A
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
L /L a m b d a
Figure 6.2: Error in % of the phase o f Sn and its gradient
6.2.2 Optimization
The initial value for the optimization is chosen to be:
?^initial / A = 1 /3
The target phase is -0.7856 radians. The initial phase o f Si i differs from the target by
146%. The results o f the optimization are displayed in Figure 6.3, which shows the error
in the obtained phase after every 3D solution. Figure 6.3 shows how the program
converged from the first iteration to a solution that is very close to the final one. This
fast convergence is because the cost function is a quadratic function o f L I X as shown in
Figure 6.4. Figure 6.4 presents the cost function in the vicinity o f the optimum solution,
starting from the initial value. This result demonstrates both the accuracy o f the gradient
49
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
calculations and the powerfulness o f the optimization technique used (quasi-Newton +
line search). The final value o f /-initial / A is 0.471.
160
140
120
100
80
60
40
20
Solution#
Figure 6.3: Error in the phase o f S1{ versus solution number
50
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0.3
0.35
0.4
0.45
0.5
0.55
0.6
L/Lambda
Figure 6.4: Cost function in the vicinity o f the optimum solution
6.3 E -p la n e M iter B e n d W a v e g u id e
The second example is the miter-bend shown in Figure 6.5. The bend is in the Eplane, i.e. the 6-dimension. The device has two ports; port one is located on a y-z plane,
while port two is located on an x-z plane. The operating mode is again the fundamental
TEio mode.
The dimensions o f the rectangular waveguide are:
a/X = 2/3
b/X = 1/3
The guided wavelength is:
Xg/X = 1.515
51
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
X is the free space wavelength at the operating frequency, and Xg is the corresponding
guided wavelength. The design parameter g is shown in Figure 6.5, and the movable
boundary related to it is the inclined wall o f the miter-bend. Increasing or decreasing the
design parameter g causes the inclined wall to move diagonally (45 degrees) outward or
inward respectively.
g
b
Figure 6.5: E-plane Miter Bend.
52
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The cost function is defined as:
(6.5)
where the cost is evaluated at a single frequency point f 0. The gradient o f the cost
function is
(6 .6)
The gradient o f the cost function with respect to g is a function o f the return loss S u
and its gradient with respect to g. In our case the desired amplitude o f S u was chosen to
be 0.
6.3.1 Verifying the Gradient
To calculate the gradient, the program needs to know in which direction each mesh
node on the inclined wall is moving when the geometrical parameter g changes. For the
inclined wall o f the miter bend there are many ways of moving the mesh nodes. The
method chosen is to divide the surface o f the inclined wall into two parts: the mesh
nodes located in the lower part are moved in the y-direction, while the ones on the upper
part are moved in the x-direction, as demonstrated in Figure 6.6.
Given this information the program can calculate the gradient o f S\ \ with respect to
the design parameter g. Since there are no analytical formulae for this problem, we use
the finite difference technique to calculate the gradient. To do that, the design parameter
53
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
g is perturbed by ▒ 0.5 x 10-3:
^^Ig+O -SxlO "311 g-O.SxlO~3
^
A g ~ 1.0 xlO -3
To be able to compare properly both the finite difference and the finite element results,
the mesh nodes o f the inclined wall m ust be displaced in the direction that was adopted
when calculating the gradient. This means that the nodes o f the lower part o f the inclined
wall are displaced in the y-direction, while the ones on the upper part are displaced in
the x-direction.
Upper Section
Lower Section
To P o r t
Figure 6.6: Displacement of the mesh nodes located on the surface o f the inclined wall
54
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Using formula (4.2), we calculated the gradient o f the amplitude o f Su for both the
finite difference and the finite element methods. The results are presented in Figure 6.7.
Very good agreement between both solutions can be stated, g/b was varied from -0.8 to
0.8. The locations o f the bend for g/b=-Q.% andg/&=0.8 are shown in Figure 6.8.
0.02
0.01
-
0.8
-02
02
0.4
o -0.01
i
-
0.02
-0.03
-0.04
Figure 6.7: % Error in the gradient of |5;/| o f the miter-bend example
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.6
0.8
6.3.2 Optimization
The initial value o fg is:
Sinitial ^
~0.8
The results o f the optimization are displayed in Figure 6.9. The optimizer converged
after 6 solutions to within 0.01 o f the final solution, i.e. 1% of the b dimension. The
optimizer performed additional solutions before it hits the geometrical tolerance {g/b
geometrical tolerance = 0.0001). Better performance can be obtained by choosing a
better initial value, g/b=0.0. The optimizer needed only 5 solutions before reaching the
geometrical tolerance.
Cbl
C al
Figure 6.8: Miter-bend at two different positions: (a) g/b=-0.8, and (b) g/b=0.8
56
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.8
0.6
0.4
02
j
0
-lnitia!=-0.8 j
-lru'tial=0.0 |
f
-02
-0.4
-
0.6
-
0.8
-1
Solution#
Figure 6.9: Miter-bend position as a function o f the solution number
Figure 6.10 presents the cost function versus solution number. The results indicate
that to within +/-0.01 o f the final solution the cost function is around -40 dB, i.e. the
return loss is 40 dB. To get better performance a more accurate g/b is required. In this
problem, up to 4th digit o f precision was needed to reach beyond -6 0 dB. When the
initial value was -0.8, the optim izer needed extra six solutions to reach -6 0 dB, but for
the better initial value, i.e. 0.0, it used only three additional solutions.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Cost (dB)
0
lnitial=-0.8
lnitial=0.0
Solution #
Figure 6.10: Cost (dB) versus solution number
The single step E-plane miter-bend is a microwave component with a single resonant
frequency. The return loss over a range o f frequencies is as shown in Figure 6.11. In
Figure 6.11, the optimum design was used to test the return loss for A /a from 1.072 to
1.875. For such a device one does not need to optimize at more than one frequency. It is
sufficient to optimize the bend at one frequency, usually at the desired center o f the
frequency-band, to get the best optimum dimensions. The miter-bend was optimized at
A /a= l.5, and hence the curve is centered at 1.5 as expected. From the graph, for 40dB
performance, the designed miter-bend can be used with 24%-bandwidth systems.
58
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
1.2
1.6
1.4
1.8
Lam bda/a
Figure 6.11: Retum-Ioss o f the optimum miter-bend
6 .4 S in g le S ta g e W a v e g u id e Im p e d a n c e T r a n s fo rm e r
The final example is the uniform rectangular waveguide transformer shown in Figure
6.12. The transformation is in both a- and b- dimensions. The device has two ports, and
the operating mode is the fundamental TEio mode.
The dimensions o f the input and output ports are
Qinput
2.0 m
binput ~
0.4 m
a output= =2.4
m
00
o
m
boutput=
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The optimization is performed over 8 frequencies equally spaced from 95 MHZ to
105 MHz: 95.000, 92.143, 94.286,96.429, 98.572, 100.715, 102.858, 105.000 MHz.
Output Port
Input Port
Figure 6.12: Single stage rectangular transformer
There are three design parameters in the model: width, height, and length o f the middle
section. For the width and height, i.e. the a-dimension and 6-dimension, there are two
movable boundaries, while for the length there is only one movable surface.
The cost function is defined as:
c =i;(is?(/.)|-|i? if.tf
(6'8)
/= i
The gradient of the cost function is similar to (6.6), with g replaced by a, b, or /. The
target value at all frequencies was chosen to be equal to 0.
60
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
6.4.1 Verifying the Gradient
Like the miter-bend example, there are no analytical formulae for the gradient o f S u
with respect to a, b, o r / dimensions o f the middle stage. Consequently, the finite
difference technique is used again to calculate the gradient. To do that, the design
parameters were perturbed by ▒ 0.5 x 10-3:
t L o5лio-) - i L-o.s*io-?_
1.0 x lO '3
Aa
A Sn
Ab
A5?
AI
_
^
1 L+O -SxIO "3 ~ <^ l l | * _ 0 3 x l 0 - 3
(6.9a)
(6.9b)
1.0 XlO'3
_
'^ ? ? L + O J x I O '3 ~ ?^ I l |/ - 0 . 5 x i ( T 3
(6.9c)
1.0 XlO'3
In Figure 6.13 both the finite difference and the finite element results are compared
for one geometrical position at eight different frequencies. Notice that the difference
between FEM results and (6.9a) and (6.9c) were magnified 100 times. Good agreement
between both solutions can be observed.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0 .5
>8
w
з
a -2 .2 C 100) !
-0.5
b= 0.6
;
UJ
-1.5
90
95
100
105
F r e q u e n c y (M H z)
Figure 6.13: % Error in the Gradient o f |Sy/|
6.4.3 O ptim ization
The initial values o f a,
b,
and / o f the middle section are:
Q initial ~ 2 . 2 IT1
hnitiai
Ml
b in itia l
0 .6
I initial
1 .0 m
was chosen to be roughly a quarter wavelength as it should be for a quarter-
wave transformer [43]. The results o f the optimization are displayed in Figure 6.14,
which shows the cost as a function o f solution number. The results demonstrate how
choosing the proper initial values leads to only few solutions before reaching the
optimum solution. In our case here, 3 solutions were enough to reach a very good
62
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
solution.
The final dimensions are:
afmai = 2.228 m
bfmai = 0.5768 m
Ifrnal
= 1.0000 m
-25
-26
-27
-29
-30
-31
-32
1
3
2
4
!
5 !
Solution#
Figure 6.14: The cost (dB) as a function o f the solution number
Figure 6.15 presents the frequency response o f the transformer before and after the
optimization. A small improvement was observed at both ends o f the band, while a
major improvement was accomplished for the center frequencies. The final frequency
response is roughly symmetric about the center frequency. Mode matching results were
also presented to validate the results obtained from the finite element solver.
63
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
50
45
40
?Initial Response
-Final Response
-Mode Matching (30 inodes)
9 35
K 30
25
2 0 -I
95
96
97
98
99 100 101 102 103 104 105
Frequency (MHz)
Figure 6.15: Return Loss (dB) for the initial and the optimized geometry
6 .5 S u m m a ry
In this chapter we presented three examples: a short-circuited rectangular waveguide,
an E-plane miter bend, and a single stage waveguide transformer. The design criterion
for the first example was to change the length o f the waveguide until the desired phase
o f S\ i is obtained, and for the second and third examples the design criterion was to
change some geometrical parameters to optimize the return loss response for one or
group o f frequencies.
For the three examples presented here, the gradients o f the amplitude or phase o f S u
were verified versus analytical or finite difference results, and the agreement were
remarkable.
64
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The optimization performed on the three cases took between 4 and 12 3D finite
element solutions to converge. It was demonstrated that good initial values would lead to
fast convergence. However, even when the initial values were bad, as in the first case for
the miter-bend, the optimizer managed in 5 solutions to get close to the final solution.
Such results would not been obtained if the gradient information was not available, or
was not calculated accurately enough.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CH APTER 7
CONCLUSION
7.1 Summary
The adjoint method provides a way o f calculating the derivatives o f scattering
parameters with respect to geometric parameters. However, straightforward application o f
the method requires additional adjoint solutions. In this thesis, a new formula for the
derivatives was presented requiring no additional solutions.
A finite element program package was created for the design of microwave devices.
The package uses a quasi-Newton method, specifically the BFGS method, to do the
optimization and makes use o f the efficiently calculated derivatives. The line search
method was used also to calculate the step size, to improve the convergence rate o f the
quasi-Newton module.
The new derivative formula was verified using three different examples: a shortcircuited waveguide, a miter bend, and a 3D transformer. The first model has analytical
solutions for both its return loss and the gradient, while for the second and third examples
a finite difference approach was used to calculate the gradient. The values calculated
from the derivative formula closely matched the values obtained from analytical or finite
difference calculations.
The three models were then optimized under different criteria. The criterion in the
short-circuited waveguide model was a specific phase o f S u , while for the miter bend and
the transformer, the criterion was minimum amplitude o f S n. For the three models, the
66
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
optimizer needed only 2-5 solutions to get very close to the final solution, and only 4-12
iterations to reach the final solution.
7 .2 O riginal C o n trib u tio n s
The original contributions o f this work are two fold:
(a) An improved version o f the adjoint formula to calculate the gradient o f the
scattering parameters. The new formula requires half the number o f solutions needed by
the standard adjoint method.
(b) A simple and inexpensive formula for the calculation o f the derivatives o f [AT] with
respect to any geometrical parameter, for the widely used tetrahedral edge element.
7.3 S u g g e s tio n s fo r f u r th e r w o rk
(a) Allowing Ports to vary: In this dissertation, the ports are assumed to be
unchangeable: the shape, the location, and the excited modes do not vary with the
geometric parameters. One could look into removing some or all o f these conditions.
(b) Antenna Problems: The approach may be extended easily to antenna problems.
In such problems, one can assume that the ports and the absorbing boundary are
unchangeable. Such an assumption does not effect in any way the design o f the antenna,
since the feed port is usually known, and the absorbing boundaries are fictitious boundary
located at a distance from the antenna.
(c) Curvilinear Elements: Instead o f using tetrahedral elements, one could study the
effect o f using curvilinear elements. These elements have the advantages of
approximating more curved boundaries more accurately [55-57].
67
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(d) Software Improvement: In the QMR matrix equation solver used, the
preconditioning was simply the diagonal terms o f the [AT] matrix. Better preconditioning
[58] can be implemented to improve the convergence rate. Also replacing both the nodalmesher and the nodal-to-edge modules, with a fast edge-mesher could improve the speed
o f the whole package.
(e) Mesh Refinement + Optimization: One interesting subject to look into is how to
simultaneously, optimize the device and refine the mesh. Such a combination is very
promising, since the refinement o f the mesh is a time-consuming step in optimization.
(г) Discontinuous Boundaries: It was observed by many authors [59], that there is
an inherent possibility o f converging to an oscillating or saw-toothed boundary as an
optimal solution. This can be overcome by constraining the geometric parameters by
smooth functions, e.g. constraining the geometrical parameters to a Bezier curve [60].
68
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX 1
Matrix Assem bly for 3D Edge Elem ents
High frequency software generally uses edge elements rather than nodal elements,
because nodal solutions suffer from the effects o f non-physical spurious solutions [I], the
inconvenience in imposing boundary conditions at material interfaces, and the difficulty
in treating conducting and dielectric edges.
For a tetrahedral element, with nodes labeled locally from 1 to 4 as shown in Figure
A l.l, one can define six basic functions N -, /= 1,
6. Using these six basic functions,
the field inside a tetrahedral element can be described in the following form:
л
E (x, y , z ) = 2 г,N ,(x, y , z)
l╗l
The definition o f N,. is,
(A l.l)
(A1.2)
Nf = C aVг6 - г 6VC0
Table A l.l presents the values o f a and b for i =1 to 6.
Table A l.l Edge Definitions for a Tetrahedral Element
Edge Number
I
1
a
b
1
2
2
1
3
3
1
4
4
2
3
5
4
2
6
3
4
69
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Given these definitions and formulae (3.23-24), the first term on the right hand side o f
formula (2.24) can be rewritten in more compact form as [1]:
6K
jkr\a
^
/{ V x N ,./J; l V x N ? } d n =
AV
jk r\ (6V ) ^
(A1.3)
~ darCbr Xpmdbu ~ d aucbu)
+ (d a r K
- b a r d h r \ d a u b bu - b a u d b u )
( b a rc b r
Ca rb b r
^ o ii^ te )]
The reader can refer to Jin?s book [1] for the definitions o f b, c, and d terms.
In Table A 1.2 we include the closed form formulae for the second term.
Table A1.2: Closed Formula for
R
u
1
1
1
1
1
1
1
2
2
2
2
2
3
3
3
3
4
4
4
5
5
6
2
3
4
5
6
2
3
4
5
6
3
4
5
6
4
5
6
5
6
6
J{fc2e N r.Nu^/Q
jk*! o_
(f2 2 -fl2 + ftl)/(3 6 0 V )
( 2 f 23 -f2 i-f,3 + fn )/(7 2 0 V )
( 2 f 24 -f2 .-f.4 + fn )/(7 2 0 V )
(f23-f22-2f,3+fl2)/(7 2 0 V )
( f 2 2 - f 24- fi 2+ 2 f u ) /( 7 2 0 V )
(f24-f23 -fl4 + f.3 )/(7 2 0 V )
(f3 3 -fl3 + fll)/(3 6 0 V )
(2 f3 4 -fi3 -fi4 + fn )/(7 2 0 V )
( f 33-f23-fl3+ 2f,2)/(7 2 0 V )
(f23-f34-fl2+ fl4 )/(7 2 0 V )
(f,3-f33-2f,4+ f34)/(720V )
(ft4 -fi4 + f. i)/(3 6 0 V )
( f 34-f24 -fl3 + f.2)/( 7 2 0 V )
( f24 -г w -2 f12+ f l4)/(7 2 0 V )
( f u - f 34 -fi4 + 2 f,3)/(7 2 0 V )
(f33-f23+ f22)/(360V )
( f23- 2 f 34-f22+ f24)/(720V )
( f34-f33-2f24+ 2f23) /(7 2 0 V )
( f22-f24+ f╗4)/(360V )
(f 24- 2 f 23- f w+ f3 4 )/(7 2 0 V )
(г w- f 34+ f33)/(360V )
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where f|j= b,- bj + c*Cj + d,- dj.
The third term which contains V(N r) is evaluated using Gaussian Quadrature [61]:
r
A
A
(A1.4)
*? I
i.e. the integral is equal to the summation o f the values o f the function at n^-point,
multiplied by weight numbers, w,. hi our case nq=3, all w, were equal to 1/3, and the
simplex values were chosen according to Table A 1.3.
Table A1.3 : The Simplex Coordinates for the Gaussian Quadrature Rule, nq=3
I
1
2
3
C.
0.16666666667
0.16666666667
0.66666666667
0.16666666667
0.66666666667
0.16666666667
0.66666666667
0.16666666667
0.16666666667
Formula (2.2S) is also evaluated using the Gaussian Quadrature rule.
71
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX 2
Quasi-Minimal Residual Method
The Problem is to solve the set o f linear equations defined as [A]{x}={b}, where [A]
is a large, sparse, symmetric, complex matrix.
The iterative solver used in the finite element package, is called a no-look-ahead Quasi
Minimal Residual Method (QMR).
The algorithm is presented below [62]:
[M x] = diag{[A \}
COMPUTE {r(o)} = {6} ?[A] {x(o)} for some initial guess {x(o)}.
{V░>} =
{rio)} ; SOLVE [Af, ] {y } =
; Pl = ||y||2 .
Y o= l >*lo= 1
FOR t=l,2, ...
if Pi = 0 or
=0
then method fails
{v(0} = {v(')}/ p i. ; M = M / p ,
SOLVE [M ,]{ y } = M
Si = {v}T{y]
IF ( i=l) THEN
{pw ) = m
ELSE
{p(,)} = m - ( P A / e 1-,){ y M)}
72
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
m =[4{p">>
г,- ? { p il)}T{ p } ; IF Si = 0 method fails
Pi = Sj / S i ; IF Pi = 0 method fail
{v<w?} = о - f t { v <0>
S O L V E
=
{V <'*"[
Pm = IW|2
e , = P m /(7 m |A |) ; Y: = I / J l + O f ; IF Yi = 0 method fails
l^-^-iP iY F K P iY M )
IF (i=l)
ELSE
{?s(l)} =
+ (Om Y,)2 {*(M)}
END IF
{jc(0} = {jc(w )} + {rf(0}
{r(,)} = {r(,-l)} ?{j(0 }
Check Convergence,
END
73
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX 3
Derivations
A3.1 D eriv in g F o rm u la (2.18)
Consider an N-port device, where port i is excited with mode m, while the other ports
are terminated by perfectly matched loads, i.e. only
= 1. The transverse part o f the Inн
field on port j can be then described as (eqn 2.15b):
H ░ '= 2 V . h.
H
(A3.1)
where the summation is over the set o f propagating and evanescent modes. But from
(2. 16),
, a c n _ (Koi _ flo>)
= 2 a ? - r .░ >
(A3.2)
Vnu) is computed from the field distribution using formula (2.17), and hence
/y > = 2 u л -ry ╗ (E л )
(A3.3)
Substitute (A3.3) in (A.3.1) to get,
H░> = г ( 2 a y h . - ( ', ? ( E л ) h . )
H
But all alnJ) are zeros except
(A3.4)
= 1, consequently,
H " ' = 2 5 , h . - ^ ? (e ? )!.,
74
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A3-5)
A3.2 Deriving Formula (3.19)
Define r ? as,
r ?
x
i=l
y
z
[ =l
1=3
(A3.6)
i.e. the x, y , or z component o f a point in 9?3. For a point inside a tetrahedral element, its
x, y , or z component can be expressed in terms o f the simplex coordinates as,
r ' = Z ryгy
(A3.7)
where r j ?s are the coordinates o f the 4 vertices o f the tetrahedron. Using the following
relation,
C4 = ! - г , - C 2- C 3
(A3.8)
the summation o f formula (A3.7) can be reduce to three terms only,
>1
( A
3 -9 )
y=i
Formula (A3.9) is a mapping function from
г 3) to (r?, r 2, r 3). Consequently,
one can derive the corresponding Jacobian matrix, whose entries are defined as,
dr
=
(A3.10)
The inverse Jacobian matrix is then,
(s-l
- * dr"
?
(A 3.ll)
75
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Since
(A3.12)
It follows that,
_____
d rlk d r3
d rlk
3 .
= - z z ( '- o ~ > (?/-* ),
/╗=? ?=>
оrk
= - f y j г ╗
A ( r p - r ^>\
p=, ,=? d rp drk
9
dr"
- it
(A3.13)
For k = l, formula (A3.13) becomes
8 dGa _ ^ V f c
e \ 5C?
drlk drn
d r 'U
9
^
= _d^_e^_
5rp dr"
(A3.14)
Similar results can be obtained for k=2, and 3. For k=4,
A ^ . =_ ^ .v r
drt'dr"
d r ' ^ 4*
44V
_ dC? y dgg
5rp г dr"
(A3.15)
Using again the following relation,
C4 = i - C , - C 2 - C 3
d г 4 = a ? , a g 2 ag,
drp
drp drp drp
in formula (A3.IS) to obtain,
76
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A3.16)
a ag. =
ar/ ar"
ag. y 4
a r ' ar"
(A3.17)
Applying the results from (A3.16) and (A3.17) to all components, i.e.
_a_ a^ a^ ^
r
ar/
(,
a r 1 ? a r 2 ? ar3 '
a^_ ac?d$k ac? ag.
(A3.18)
d r 1 d r 1 ? a r1 d r 2 ? ar' ar3
which can be rewritten in more compact form as
-^ -V г a = -V г*
ar*
?
4 d r1
(A3.19)
A3.3 D eriv in g F o rm u la (3.20)
The volume o f a tetrahedral element can be written as:
(A3.20)
Then formula (3.19) can be written as:
J | = ^ ( л v c 2 - ( v c ,* v c , r
(A3.21)
But,
| ^ = -(6V C i ,(V?J xV C 1))- 2
drk
6^
i .( '7 f ) x V C ,)+
. ar*
6VC2.
av^3
r - x v Ci +
i, a *
6VC2 VC3 x
av^L
dr? J
Substitute (3.19) in (A3.22) to get,
77
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A3.22)
av
i- 2
^ - = ( 6 V л I .(VCJ x V C ,╗
* ? ? ,) +
6 V C , f v л , ^ f - x V C 1] +
6VC 1/
vcj xV
C , |г l '
(A3.23)
For k=l formula (A3.23) becomes,
^ j = ( 6 V C I .(VCJ xVC1╗
1-2
6 V ? ,^ f .( V < ,x V C ,)
['
s v c .- f v c .^ - x v c ,
6V г J V f.xV C 3 0 '
1 Sr7
J
(A3.24)
The first and second terms on the right hand-side o f formula (A3.24) are equal to zero.
Consequently,
dV
^ j = (6VC! .(VCJx V C ,r
( e v c ^ x v c .f ^
= (6VC:.(VCJ x V f , ) ) - ' ^ j
(A3.25)
A similar result is obtained for k=2 and k=3. For k=4,
78
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
avr = (6VC2.(V C,*V C,))-
Sri
6 v c 2. ( v c 4| г ) - x v c , l
6 V г 2. VC3xV C 4
(A3.26)
S r'
However,
vC
4 = -vc, -
(A3.27)
v c 2- v c 3
Substitute (A3.27) in (A3.26),
| ^ = ( 6 v c ! . ( v ^ x v<;,))-i
dr4
- 6 ^ r ( VC, .(V f, x V f , ) + V f 2.(Vг2 x V ? ,) + VC, .(VC, X V C ,))
- 6^HvC2'(vc, x vc,)+vCj.(vc2xvc,)+vc,.(vCj x vc,))
-л ^ f(v c 2-(vcJ xvc,)+vc2.(vcj xVCj)+vC;.(vCj xvc,))l
(A3.28)
The first, third, fourth, fifth, eighth, and ninth terms are all equal to zero. So formula
(A3.28) can be rewritten as
dV
i - F = (6VC! .(V C ,x V f1))-!
( - 6 V C i .(VCJ x V C , ( ^ + f f + |
f )]
= (6VC2.(VC j x VC,))-!
-6 V C 2.(VC, x V C , / - |L f
(A3.29)
=v { ^ )
\d r? )
79
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A3.5 Deriving Formula (3.21)
The basic function N ,is equal to,
(A3.30)
N r = C av
Applying the curl operation,
V x N , = V x f c .V < ,) - V x f c .V C .)
(A3.31)
Consider the following identity operation,
V x(i^a) = V i ^ x A +\p V x A
(A3.32)
If A is the gradient o f a scalar function, as in (A3.31), then the second term in (A3.32)
can be eliminated using the identity operation,
V x (V \f/) = 0
(A3.33)
Consequently, (A3.31) becomes,
V xN r =VCa xVC6 -VC*xVCa
(A3.34)
V x N , =2(V C ?xV <;4)
(A3.35)
Differentiating with respect to r[ ,
avxN
drl
SVxN
a r;
which is formula (3.20).
=
(A3.36)
2
^ x v c
dr'
o rI
(A3.37)
d r?
~
A3.2 D eriv in g F o rm u la (3.22)
The basic function N r is equal to,
80
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A3.38)
N r = C aV г , - c * v c a
Differentiating with respect to r / ,
(A3.39)
= rS a ^ - L / - rS . *
Sr/
^
?
о Ll
ar/
E г i . +r r
dr'
^
dr?
╗ r
ar/
5Vг5__
a r'
^
d r1
avc,
a r' J
(A3.40)
which is formula (3.21).
A 3.6 D eriv in g F o rm u la (4.1)
Starting first with the derivative o f the amplitude term, we find that
(A3.41)
dG
'd G
1 1
However,
e\s\_dM,Y +(s,Y
dG
dG
4s,f)l M
= t/
dG
m
_
y)
dG
.y m y
s m
+ S iM i
r dG
1 dG
MY MY
where S ? S r + j S i . Substitute (A3.42) in (A3.41) to obtain the first term o f (4.1).
Next is the phase term,
81
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A3.42)
However,
dZS
dG
5 ta n ??(? )
S.
dG
r S, '
1
2 dG
'r J
5,. a f e ) | 1 d jS ;))
v S r2 dG
S r dG
-S
Sl so
% ?,)
S?~an
W'
where S = S r + jS i . Substitute (A3.44) in (A3.43) to obtain formula (4.3).
82
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
REFERENCES
1. Jin J., The Finite Element M ethod in Electromagnetic. New York: Wiley, 1993.
2. Lee H., and Itoh T., ?A Systematic Optimum Design o f Waveguide-to-Microstrip
Transition,? IEEE Trans. on M icrowave Theory and Tech., Vol. 45, No. 5, May 1997,
pp. 803-809.
3. Garcia P., and Webb J. P., ?Optimization o f Planar Devices b y the Finite Element
method,? IEEE Trans, on M icrowave Theory and Tech., Vol. 38, No. 1, Jan. 1990,
pp. 48-53.
4. Akel H., Webb J.P., ?Design Sensitivities for scattering-matrix calculation with
tetrahedral edge elements, ? presented at the Conference on the Computation o f
Electromagnetic Fields (COMPUMAG), Sapporo, Japan, October 25-28, 1999.
5. Bandler J. W., ?Optimization Methods for Computer Aided Design,? IEEE Trans.
M icrowave Theory Tech., Vol. 17, August 1969, pp. 533-551.
6. Garcia P., Optimization o f H -plane Junctions Using Finite Elem ents, Master Thesis,
McGill University, Canada, 1989.
7. Rechenberg I., Evolutionsstrategie,
Prinzipien
der Brologishen
Optimierung Technischer
Evolution.
Stuttgard-Bad
Systeme Nach
Cannstacht.
Germany:
Frommann Holzboog, 1973.
8. Tucholke U., Arndt F., Weirdt T., ?Field Theory Design o f Square Iris Polarizers,?
IE EE Trans. Microwave Theory Tech., Vol. 34, N o.l, Jan 1986, pp. 156-159.
83
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
9. Papziner U., Arndt F., ?Field Theoretical computer-aided design o f rectangular and
circular iris coupled rectangular or circular waveguide cavity filters,? IEEE Trans.
Microwave Theory Tech., Vol. 41, No. 3, Mar. 1993, pp. 462-471
10. Keller R. Arndt F., ?Rigorous modal analysis o f the symmetric rectangular iris in
circular waveguides,? IEEE Microwave Guided Wave Letter, Vol. 3, June 1993, pp.
1850-187.
11. Krauss P., and Arndt F., ?Rigorous mode-matching m ethod for the modal analysis o f
the T-junction circular to sidecoupled rectangular waveguide,? M TT-S Int.
Microwave Symp. D ig., Orlando, FL, May 1995, pp. 1355-1358.
12. Sieverding T., Bomemann J., and Amdt F., ?Rigorous design o f sidewall aperture
couplers,? M TT-S Int. M icrowave Symp. Dig., Vol. 2, June 1993, pp. 761-764.
13. Dittloff J., and A m dt F., ?Rigorous field theory design o f millimeter wave E-plane
integrated circuit multiplexers,? IEEE Trans. M icrowave Theory Tech., Vol. 37, Feb.
1989, pp. 340-350.
14. Am dt F., Tucholke U., and Wriedt T., ?Computer-Optimized multisection
transformers between rectangular waveguides o f adjacent frequency bands,? IE EE
Trans. M icrowave Theory Tech., Vol. 32, Nov. 1984, pp. 1479-1484.
15. Am dt F., Koch B., Orlok H. J., and Schroeder N., ?Field Theory design o f
rectangular waveguide broad-wall metal-insert slot couplers for millimeter-wave
applications,? IE EE Trans. M icrowave Theory Tech., Vol. 33, Feb. 1985, pp. 95-104.
16. Am dt F., Sieverding Th., W olf T., Papziner U., ?Optimization oriented design o f
rectangular and circular waveguide components using efficient mode-matching
84
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
simulators in commercial circuit CAD tools,? Int. J. M icrowave MM-Wave Comp.
Aided Eng., Vol. 7, Jan. 1996, pp. 37-51.
17. Rieter J. M., Arndt F., ?Rigorous analysis o f arbitrarily shaped H- and E-plane
discontinuities in rectangular waveguides by a fiill-wave boundary contour modematching method,? IEEE Trans. Microwave Theory Tech., Vol. 43, Apr. 1995, pp.
796-801
18. Amdt F., Beyer R., Sieverding Th., and W olf T., ?Automated Design o f waveguide
components using hybrid mode-matching/numerical em
building blocks
in
optimization-oriented CAD frame works- state-of-the-art and recent advances,? IEEE
Trans. M icrowave Theory Tech., Vol. 45, No.5, May 1997, pp. 747-759.
19. Zhang M., Weiland T., ?Automated Optimization o f a Waveguide Microstrip
Transition Using an EM Optimization Driver,? IEEE Trans. M icrowave Theory
Tech., Vol. 45, No. 5, M ay 1997, pp. 861-864,.
20. Holland J. H., ?Genetic Algorithm,? Scientific American, July 1992,44-50
21. John A., Jansen R. H., ?Evolutionary Generation o f M(MIC) Component Shapes
Using 2.5D EM Simulation and Discrete Genetic Optimization,? IE EE MTT-S Intern.
Microwave Symp. Digest, V ol.l, 1996, pp. 745-748.
22. Jones E.A., Joines W.T., ?Design o f Yagi-Uda Antennas Using Genetic Algorithms,?
IEEE Trans, on Antennas and Propagation, Vol. 45, No. 9, Sep. 1997, pp.1386-1392.
23. Altshuler E. E. and Linden D.S., ?Design o f Loaded Monopole Having Hemiн
spherical Coverage Using a Genetic Algorithm,? IEEE Trans. On Antenna and
Propagation, Vol. 45, No. 1, January 1997, pp. 1-4.
85
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
24. Okoshi T., Imai T., Ito K., ?Computer-Oriented Synthesis o f Optimum Circuit Pattern
o f 3-dB Hybrid Ring by the Planar Circuit Approach,? IEEE Trans. M icrowave
Theory Tech. Vol. 29, No.3, March 1981, pp. 194-202.
25. Miyoshi T., Shinhama T., ?Fully Computer-Aided Synthesis o f Planar Circulator,?
IEEE Trans. M icrow ave Theory Tech., Vol. 34, No.2, Feb 1986, pp. 294-297.
26. Aloss J. T., Guglielmi M ., ?Simple and Effective EM-Based Optimization Procedure
for Microwave Filters,? IEEE Trans. M icrowave Theory Tech., Vol. 45, No.5, 856859, May 1997.
27. Hadamard J., Lecons sur Le Calcul des Variations, Librairie Scientifique A. Hermann
et Fils, Paris 1910.
28. Lee H., Lee S., Juang H., Hahn S., ?An Optimum Design Method for Eigenvalue
Problem,? IEEE Trans. M agnetics, Vol.32, No.3, M ay 1996, pp. 1246-1249.
29. Rekanos I. J., Tsiboukis T. D., ?Electromagnetic Field Inversion Using QuickProp
Method,? IEEE Trans. M icrowave Theory Tech., Vol. 33, No. 2, Mar 1997, pp.
1872-1875
30. Suzuki M., Hosono T., ?Optimum Sectional Shape o f Dominant Mode Waveguide,?
IEEE Trans. M icrowave Theory Tech., Vol. 34, N o .l, Oct 1983, pp. 156-159,
31. Zhang M., ?The Gradient Associated Conjugate Direction Method,? Appl.
Computational Electrom agnetics Soc. J., Nov. 1996.
32. Collin R.E., Foundations fo r Microwave Engineering, McGraw Hill, New York,
1966.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
33. Webb J. P., and Kanellopoulos, ?Absorbing Boundary Conditions for the finite
element solution o f the vector wave equation, Microwave Opt. Tech. Letter, Vol. 2,
Oct. 1989, pp. 370-372
34. Jin J. M., Volakis J. L., and Liepa V. V., Fictitious Absorber for truncating finite
element meshes in scattering, IEE Proc. H, Vol. 139, Oct. 1992, pp. 472-476.
35. Peterson A. F., ?Accuracy o f 3D radiation boundary conditions for use with the
vector Helmholtz equation,? IEEE Trans. Antennas and Propagation, Vol. 40, Mar
1992, pp. 351-355.
36. Berenger J. P., ?A Perfectly matched layer for the absorbing o f electromagnetic
waves,? J. Comp. Phy., Vol. 114, O ct 1994, pp.185-200.
37. Katz D. S., Thiele E.T., and A. Taflove, ?Validation and extension to threedimensions o f the Berenger PM L absorbing boundary conditions for FD-TD
meshes,? IEEE Microwave and Guided Wave Letters, Vol. 4, No. 8, August 1994,
pp.268-270.
38. Zhao L., and Cancellaris, ?A general Approach for the Development o f Unsplit-Field
Time Domain Implementation o f Perfectly Matched Layers for FD-TD Grid
Truncuation,? IEEE Trans. Microwave Guided Wave Letter. Vol. 6, 1996, pp.209211 .
39. Kingsland D.M., Gang J., Volakis J.L., and Lee J.F. ?Anistropic artificial absorber
for truncuating finite-element meshes, ?IEEE Trans. Antennas and Propagation, Vol.
46, July 1996, pp. 975-982.
with permission of the copyright owner. Further reproduction prohibited without permission.
40. Stupfel B., ?Vector wave equation, second and third order absorbing boundary
conditions, ? IEEE Trans. Antennas and Propagation, Vol. 45, March 1997, pp. 487492.
41. Gavrilovic M., Webb J.P., ?A port element for p-adaptive S-parameter calculation,?
IEEE Trans. Magnetics, Vol. 35, No. 3, May 1999, pp. 1530-1533.
42. Silvester P.P., and Ferrari R. L., Finite Elements for Electrical Engineers, Cambridge
University press, 1983.
43. Uher J., Bomemann J., and Rosenberg U., Waveguide Components fo r Antennas
Feed Systems: Theory and CAD, Norwood, MA: Artech House 1993.
44. Marcuvitz N., Waveguide handbook. New York: McGraw-Hill, 1951, pp.80-84.
45. Avriel M., Nonlinear Programming: Analysis and Methods, Prentice Hall, 1976.
46. Dixon L.C.W., Nonlinear optimization, The English Universities Press Ltd, 1972.
47. Fletcher R., Practical M ethods o f Optimization ? Unconstrained Optimization, V ol.l,
John Wiley & Sons, 1980.
48. Scales L.E., Introduction to Non-Linear Optimization, Springer-Verlag New York
Inc., 1985.
49. Ledermann W., Churchhouse R.F., Handbook o f Applicable Mathematics ?
Numerical Analysis, Vol. HI, John Wiley & Sons, 1981.
50. Biggs M.C., ?A Note on Minimization Algorithms which make use o f Non-Quadratic
Properties o f the Objective Function,? Vol. 12, J. Inst. M ath. Applies., 1973, pp. 337338.
51. Pierre D.A., Optimization Theory with Applications, New York, Dover, 1986.
52. Walsh G. R., Methods o f Optimization, John Wiley & Sons, 1985.
88
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
53. Powell
?How Bad are the BFGS and DFP Methods W hen the Objective
Function is Quadratic,? vol.34, M ath. Prog., 1986, pp. 34-47.
54. Choufani M.P., A 3D Extrusion M esh Generator fo r P-Type F inite Elem ent Methods,
Master Thesis, McGill University, Canada, 1997.
55. Welt D., Webb J., ?Finite Element Analysis o f Dielectric Waveguides with Curved
Boundaries,? IEEE Trans. M icrowave Theory Tech., Vol. 33, No. 7, July 1985, pp.
576-585.
56. Zienkiewicz O.C., The Finite Element Method, McGraw-Hill, London, 1977.
57. Crowley C., Silvester P.P., Hurwitz H. Jr., ?Covariant Projection Elements for 3D
Vector Field Problems,? IEEE Trans. M icrowave Theory Tech., Vol. 24, N o .l, Jan.
1988, pp. 397-400.
58. Press W.H., Teukolsky S.A., Vetterling W.T., and Flannery B.P., Num erical Recipes
in C, Cambridge University Press, New York, 1988.
59. Marrocco A., Pironneau O., ?Optimum Design with Lagrangian Finite Elements:
Design o f an Electromagnet, ? Comp. Methods. Appl. Mech. Eng., Vol. 15, 1978, pp.
277-308.
60. Shyy Y.K., Fleury, ?Shape Optimal Design Using High-Order Elements,? Computer
M ethods in Applied M echanics and Engineering, Vol. 71,1988, pp. 99-116.
61. Dunavani D. A., ?High Degree Efficient Symmetrical Gaussian Quadrature Rules for
the Triangles,? Int. Journal fo r Num erical M ethods in Engineering, Vol. 21, 1985,
pp. 1129-1148.
62. Barrett R., Berry M, Chan T. F., Demmel J., Donato J., Dongarra J., Eijkhout V.,
Pozo R., Romine C., and Van der Vorts H, Templates fo r the Solution o f Linear
89
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
S ystem s: B uilding Blocks fo r Iterative Methods, 2nd edition, SIAM, Philadelphia, PA
, 1994.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
we get
= 2V gt x
drl
aw .
U r'
f
(3.21)
Sa
?
- ^ r ^ C ,b
3 r#
S*
ac╗
where N r and N uare,
27
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.22)
Nr - c .v c ,- e * v c .
(3.23)
(3.24)
Introducing the following definitions:
(3.25)
C a, =
dr1
(3.26)
d+
(3-27)
= 5 > ? ,c ?
/= I
(3.28)
e * * = 4 (V f. x V f 4) .( V ^ x V f ? ) = 4
- <*?</?)
formula (3.18) can now be rewritten in simple form as,
SAT
dr/
I
6V
A rja
(3.29)
'abed
?
г r ( ^ o c C^bd
I ad ^ b e
^b e ^ ad + ^ b d ^ ac )
6V
+
r*r
( C bl e kacd
C al e kbcd + C dl e kcab
C cl e kdab)
?6 V k a er(cal Ibc d u ?ca[ I m d kc ?cbl Iac d kd + cbl I ad d kc
+ Cct Ida d u, ?ccl I db d ka ?cdl I ca dkb + cdl Icb d^ )]
Using (3.17) and (3.29), it is straightforward and inexpensive to obtain the derivatives
o f [K\ for any geometrical parameter g.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
Optimization of Microwave Devices
4.1 In tro d u c tio n
The main purpose o f optimizing microwave devices is to meet design objectives such
as return loss, insertion loss, and coupling. There are many techniques for optimizing
microwave devices, but, since we have gradient information available at very little cost,
only optimization techniques that make use o f gradient information are o f interest.
In section 4.2 we present a general formula for a cost functions used in microwave
design, and the corresponding gradient formula. In section 4.3, we describe the conditions
for a point to be a local minimum, then in section 4.4 we discuss, in general, the theory o f
the gradient methods. Different kinds o f optimization techniques are discussed in section
4.5, while in section 4.6 we present the Quasi-Newton methods, which are the methods
used in this research. In section 4.7 the Line Search method is presented in details.
4 .2 C o s t F u n c tio n f o r M icrow ave D e v ic e s
A general cost function o f an N-port microwave device is [6]:
(4.1)
where a ',
and P ~ are weight functions, S - { f t ) are the desired scattering parameters,
and the frequency range has been discretized into Ah frequencies. The a ~ and p~
weights control the emphasis on the magnitude and phase repectively o f the scattering
parameters. From this general formula, one can also deduce the formula o f the gradient o f
29
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the cost function with respect to any geometrical parameter, say g. The derivative o f the
cost function is a function o f the scattering parameters and the derivative o f the amplitude
and/or the derivative o f the phase o f the scattering parameters. The derivative o f the
amplitude is equal to:
r f |W < ) |
dg
R efc, ( / , )}R e(8,?; ( / ' *1 +
dg
( / , )}lro(35╗ ( / ' }
dg
(4.2)
? W /)I
while the derivative o f the phase is:
(4.3)
<US,
dg
In terms o f these, the general formula for the derivative o f the cost function is:
f-? ? ?
(J , )
dg
f c ,( .f, )| - | s , ( / , )|)+
d Z S jX fi) (
*
(4'4>
\
4 .3 C o n d itio n s fo r a L o cal M inim um
In general, the cost function C can be written as a function o f a vector {g} o f
geometrical parameters. Let {g}={g*} be a local minimum, and C be continuously
differentiable at {g*}. Then the following condition must be satisfied [45-48]:
V C \g* = 0
(4.5)
i.e. the slope o f C is zero at {g*}. However, satisfying (4.5) does not necessarily
mean that {g*} is a local minimum; it could be a local maximum. To obtain minima only,
an additional criterion is needed. The new criterion is derived from the second derivative
30
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
o f the cost function C.
Let C be twice continuously differentiable at {g*}. Then, one can expand C in the
neighborhood o f {g*} to obtain the formula:
c({g-}+ {a})= cfe*})+
where
)
(46>
is a very small perturbation vector, and the [H] matrix is the Hessian o f C:
(4.7)
>J
Sgi8gj
g
*
Higher order terms are neglected in formula (4.6). For {g*} to be minimum, the left-hand
side o f (4.6) should be greater than the cost function at {g*}. Consequently, the additional
required criterion for minima is:
(4.8)
{5}r [tf]{5 }> 0
V {S }* 0
Formula (4.8) can be satisfied if [//] is a positive definite Hessian matrix.
4 .4 M ultivariate G ra d ie n t M e th o d s o f M in im izatio n
In all gradient methods, a vector {g} o f geometrical parameters is iteratively updated
until the objective cost, C, becomes lower than a prespecified threshold. The A+1st update
is o f the following form:
(4.9)
where {p(kJ} is certain search direction, and
is a step size in the direction o f {p(k)}-
Gradient methods use the slope (gradient) o f C at {g?k)} to construct a suitable search
direction. Some gradient methods do not require derivative information, like Powell?s
31
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
method. However, we are only interested in those techniques that make use o f the
derivative information such as, steepest descent, Newton, quasi-Newton, and the
quickprop method.
The gradient methods differ from each other in the w ay the vector {p(k>} is updated.
The value v(k) can be constant, or it can be calculated using single-variable optimization
(line search) technique.
In the steepest descent method, developed by Cauchy, the vector {p(k>} is given by:
{/>"}=-vcrfe};
(4I0)
i.e. the direction o f search is along the line o f maximum rate o f change o f C. The steepest
descent method has a linear convergence rate for a quadratic function. This is due to the
fact that the second derivative is not used. Thus, the method is not usually recommended
for general applications.
On the other hand, the Newton Raphson method uses the second derivative
information to update the search direction vector {p(k>}- The formula used is:
{p'l))= [H m y 'c u g } )
<4 -n>
where [//] is the Hessian matrix o f equation (4.7) which includes the second derivative
information.
The main disadvantage o f the Newton-Raphson method is the cost o f calculating [//],
which rises dramatically as the number o f variables increases. In addition, singular
Hessians, indefinite matrices, and overshooting could be possible sources o f problems to
this technique.
For conjugate gradient methods, the search direction {p?}, satisfies the conjugacy
condition [48]:
32
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
m ^n
(4.12)
One can show that {p(k)} vector can be generated iteratively without constructing the
Hessian matrix [//], using formula:
(4.13)
{p(m,}= -W C (k>+
Conjugate gradient methods are highly efficient when the number o f variables is large,
and no information is available on the second derivatives, or it is very difficult to obtain
them. An improved convergence rate can be accomplished through exact line searches to
obtain the optimum step size v(k) in the direction o f {p(k)}.
In the quickprop method [29], the cost function is approximated by a parabola which
has arms open upwards. The assumption made is that the slope o f the cost with respect to
a variable is not affected by the change in the values o f all other variables in the model o f
concern. Then, one can easily show that:
(4.14)
Rekanos showed that the QP technique is 3-4 times faster than the steepest descent
method. No comparisons with Newton or quasi-Newton techniques are available yet.
4 .5 Q u a si-N e w to n M e th o d s
To avoid the calculations o f the exact Hessian matrix in the Newton-Raphson method,
the quasi-Newton methods use an iterative scheme to approximate it. The gradient and the
33
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
cost function are used in the calculations, which take place at each iteration.
The search direction {p(k)} is found by:
(4.15)
{pw } = -[A fa ) ]{vCa ) }
and the [M\ matrix is updated using the following formulae [45][47-51][52]:
(4.16a)
(4.16b)
W
J
W
w
\ r )
(4.16c)
2(2(5+ l - 3 i j )
Where
{5(i)}r { v c (i_I)} ?
C (t) - C (*~l)
{/*>}= {vc(t)}- {vc(?_1)}
and,
Formula (4.16) is known as the Broyden-Fletcher-Goldfrab-Shanno (BFGS) method.
There is another quasi-Newton method called the Davidon-Fletcher-Powell (DFP)
method. Powell [53] showed that BFGS performs better than DFP in minimizing
quadratic functions o f two variables, which suggests that BFGS should be used also for
34
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
general non-linear functions.
To improve the overall performance o f the algorithm, a good initial value for [M] is
necessary [54]. A convenient choice is to set the initial value to v(0)[/o], where [70] is the
identity matrix. v(0) is the initial step size.
When using a line search, the value o f the step size J k) converges to 1.
The algorithm o f the quasi-Newton module has the following steps for minimization:
1. Begin with [M]=[/a].
2. Using the initial geometry, calculate the cost and its derivative.
3. Call the line search module, and perform an approximate line search along
{/>(0)}= -{V C (0)}. The line search module returns: the new step v(0), the Cost and its
derivative at v(0)
4. Set [A/0)]=
v (0)[ / o].
5. Evaluate formulae (4.16).
6. Evaluate formula (4.15).
7. Update the geometrical parameters using formula (4.9).
8. Call Line Search which returns under three conditions:
8.1 If the cost or its derivative are better than the target values, return.
8.2 If the geometrical parameters were not changed, return.
8.3 Otherwise proceed to step 9.
9. Redo 5-8.
4 .6 L ine S e a rc h
To improve the convergence rate o f any gradient method, a line search is used to
update the value o f the step size v/*; at each iteration. The goal o f the line search is to find
35
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the v╗ which makes C({g(k)}+ \fk> {p(k)}) locally minimal, i.e.
(4.17)
Using the chain rule, this becomes
{p(i)}r { v c (*+l,} = 0
Formula (4.18) represents the condition for the termination o f an exact search.
(4.18)
The line search module is composed o f two consecutive steps: establishing bounds for the
search, and allocating a minimum along the {p(k)} vector.
4.6.1 Establishing Bounds for the Search
In order to find a minimum, the line search module needs to identify an interval
bounding at least one local minimum. Making use o f the derivative information,
Davidon?s polynomial extrapolation method is used to allocate the upper or lower bound
o f this interval.
Let a be a given lower or upper bound for v, depending on the sign o f C '. Then an
estimate for the minimum is:
(4.19)
b =a +2
where C? is the cost at a, and C, is the target cost (a user-defined value). C ' is the slope
o f the cost function with respect to v at point a.
Formula (4.19) can be applied iteratively : the program starts with a given a, Ca, and
C ', then estimates the upper bound o f the interval i f C ' is negative, or the lower bound
o f the interval if C ' is positive. If the sign o f the slope o f the cost function at the new
position b is the same as at a, the program changes the value o f a to b, and seeks another
value for b. The process stops when the slope changes its sign.
36
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The algorithm o f the interval step includes the following steps:
1. Start with a=0.
2. I f the slope o f the cost function with respect to v is negative then,
2.1. Evaluate b using formula (4.19).
2.2. Evaluate formula (4.9) using v=b.
2.3. Solve the problem and calculate the cost function and its derivative.
2.4. I f the cost function is below the target cost function Ct, return.
2.5. I f the derivative o f the slope is negative, then,
Set a~b, Ca=Cb, and C ' = C 'b , then go to 2.1
2.7 Else return.
3. If the slope o f the cost is positive then,
3.1. Set b=a, Cb=Ca., C'b = C'a
3.2. Evaluate a using formula (4.19), where a is replaced w ith b and b with a.
3.3. Evaluate formula (4.9) using v=a.
3.4. Solve the problem and calculate the cost function and its derivative.
3.5. If the cost is below the target cost function C,, return.
3.6. If the derivative o f the slope is positive then,
Set b=a, Cb=Ca, and C'b = C ', then go to 3.1
3.7. Else return.
4.6.2 Locating a Minimum
Once an interval bounding a minimum has been determined, polynomial interpolation
methods can be used to approximate the cost function by a polynomial, from which a
37
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
possible minimum is calculated. In this thesis, cubic interpolation is used. The method
exhibits quadratic convergence.
It can be shown that if the cost function and its derivatives are known at the bounds o f
the interval ([a,b]), the minimum o f a cubic polynomial fitting these values is given by:
(4.20a)
(4.20b)
(4.20c)
w
Formulae (4.20) are used in an iterative process to obtain the real m inim um o f the cost
function within the given interval. For each new ?cubic? m inim um , the cost function and
its derivative are evaluated. I f the derivative o f the cost function is not small enough, then
the calculated point is not a real minimum. However the new information is not wasted,
but used to narrow the interval from one o f the two sides (lower or upper) depending on
the slope o f the cost function at v. This process is repeated for the new interval, until the
slope o f the cost is below certain threshold, or the interval is less than a given intervalthreshold.
The algorithm for the allocation o f a minimum is:
1. Calculate v using formulae (4.20).
2. Evaluate formula (4.9).
3. If the differences between the new and old geometrical parameters are below the
geometrical tolerance, return
4. Adjust the geometry and solve the new model.
5. Calculate the new cost function and its derivative.
38
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
6. If the new cost function or its derivative are better than the target cost Ct return.
7. If the derivative is negative then set a=v, Cfl=Cv, and C ' = C ', and go to 1.
8. Else set b ?v, Cb~Cv, and C ?b = C ', and go to 1.
39
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 5
Computer Program Implementation
5.1 In tro d u c tio n
Using the theory discussed in the previous chapters, a Fortran-99 program was
developed.
The program consists o f two main parts: the optimizer and the 3D finite element
package. The optim izer is made o f two modules: the Quasi-Newton and the Line Search
modules. The 3D finite element solver is made o f four modules: the m esh generator, the
node-to-edge converter, the 3D finite element solver, and the remodeller. The program
interface is via input and output files. The program needs four input files which are :
Geo.dat, 2D.dat, Input.dat, Ports.dat, Mats.dat, and Param.dat. It generates three
temperory files: Tets.dat, Edges.dat, and Results.dat, and the output o f the program is one
file: Tracing.log.
In Section 5.2 we describe the six modules o f the program, then in section 5.3 we
present the different data files used during the run.
5.2 P ro g ra m S tr u c tu r e s
The flow chart shown in Figure 5.1 demonstrates the interconnection between the six
main modules o f the program.
5.2.1 Quasi-Newton Optimize
The method implemented here uses the BFGS technique described in the previous
chapter. The efficiency o f this optim izer depends on how fast and accuratly the cost
40
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
function and the gradient information are obtained.
The Quasi-Newton method is responsible for generating new geometric param eters.
The Quasi-Newton calculates the search direction p, and calls the Line Search module to
obtain the minimum along p.
Start ^
E valuate Cost Function
Call
Initialization
\
?
f
S
Remodeller
j
Return
Mesh
Generator
/?
V.
( a d
'
y
/
Quasi-Ne^on
Module
C all
Nodal-to-Edge
Converter j
Return
Line Search
Module
Call
FEM Solver
Return
Figures. 1 Representative flowchart o f the relation between the six main modules : quasi-Newton, line
search, remodeller, mesh generator, nodal-to-edge converter, and the 3D FEM solver.
5.2.2 Line Search
The Line Search module is composed o f two steps: in the first step Davidson?s quadratic
polynomial extrapolation is used to construct an interval bounding at least one local
41
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
minimum, and in the second one, the cubic polynomial interpolation is implemented to
minimize the cost function along a search direction p.
The Line Search may not perform the two steps all the time, as discussed in the previous
chapter. The user may also choose to discard the use o f the Line Search option.
5.2.3 Mesh Generator
The mesh generator was written by Marc P. Choufani (Master Project) [54]. The user
supplies two data files: a module file and a design parameter file. The program generates
a data file which contains the (x,y,z) coordinates o f all the points, elements material labels
and nodes, and the boundary conditions o f each o f the four faces o f all the elements. The
program uses the Delauny triangulation method to construct the 2D mesh, then uses the
extrusion method to make the final 3D-nodal mesh.
5.2.4 Node-to-Edge Converter
The 3D mesher generates a 3D mesh in which each node is assigned a global number.
However, the 3D finite element solver is a 3D-edge solver. The Node-to-Edge converter
module assigns a global number to each edge in the 3D mesh. The output file, Edges.dat,
contains information on all edges and their associated nodes (two nodes per edge), and all
elements and their associated edges (6 edges per element).
Some elements have the number zero assigned to some o f their edges, this indicates
that theses edges are located on perfectly conducted surfaces, i.e., homogeneous Dirichlet
Boundaries.
5.2.5 3D Finite Element Solver
This module performs four different operations: constructing the [K] and [R] matrices,
solving the linear set o f equations, calculating the scattering parameters, and finally
calculating the derivat
Документ
Категория
Без категории
Просмотров
0
Размер файла
2 955 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа