# 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

1/--страниц