Novel Techniques for the Design and Characterization of Electromagnetic Devices with Application to Multilayer Structures and Waveguide Filters Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classified information. Daniel Lee Faircloth Certificate of Approval: Sadasiva M. Rao Professor Electrical and Computer Engineering Michael E. Baginski, Chair Associate Professor Electrical and Computer Engineering Lloyd S. Riggs Professor Electrical and Computer Engineering Soo-Young Lee Professor Electrical and Computer Engineering Stephen L. McFarland Dean, Graduate School Novel Techniques for the Design and Characterization of Electromagnetic Devices with Application to Multilayer Structures and Waveguide Filters Daniel Lee Faircloth A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama May 11, 2006 Novel Techniques for the Design and Characterization of Electromagnetic Devices with Application to Multilayer Structures and Waveguide Filters Daniel Lee Faircloth Permission is granted to Auburn University to make copies of this thesis at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Vita Daniel Lee Faircloth, son of Sam A. Faircloth and Carlene K. Wade, was born August 13, 1981 in Opelika, Alabama. He graduated from the Alabama School of Fine Arts in 1999. Beginning in the Fall of 1999, he attended Auburn University on a number of scholarships and graduated in August 2002 with a Bachelors of Electrical Engineering degree. In the same month, he entered graduate school at Auburn University after accepting a NASA Graduate Student Researchers Program Fellowship supported through the Electromagnetics Research Branch of Langley Research Center in Hampton, Virginia. In December 2003, he received the Master of Science degree in Electrical Engineering. He is married to Lisa P. Faircloth with whom he has a son Joseph A. Faircloth. iv Dissertation Abstract Novel Techniques for the Design and Characterization of Electromagnetic Devices with Application to Multilayer Structures and Waveguide Filters Daniel Lee Faircloth Doctor of Philosophy, May 11, 2006 (M.S., Auburn University, 2003) (B.S., Auburn University, 2002) 158 Typed Pages Directed by Michael E. Baginksi In this work, optimization methods are presented for the analysis of inverse and design problems in electromagnetics. The current and future demands placed on design engineers necessitate the hybridization of forward solvers (e.g., computational and empir- ical methods) with novel optimization algorithms to yield engineering development and analysis tools which are largely autonomous and not restricted to the familiar principles employed in traditional design work. Two specific problems are addressed to demonstrate the capabilities of the proposed methods. First, several optimization algorithms (i.e., Sequential Quadratic Programming, the Genetic Algorithm, and Particle Swarm Optimization) are presented for the estimation of complex constitutive parameters of multilayered materials. Using X band waveguide S-parameter measurements, the complex constitutive parameters of each individual layer are extracted. The results are compared to measurements as well as those of single layer techniques which estimate the constitutive parameters of individual materials. v The second problem addressed is the automated design of waveguide filters. Initially, a study is conducted into whether dielectric slabs randomly doped with conducting in- clusions such as short, thin wires or thin patches could yield useful frequency dependent reflection and transmission behaviors when placed inside a waveguide. Results obtained by placing conducting patches on the slab?s surface were found promising. Therefore, op- timization techniques were then employed to find the appropriate arrangement of patches of the dielectric?s surface so that the resulting transmission response closely matched the response specified by the user. Results of this study were verified by fabrication and measurement for X band filters and, in all cases, found to be in excellent agreement. vi Acknowledgments The author would like to express his sincere gratitude to Dr. Michael E. Baginski and Dr. Manohar D. Deshpande for supporting him throughout his research. He would also like to extend his thanks to Dr. Sadasiva M. Rao, Dr. Lloyd S. Riggs, Dr. Stuart M. Wentworth, and Dr. Soo-Young Lee for their friendship and guidance. To his parents, he thanks them for always pushing him in the right direction and for all the years of encouragement. He thanks his wife Lisa and son Joey for always understanding the late nights and giving him the constant love and support he needed. Finally, to the Lord Jesus Christ, the author thanks Him for the blessings and gifts he has been given. vii Style manual or journal used IEEE Transactions on Microwave Theory and Techniques (together with the style known as ?auphd?). Computer software used The document preparation package TeXnicCenter and LATEX with the style-file auphd.sty. Illustrations prepared using Canvas X. viii Table of Contents List of Figures xi List of Tables xiv 1 Introduction 1 1.1 Electromagnetic Optimization . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Overview and Historical Perspective . . . . . . . . . . . . . . . . . . . . . 5 1.3.1 Optimization Techniques . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.2 Constitutive Parameter Extraction . . . . . . . . . . . . . . . . . . 11 1.3.3 Waveguide Filter Optimization . . . . . . . . . . . . . . . . . . . . 13 1.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2 Complex Permittivity Extraction 15 2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Error Function Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4 Optimization Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.1 Sequential Quadratic Programming . . . . . . . . . . . . . . . . . . 20 2.4.2 Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4.3 Particle Swarm Optimization . . . . . . . . . . . . . . . . . . . . . 27 2.5 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.6 Measurements Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.6.1 Single Layer Measurements . . . . . . . . . . . . . . . . . . . . . . 36 2.6.2 Two Layer Measurements . . . . . . . . . . . . . . . . . . . . . . . 44 2.6.3 Three Layer Measurements . . . . . . . . . . . . . . . . . . . . . . 48 2.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3 Complex Constitutive Parameter Extraction 52 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.2 Extraction Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.2.1 Wolfson-Wentworth Method . . . . . . . . . . . . . . . . . . . . . . 53 3.2.2 CP Extraction Method Modifications . . . . . . . . . . . . . . . . 53 3.3 Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.1 Single Layer Measurements . . . . . . . . . . . . . . . . . . . . . . 57 3.4.2 Two Layer Measurements . . . . . . . . . . . . . . . . . . . . . . . 65 ix 3.4.3 Three Layer Measurements . . . . . . . . . . . . . . . . . . . . . . 69 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4 Waveguide Filters Using Conductor Doping of Dielectrics 71 4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 Numerical Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3 Finite Element Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.1 Absorbing Boundary Condition . . . . . . . . . . . . . . . . . . . . 79 4.3.2 Mode-Matching Domain Truncation . . . . . . . . . . . . . . . . . 82 4.3.3 Scattering parameters . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.3.4 Asymptotic Waveform Evaluation . . . . . . . . . . . . . . . . . . 88 4.4 Method Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5 Waveguide Filter Optimization Using Surface Patches 106 5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.2 Numerical Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.2.2 Optimization Technique . . . . . . . . . . . . . . . . . . . . . . . . 108 5.2.3 Practical Considerations . . . . . . . . . . . . . . . . . . . . . . . . 109 5.2.4 Forward Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.3.1 Notch Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.3.2 Bandpass Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.3.3 Low Pass Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.3.4 High Pass Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6 Conclusions 126 6.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.1.1 Constitutive Parameter Extraction . . . . . . . . . . . . . . . . . . 128 6.1.2 Waveguide Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.1.3 Extension To Antenna Design . . . . . . . . . . . . . . . . . . . . . 129 Bibliography 131 x List of Figures 1.1 Rastrigrin?s function which has many local minima and a global minimum located at (0,0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Error surface well-suited to local optimization problems. . . . . . . . . . . 8 2.1 Rectangular waveguide loaded with n-layer sample. . . . . . . . . . . . . . 17 2.2 Flowchart illustrating the Genetic Algorithm procedure. . . . . . . . . . . 23 2.3 Crossover between two parents using binary string encoding. . . . . . . . 25 2.4 Comparison of convergence for the traditional GA and GA?s using three redundancy removal methods. The convergence rate of PSO is also shown. 28 2.5 Illustration of the neighborhood structure implemented in PSO. Each par- ticle is connected to its two neighbors to the left and right. . . . . . . . . 30 2.6 Flowchart illustrating the PSO procedure. . . . . . . . . . . . . . . . . . . 33 2.7 Magnitude of S11 and S21 for single layer Bakelite sample. The calculated S-parameters are generated using the GA with error function (2.9). . . . . 37 2.8 Magnitude of S11 and S21 for single layer Garlock Rubber sample. The cal- culated S-parameters are generated using the GA with error function (2.9). 38 3.1 Comparison of MPSQP and Wolfson-Wentworth method extracted com- plex permittivity values for F125 sample. . . . . . . . . . . . . . . . . . . 61 3.2 Comparison of MPSQP and Wolfson-Wentworth method extracted com- plex permeability values for F125 sample. . . . . . . . . . . . . . . . . . . 62 3.3 Measured and generated S-parameter magnitudes for the F125 sample. . . 63 3.4 Measured and generated S-parameter phases for the F125 sample. . . . . 64 3.5 Measured and generated S-parameter magnitudes for the F125/Teflon sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 xi 3.6 Measured and generated S-parameter phases for the F125/Teflon sample. 68 4.1 Loaded waveguide geometry with doped slab. Boundary notations for FEM formulation are indicated. . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2 Discretized slab illustrating wires lying along tetrahedral edges. Notice some edges join to form long, erratically shaped wires. Wires may also be located on interior edges but are not shown here. . . . . . . . . . . . . . . 75 4.3 Tetrahedral element showing edge and node numbering. . . . . . . . . . . 78 4.4 Waveguide mesh used for method verification. . . . . . . . . . . . . . . . . 92 4.5 Comparison of FEM and theoretical |S11| for the loaded waveguide. . . . . 93 4.6 Comparison of FEM and theoretical |S21| for the loaded waveguide. . . . . 94 4.7 Comparison of magnitude of S21 for triangular patch doping using present method and HFSS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.8 A wide variation of S21 responses is shown for the cases of randomly embedded wires. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.9 Illustration of the wide range of notch behaviors observed even when wire connections within the dielectric are not allowed. . . . . . . . . . . . . . . 100 4.10 Magnitude of S11 and S21 for typical case of a dielectric slab doped with 50 2 mm wires. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.11 Comparison of magnitude of S21 for typical cases of a dielectric slab doped with 10, 20, and 35 3 mm wires. . . . . . . . . . . . . . . . . . . . . . . . 103 4.12 Comparison of magnitude of S21 for cases of a dielectric slab doped with 20, 30, and 75 triangular patches. . . . . . . . . . . . . . . . . . . . . . . . 104 5.1 Rectangular waveguide with optimized filter. . . . . . . . . . . . . . . . . 107 5.2 Gridded substrate face showing metallization for fabrication and FEM simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.3 Flowchart illustrating the GA procedure with geometry refinement, re- dundancy removal, and high mutation procedures. . . . . . . . . . . . . . 111 xii 5.4 The lower left quarter of the waveguide filter is shown illustrating the rules governing the geometry refinement process. The inset illustrates the problem with corner-connected patches. . . . . . . . . . . . . . . . . . . . 114 5.5 Comparison of ideal data, MM/FEM, HFSS, and measurements. The MM/FEM matches well with the ideal (requested) response. However, the geometry generated during the process yields a different response in actuality as given by HFSS and measurements. . . . . . . . . . . . . . . . 116 5.6 HFSS model used for optimization . . . . . . . . . . . . . . . . . . . . . . 118 5.7 Fabricated samples produced using the patterns generated by the opti- mization procedure. A dime is shown for size reference. . . . . . . . . . . 119 5.8 Comparison of S21 responses of ideal, simulated, and measured notch filter.120 5.9 Comparison of S21 responses of ideal, simulated, and measured bandpass filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.10 Comparison of S21 responses of ideal, simulated, and measured low pass filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.11 Comparison of S21 responses of ideal, simulated, and measured high pass filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 xiii List of Tables 2.1 Complex Permittivities and Thicknesses for the 1-, 2-, and 3-Layer Com- puter Generated Data Sets . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2 Results using ideal S-parameters . . . . . . . . . . . . . . . . . . . . . . . 35 2.3 Extracted Permittivity & Error for Bakelite . . . . . . . . . . . . . . . . . 40 2.4 Extracted Permittivity & Error for Ceramic . . . . . . . . . . . . . . . . . 41 2.5 Extracted Permittivity & Error for Garlock Rubber . . . . . . . . . . . . 42 2.6 Extracted Permittivity & Error for Nano Material . . . . . . . . . . . . . 43 2.7 Extracted Permittivity & Error for Garlock/Bakelite . . . . . . . . . . . . 45 2.8 Extracted Permittivity & Error for Garlock/Ceramic . . . . . . . . . . . . 46 2.9 Extracted Permittivity & Error for Garlock/Nano Material . . . . . . . . 47 2.10 Extracted Permittivity & Error for Nano Material/Garlock/Garlock . . . 49 3.1 Single Layer Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 MPSQP Error for Single Layer Sample . . . . . . . . . . . . . . . . . . . . 60 3.3 Performance Comparison of GA and MPSQP . . . . . . . . . . . . . . . . 60 3.4 Results for Two Layer Samples . . . . . . . . . . . . . . . . . . . . . . . . 66 3.5 Results for F125/F40/Teflon Sample . . . . . . . . . . . . . . . . . . . . . 69 xiv Chapter 1 Introduction Over the last several years, both consumer and military demands for better per- forming yet smaller and cheaper electromagnetic (EM) devices have overcome design engineers? ability to rely solely on traditional design techniques [1]. Specifically, per- formance requirements for EM devices have created the need for new techniques and technologies to be developed in order to meet the ever-increasing demand. Prior to the early 1980?s, much of the design work was conducted via trial and error using empirical or semi-empirical formulations. However, with the significant gains in computational power and simultaneous development of accurate numerical electromagnetic solution techniques in the 1980?s, designers could now reliably simulate many designs before ever moving to the fabrication and measurement stage. Although this development was a significant milestone in the evolution of EM design, the actual design process still relied heavily on trial and error since engineers continued to use their intuition for design mod- ifications. Of late, computational analysis tools and optimization methods have become hybridized which has allowed designers to realize novel devices that would otherwise have been impossible to conceive through traditional methods. This union of analysis and optimization tools has also opened possibilities for reverse engineering processes (inverse problems) such as characterization of complex manufactured devices and novel multilayered materials to be used in Monolithic Microwave Integrated Circuits (MMIC). 1 1.1 Electromagnetic Optimization Historically, optimization has been applied to electromagnetics in two ways: viz. inverse problems [2] and design problems [3]. An inverse problem is the task of using limited information such as scattered fields or S-parameters to obtain the physical char- acteristics of the scatterer. The parameters being solved for often include geometric shape and dimensions, constitutive parameters (CP), position, and orientation. Opti- mization algorithms typically require some initial estimate for the unknown parameters and iterate until values are found which yield within the tolerance limits the same quan- tity as that which was measured. The solution of an inverse problem is typically quite difficult since limited information is available and, therefore, not sufficient for satisfying uniqueness. Hence, there is the possibility that many solutions may return nearly iden- tical information. To help curtail this problem, constraints are usually placed on the physical characteristics of the object thereby reducing the chance of encountering mul- tiple solutions. As an example, consider trying to determine the physical characteristics of a scatterer using only a small sample of the radiated field. Without any assumptions made, this problem lacks sufficient information to yield a meaningful solution. However, if the scatterer shape is assumed known, then the CP may be calculated to a given degree of accuracy. Design optimization problems closely resemble their inverse counterparts with the exception that the information provided to the optimizer is different. Inverse solvers typically utilize measurements, whereas design optimizers are fed design goals as specified by the design engineer. For example, a certain bandwidth characteristic may be desired 2 for a patch antenna, and an optimization can be performed to determine the shape of the patch that will yield the required results [4,5]. 1.2 Problem Statement This research focuses on the development of novel optimization tools to automate the design and characterization process for a wide range of EM engineering problems. Toward this end, a series of specific problems are addressed even though the methods presented may easily be adapted to a variety of similar tasks. Initially, the inverse prob- lem technique is presented in which the complex permittivity of a single layer structure or each individual layer of a multilayer structure are extracted using S-parameter data from a loaded waveguide. Using this method, the complex permittivity of each layer in the structure is extracted using a variety of optimization techniques. One of the fundamental challenges associated with optimization techniques is the proper balance of convergence rate and robustness (i.e., their ability to adequately search the entire solution space). Therefore, the initial focus is placed on methods that improve the convergence character- istics of the algorithms implemented here without sacrificing robustness. In addition to considering different optimization techniques, a thorough investigation is conducted into the effects of different error functions on the convergence characteristics and accuracy of the optimization algorithms. By varying the mathematical form of the error functions as well as the quality and quantity of S-parameter information included, conclusions may be drawn about the maximum number of layers for which a specific error function is effective. Some consideration is also given to the calibration effort required to make 3 accurate S-parameter phase and magnitude measurements since the availability and ac- curacy of this information will have a dramatic effect on each algorithm?s performance. After this analysis, the method is further modified to include the extraction of complex permeability. Additional discussion of the optimization algorithms is given as this higher dimensionality problem requires a more robust procedure. The second primary area of focus pertains to the use of the aforementioned op- timization algorithms to automate the design of EM devices. Specifically, this work addresses the feasibility of obtaining (optimizing) reflection and transmission properties using dielectric slabs as waveguide filters. The filters are analyzed using a hybrid Mode- Matching/Finite Element Method (MM/FEM) which acts as the forward solver within the context of the optimization algorithm and has been shown to be an efficient algo- rithm for loaded waveguide problems [6,7]. Two analyses are considered with regard to the geometry of the filters. In the first, an attempt is made to characterize the feasibil- ity of obtaining useful filter responses by embedding randomly oriented thin conducting wires and other conducting objects (e.g., thin patches and small volumetric inclusions) within the dielectric slab. Such a material would be quite novel since it would possess designable filter characteristics and yet maintain similar visual and structural proper- ties. For the practical applicability of this type of filter, the responses of the randomly doped slabs would be required to maintain a certain level of consistency for given doping levels (densities) and inclusion sizes. Consequently, a thorough analysis is conducted into whether or not this method is both computationally and physically practical. In the second part, desired waveguide filter responses are found by optimizing the shape 4 of conducting patch on the surface of the dielectric slab. This method has practical advantages since fabrication of such a filter is relatively simple. 1.3 Overview and Historical Perspective 1.3.1 Optimization Techniques Optimization refers to the solution of a problem that can be cast in the form f (vectorx) = 0 where vectorx is an n-dimensional vector of parameters, and f is referred to as the error, objective, or fitness function. Depending on the problem to be solved, each of the parameters xn may be continuous or discrete and may be bounded or unbounded. The problem may also be subject to equality constraints of the form ci (vectorx) = 0 and/or inequality constraints of the form gj (vectorx) ? 0. Within the field of optimization, there are generally two classes of methods: local and global (the classes are also referred to as deterministic and stochastic, respectively). Local Optimization Techniques Local optimization methods utilize information about an error function in the vicin- ity of the current iteration?s parameter values. Chronologically, local methods were the first to be developed and utilized for engineering development since they can usually locate the optimal solution with a minimal number of error function evaluations depend- ing on the initial guess and complexity of the algorithm. Generally, these techniques are subdivided into two categories: direct search methods and gradient-based methods [8]. Direct search methods were first developed in the 1950?s [9] and received much attention due to the fact that they do not explicitly utilize information about the gradient 5 of the error function. At that time, this attribute was particularly beneficial since the numerical computation of gradients was difficult due to computational limitations. These were also the only methods available for optimization of discontinuous search spaces. However, they were all but abandoned by the optimization community beginning in the early 1970?s for a number of reasons outlined by Swann including the fact that there were no mathematical proofs that such methods would ever converge to a minimum [10]. Within the last 10 years, however, direct search methods have regained popularity as certain proofs are now available guaranteeing convergence [11?14]. Gradient-based optimization, as the name implies, relies on the gradient of the er- ror function to formulate a search direction from which the function minimum can be found. The Steepest Descent Methods [15] and others like them, in their simplest forms, require the search space be continuous in its first derivative. Higher order methods such as Newton?s Method [15] require that the Hessian (second-order derivative) be available. Other methods include Quasi-Newton methods, Levenberg-Marquardt (LM) [16], Se- quential Quadratic Programming (SQP) [15,17?21], etc. These algorithms often obtain the minimizing parameters with high accuracy and low computational cost provided that the gradient and Hessian are easily and accurately calculable. When these methods are applied to problems where the derivatives are not analytically available and must be cal- culated by numerical means or where the search space is discontinuous, their accuracy generally suffers or, in the worst case, may fail altogether. A similar problem occurs when the search space has many minima such as in the case of Rastrigin?s function, a 6 function often used to benchmark optimization algorithms (see Fig. 1.1) [22]. An ex- ample error surface for which local methods are particularly well-suited (surfaces having only one minimum) is shown in Fig. 1.2. X Y ?1 ?0.5 0 0.5 1?1 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 0.8 1 Figure 1.1: Rastrigrin?s function which has many local minima and a global minimum located at (0,0). Of the gradient-based methods, SQP is currently considered to be one of the most robust nonlinear programming (NP) methods available for optimization problems [20]. NP refers to the class of optimization problems in which the error functions are non- linear in their parameters and constraints. Although SQP is considered to be a local optimization method, it utilizes features which enhance its ability to avoid local min- ima. By taking a quadratic expansion of the nonlinear error function, the Hessian and Lagrangian operators can be efficiently approximated and used to create a Quadratic 7 X Y 0 0.2 0.4 0.6 0.8 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 1.2: Error surface well-suited to local optimization problems. 8 Programming subproblem. The solution of the subproblem results in the selection of the line search direction. Depending on the accuracy of the initial guess, SQP has been shown to exhibit second order convergence toward a local minimum. Therefore, when local methods are needed in the problems presented in this dissertation, SQP is utilized. For further discussion of SQP, the interested reader is referred to [15,17?19,21]. Global Optimization Techniques Global optimization refers to the task of finding the absolute best set of parameters, usually over a very large search space, which will minimize the error function. A large number of methods have been developed to handle these tasks which include techniques well-suited for combinatorial problems, discontinuous problems, constrained problems, etc. Over the years, the methods have been grouped according to their approach to the optimization. These groups include Branch and Bound Methods [23], Clustering Methods [24], Evolutionary Algorithms [25], Statistical Methods [26], and miscellaneous and hybrid techniques [27]. Of late, Evolutionary Algorithms including the Genetic Al- gorithm (GA) [28], Particle Swarm Optimization (PSO) [29], and Differential Evolution (DE) [30] have received considerable attention. These methods generally rely on stochas- tic procedures to prevent the local minimum trapping that may hinder local methods. In addition, they are particularly well-suited for problems having search spaces which are discontinuous, constrained, multi-dimensional, noisy, and multi-extremal. Unfor- tunately, this robustness results in a greater computational expense since many error function evaluations are typically required to reach a good solution. Also, the methods 9 are considered inexact because their convergence usually suffers once a solution is found very close to the global minimum. Of GA?s, PSO, and DE, GA?s are chronologically the first to have been devel- oped [31?33]. GA?s are based on the concepts of evolutionary biology in which a ?sur- vival of the fittest? paradigm is adopted. The GA begins with selection of an initial population of parameter values and selectively evolves the population toward the global minimum of the search space. This method has been extensively researched and applied to many problems found within the EM community since the error functions encoun- tered here typically possess many of the qualities of difficult NP problems previously mentioned in addition to the fact that the error functions are computationally expensive to evaluate [34]. In fact, GA?s have been applied successfully to nearly all areas of elec- tromagnetics including antenna design [35], phased arrays [36], and radar absorbers [37]. Thorough discussions regarding applications of GA?s to electrical engineering problems can be found in [34,38,39]. PSO is similar to the GA in that it is a population-based technique which utilizes information from the population to move toward a globally optimal solution. PSO, developed by Eberhart and Kennedy in 1995, is based on social theories related to bird flocking, insect swarms, and fish schools [29]. PSO begins with an initial population of ?particles? each of which has a particular location and velocity within the search space. Utilizing the information of their ?neighbors?, each particle can adjust its flight direction and speed until all particles arrive at the global minima of the search space. Advanced algorithms employ many features such as fully connected neighborhoods, particle inertia, etc. which can accelerate convergence for particular applications [40,41]. Recently, 10 PSO has been applied to several EM problems including antenna design [42] and array synthesis [43]. DE is also a population-based technique and was developed by Price and Storn to solve the Chebychev Polynomial fitting problem [30]. It relies on simple vector algebra operations on the population members to quickly generate new population members which are closer to the global minimum. This technique has been shown in certain instances to converge faster than other evolutionary algorithms [44] and is also well-suited for parallelization. Although DE has yet to receive much attention in EM optimization, some research has been conducted on antenna arrays [45] and inverse scattering [46]. Of the evolutionary algorithms discussed, the GA and PSO are utilized in this work with the GA receiving a more extensive consideration. Details of the specific implementations and modifications of the algorithms are given in later sections. 1.3.2 Constitutive Parameter Extraction Multilayer substrate materials are currently used for many practical applications that include Microwave Integrated Circuits (MIC), MMIC [47], radomes, spatial filters for antenna beam shaping [48], and Frequency Selective Surfaces (FSS) [49]. By choosing the appropriate thicknesses and material parameters for the layers, it is possible to synthesize composite structures with novel electromagnetic properties otherwise not found in a single material [50]. Recently, attention has focused on non- destructive methods of determining the constitutive parameters of each individual layer of the substrate [51,52]. 11 There are a large number of methods for determining the permittivity or perme- ability of a single homogeneous sample or the effective ?bulk? properties of a layered material. These methods include split-cylinder resonators, cavity resonators, TE10 split- post dielectric and magnetic resonators, whispering-gallery resonators, transmission-line and waveguide techniques, etc. [53?55]. For measurements of the complex permittivity and permeability over broad frequency bands (e.g., X or L bands), transmission line or waveguide techniques are generally preferred even though the achievable accuracy is reduced due to unavoidable measurement errors [50,56]. Some recent methods used to accurately estimate the complex permittivity of in- dividual layers of a multilayered or inhomogeneous structure are given by Sanadiki and Mostafavi [57], Zwick et al. [52], and Deshpande and Dudley [50]. Sanadiki and Mostafavi provide a method of solving the inverse scattering problem using a least squares error approach. This method is only tested against computer generated data and may be sensitive to errors associated with measurements. Zwick et al., utilize a GA to find the complex constitutive parameters for a multilayered sample by an evolutionary process. Their method requires measurements be obtained over a frequency range or as a func- tion of incidence angle for a given frequency and does not require phase information for the transmission or reflection coefficients. Deshpande and Dudley?s algorithm employs SQP and utilizes both magnitude and phase information of the measured S-parameters which, in turn, requires very accurate system calibration. To the author?s knowledge, no published work has address the issue of extracting both complex permittivity and permeability from individual layers of a multilayered structure. 12 1.3.3 Waveguide Filter Optimization In recent years, there has been a steady increase in the demand for fast and effi- cient passive waveguide filter design tools applicable to E-plane filters [58,59], H-plane filters [60], dual-band filters [61], etc. Thus, many researchers have begun exploiting optimization techniques to alleviate the burden of traditional analytical design meth- ods [62?67]. Although the majority of earlier work in filter design is focused on the de- velopment of longitudinal filters (i.e., structures which are oriented longitudinally within the waveguide) [58,59,61?67], of late, some work has been conducted in the design and optimization of transverse filters (i.e., structures which have their prominent dimension oriented parallel to the waveguide face) [68?70]. Transverse filters have the advantage of being lightweight, low profile, and easy to manufacture. Additionally, transverse filters do not require modification of the supporting waveguide structure as do some longitudi- nal filters [67]. Lockyer and Vardaxoglou [68] and Seager et al. [69] have investigated transverse filter performance using two-layer aperture arrays within the waveguide [68,69] and were able to realize narrow bandpass filters using structures ?1 mm thick. They have shown that, for narrowband applications, a transverse filter can be constructed quite easily. Alternatively, Monorchio et al. [70] developed an optimization technique to design a wideband transverse filter for square waveguides based on a GA and Method of Moments (MoM) solver [71]. Their approach optimized a FSS screen which transmitted K and Ka bands while reflecting X and Ku bands. Since transverse waveguide filters are basically bounded FSS, a brief history of FSS optimization is also warranted. This topic has received significant attention due the 13 complexity of the methods involved in creating these structures [71?78]. Manara et al. [71] utilized a GA optimization procedure to refine the geometry of a multiband, single layer FSS structure which could efficiently transmit L and S bands while reflecting the Ka band. Parker et al. [73] also employed a GA to reduce the angular dependence of the reflection band for a FSS. Bozzi et al. [75] and Li et al. [78] applied fast solvers to the optimization problem thereby making the analysis of realistically complex structures feasible. Chakravarty and Mittra [74,76,77] demonstrated the Micro-GA?s ability to optimize multilayer spatial filters efficiently when several design parameters were allowed to vary over realistic ranges. Finally, Ohira et al. [72] utilized the GA to optimize mask shape for FSS elements using printed circuit technology. 1.4 Thesis Outline The remainder of this dissertation is organized in the following manner. Chap- ters 2 & 3 present the formulation, optimization algorithm, and results for the complex constitutive parameter extraction technique. Chapter 4 presents the numerical technique and results pertaining to the embedded conductor filters. In Chapter 5, modifications of the optimization algorithms presented thus far are discussed as well as results for the waveguide filters generated using surface patches. Finally, Chapter 6 presents conclusions and suggestions for future work. 14 Chapter 2 Complex Permittivity Extraction 2.1 Overview In this chapter, the performance of several optimization techniques is evaluated when applied to the multilayer complex permittivity extraction problem. S-parameter measurements obtained from a loaded waveguide serve as an input to the inversion algo- rithms. The methods to be considered are Sequential Quadratic Programming (SQP), Genetic Algorithm (GA), and Particle Swarm Optimization (PSO). Their sensitivity and performance relative to the choice of error function is also investigated. Simulations are performed using computer generated S-parameters to quantify the performance of each algorithm under ideal conditions. In order to determine the robustness of these methods, the extracted permittivities determined by the algorithms are used to generate S-parameter data sets for comparison to the measured S-parameters. S-parameter X band waveguide measurements, provided by NASA Langley Re- search Center [50], were made using the following dielectric materials: Bakelite, Ceramic, Garlock-Rubber, and Nano Material. All of the materials are low-loss, non-magnetic, and were found to have complex permittivities that remained nearly constant over the X band [50]. They were used to create planar, single layer samples of various thicknesses, with all samples having X band waveguide cross-sectional dimensions. The multilayer di- electric structures were fabricated by placing single layer samples adjacent to one another (see Fig. 2.1). All measurements were obtained using an HP-8510C Vector Network An- alyzer for frequencies of 8.2 - 12.4 GHz (X band). Additionally, the complex permittivity 15 of each single layer sample was determined using the Agilent 85071 Materials Measure- ment Software [79,80] and compared to the values obtained from the algorithms. 2.2 Theory Assuming that only the dominate TE10 mode propagates in the loaded waveguide, as shown in Fig. 2.1, the formulation of the S-parameters can be expressed in terms of each layer?s thickness (dn) and unknown permittivity (??n) using ABCD-parameters as shown below [81,82]: ? ??An Bn Cn Dn ? ??= ? ?? cosh(?ndn) Zn sinh(?ndn) sinh(?ndn)/Zn cosh(?ndn) ? ?? (2.1) ? ??A B C D ? ??= nproductdisplay i=1 ? ??Ai Bi Ci Di ? ?? (2.2) where ?n and Zn are the propagation constant and wave impedance of the nth layer, respectively, and can be expressed as Zn = j??? n (2.3) and ?n = radicalBig (?/a)2 ??2???c (2.4) where a is the longer dimension of the waveguide. 16 Figure 2.1: Rectangular waveguide loaded with n-layer sample. 17 The ABCD-parameters can be directly converted to S-parameters using the following equations S11 = (A+B/Z0 ?CZ0 ?D)/X S12 = 2(AD?CB)/X S21 = 2/X S22 = (?A+B/Z0 ?CZ0 +D)/X X = A+B/Z0 +CZ0 +D (2.5) where Z0 is the empty waveguide impedance. 2.3 Error Function Selection The error functions used in the investigation are of the basic form given by Err = radicaltpradicalvertex radicalvertexradicalbt 1 N Nsummationdisplay i=1 bracketleftBig f parenleftBig [S]f,i parenrightBig ?f parenleftBig [S]m,i parenrightBigbracketrightBig2 (2.6) whereN isthetotalnumberoffrequencypoints, [S]f,i aretheformulatedS-parameters(2.5) evaluated at frequency point i, and [S]m,i are the measured S-parameters at frequency point i. Error functions can generally be grouped into two categories: (1) error func- tions used to minimize differences in both phase and magnitude of the measured and formulated scattering parameters, and (2) error functions that minimize only magnitude variations of the measured and formulated scattering parameters. For the purposes of this study, three error function definitions were used and found to accurately obtain the permittivity of single and multilayer structures: 18 Err1 = radicaltpradicalvertex radicalvertexradicalvertex radicalvertexradicalvertex radicalbt 1 N Nsummationdisplay i=1 ? ?? ? parenleftBig Re parenleftBig [S]f,i parenrightBig ?Re parenleftBig [S]m,i parenrightBigparenrightBig2 + parenleftBig Im parenleftBig [S]f,i parenrightBig ?Im parenleftBig [S]m,i parenrightBigparenrightBig2 ? ?? ? (2.7) Err2 = radicaltpradicalvertex radicalvertexradicalbt 1 N Nsummationdisplay i=1 parenleftBigvextendsinglevextendsingle vextendsingle[S]f,i vextendsinglevextendsingle vextendsingle? vextendsinglevextendsingle vextendsingle[S]m,i vextendsinglevextendsingle vextendsingle parenrightBig2 (2.8) Err3 = radicaltpradicalvertex radicalvertexradicalbt 1 N Nsummationdisplay i=1 Mi +P1,i +P2,i (2.9) where Mi = parenleftBigvextendsinglevextendsingle vextendsingle[S]f,i vextendsinglevextendsingle vextendsingle? vextendsinglevextendsingle vextendsingle[S]m,i vextendsinglevextendsingle vextendsingle parenrightBig2 Pf1,i = 1?|S11|2f,i ?|S21|2f,i Pf2,i = 1?|S22|2f,i ?|S12|2f,i Pm1,i = 1?|S11|2m,i ?|S21|2m,i Pm2,i = 1?|S22|2m,i ?|S12|2m,i P1,i = (Pf1,i ?Pm1,i)2 P2,i = (Pf2,i ?Pm2,i)2 Err1, given in (2.7), includes both phase and magnitude information. It is a slightly modified form of Deshpande and Dudley?s error function [50] and is applicable to materi- als where the permittivity remains approximately constant over the measured frequency range. The inclusion of measurement information over the entire frequency range tends to minimize the effects of instrumentation error in the calculations. Err2, given in (2.8), requires only magnitude information of the scattering parameters and is representative 19 of the fitness (error) function used by Zwick et al. [52] and Queffelec and Gelin [54] . The third error function Err3, given in (2.9), is unique in that it uses only magnitude information but includes terms accounting for dissipated power. It was initially believed that the inclusion of the power terms would increase the algorithms? ability to accurately determine the imaginary part of the permittivity. It should be noted that Eqns. (2.8) and (2.9) are calculated in decibels which was found to decrease the time required to find the global minimum in the solution space [52]. 2.4 Optimization Algorithms 2.4.1 Sequential Quadratic Programming Since publication of a series of expository papers in the 1970?s, SQP has become the most popular method for solving Nonlinear Programming (NP) problems [83?85]. SQP, itself, is not a specific NP algorithm but a framework under which many implementations have evolved. A thorough treatment of the mathematical details of SQP and its sub- algorithms is beyond the scope of this dissertation, and interested readers are referred to the foundational papers on SQP [18,86,87]. Instead, a brief and general overview is favored since specific implementations are readily available and deserve individual treatment [88]. As mentioned in Section 1.3.1, the goal of NP is to solve the problem f (vectorx) = 0 (2.10) subject to the constraints 20 ci (vectorx) = 0 gj (vectorx) ? 0 (2.11) where vectorx is the vector of optimization parameters, ci (vectorx) are the set of equality constraints on vectorx, and gj (vectorx) ? 0 are the set of inequality constraints on vectorx. SQP?s novelty lies in its ability to solve problems efficiently which are nonlinearly constrained. This is possible through the quadratic approximation of the Lagrangian function. The Langrangian function, in its unmodified form, is given by L = f(vectorx)+vectoruTvectorc(vectorx)+vectorvTvectorg(vectorx) (2.12) where vectoru and vectorv are the Lagrangian multipliers. The reasoning behind forming the quadratic subproblem is that a number of techniques are available for solving error functions which are quadratic and subject to linear constraints. Therefore, the approx- imation is formed as a quadratic function around the current iterate?s location, and the constraints are linearized. Once the quadratic approximation is formed, this can be solved by a number of Quadratic Programming (QP) methods [15,89,90] to form a line search direction. This eventually leads to the selection of the next evaluation point which will converge to the minimum of f. The basic SQP outline is given as follows: 1. Obtain initial parameters (i.e., initial guess). 2. Solve the QP problem at this point to obtain the line search direction. 3. Choose a step length such that the error is reduced sufficiently. 21 4. Update the current iterate and Lagrangian multipliers for the next iteration. 5. Check stopping criteria (convergence). 6. Proceed to next iteration. 2.4.2 Genetic Algorithm The GA, the foundations of which are found in evolutionary biology [25,31?33], is based on the concept of evolution. Therefore, much of the terminology used to explain the algorithm is biological in nature. GA?s, like the other global optimization algo- rithms used in this work, are somewhat immune to local minimum trapping since they employ stochastic processes. GA?s find application in a wide variety of fields such as bio- science [91], neural networks [92], economics [93], and many more [94]. The remainder of the section outlines the GA implemented in this work. The GA begins with the creation of a population of parameter values randomly dis- tributed throughout the solution space (see Fig. 2.2). Each member of the population is referred to as a ?chromosome? and contains one ?gene? storing a value for each parame- ter in the function to be minimized. In traditional GA?s, each gene is encoded as an n-bit binary string which is used for creation of the next population. Upon creation of the initial population, the so-called fitness of each chromosome is evaluated by determining the error function value associated with that chromosome?s list of parameters. After fitness evaluation, a new population (generation) is created. Construction of the new population begins by selecting a subset population representing the ?best? chromosomes (i.e., those chromosomes having the lowest error function value). These 22 X59X65X73 X4EX6F X4EX6F X59X65X73 X44X65X74X65X72X6DX69X6EX65X20 X45X6CX69X74X65X20 X43X68X69X6CX64X72X65X6E X49X6EX69X74X69X61X6CX69X7AX65X20 X50X6FX70X75X6CX61X74X69X6FX6E X45X76X61X6CX75X61X74X65X20 X46X69X74X6EX65X73X73 X53X65X6CX65X63X74X69X6FX6E X43X72X6FX73X73X6FX76X65X72 X4DX75X74X61X74X69X6FX6E X50X6FX70X75X6CX61X74X69X6FX6EX20 X52X65X64X75X6EX64X61X6EX63X69X65X73X3F X52X65X6DX6FX76X65X20 X43X6FX70X69X65X73 X52X65X69X6EX69X74X69X61X6CX69X7AX65X20 X54X6FX20X46X69X6CX6CX20 X50X6FX70X75X6CX61X74X69X6FX6E X45X78X69X74X20X43X6FX6EX64X69X74X69X6FX6EX3F X45X78X69X74 Figure 2.2: Flowchart illustrating the Genetic Algorithm procedure. 23 ?elite? chromosomes are inserted into the new population ensuring that these solutions will not be eliminated from the new generation. The crossover stage is the defining phase of the GA. First, two ?parent? chromo- somes are selected from the current generation using the binary tournament selection method [38]. In binary tournament selection, two chromosomes are picked at random. The chromosome with the best fitness is chosen as the first parent. This process is re- peated for the second parent. Following selection of the parent chromosomes, the GA generates a random number p ? [0,1]. If p > pcrossover, the crossover rate, the two parents are copied directly into the new population. Otherwise, a random position in the parents? binary string of encoded information is selected (see Fig. 2.3). The first ?child? is formed using the string of information to the left of the crossover point in the first parent and the information to the right of the crossover point in the second parent. The second child is formed using the remaining genetic information from each parent. This process of selection and crossover is repeated until the new population is completed. Using this concept, chromosomes with poor fitness that possibly possess ?good? genetic information have the opportunity to pass their traits on to future generations. The next major phase of the GA is the mutation stage. For each child in the new population (elite children excluded), a random number m ? [0,1] is generated. If m > pmutation, the mutation rate, no change is made to the chromosome. However, if m < pmutation, one or two random bits of the chromosome are transposed (i.e., 1 ? 0). This reduces the likelihood of local minima trapping. The GA will discontinue if a chromosome is found to have a fitness value below 10?5 or if the maximum number of generations has been reached. Depending on the number 24 1 0 1 0 1 0 1 0 1 0 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 Crossover Point 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1 Parent 1 Parent 2 Child 1 Child 2 Figure 2.3: Crossover between two parents using binary string encoding. 25 of layers being characterized, the maximum number of generations is allowed to vary. The specific details of this parameter?s selection will be discussed in Section 2.5. GA Redundancy Removal At each generation, during the creation of the new population, some chromosomes may be duplicated. Several novel redundancy removal schemes were developed and their effect on convergence rate was analyzed. The first scheme developed is referred to as the Complete Elite Redundancy Removal GA (CERRGA). After a significant number of generations, the GA converges to a population consisting of identical chromosomes representing the best obtained solution to the fitness function (with the exception of the occasional mutation). At this point, a traditional GA exits since simple mutation alone is a very inefficient way of searching for new genetic information that will more efficiently minimize the error function. In the CERRGA method, the GA checks for this redundancy in the population for each new generation. When this situation occurs, one copy of this redundant chromosome is kept as well as any mutations. The rest of the population slots are filled with random chromosomes in the same manner in which the initial population was generated, and the algorithm is restarted at the current iteration. The next redundancy removal scheme employed was termed Incomplete Elite Re- dundancy Removal GA (IERRGA). This method is very similar to CERRGA in that elite chromosome redundancy is eliminated. However, IERRGA allows a factor to be set that controls how much of the population is filled with elite chromosomes before the re- dundancies are eliminated. In the results presented shortly, a 90% IERRGA scheme was used. Therefore, the algorithm monitored the number of elite chromosome redundancies 26 present in the population. If that number exceeded 90%, only one copy was kept and the remaining population slots were filled with random chromosomes. This differs from CERRGA in that the new random chromosomes have a greater chance of reproducing with one of the elite copies. Therefore, in theory, there may be a greater number of ?good? chromosomes in the population at any one time than with CERRGA. The last redundancy removal method monitors the population at each generation and removes all redundancies. Therefore, only one copy of any chromosome is present in the population. All emptied population slots are, as always, filled with random chro- mosomes. This scheme is termed the Total Redundancy Removal GA (TRRGA). These methods were tested using a computer generated S-parameter data set for one layer (see Section 2.5). Fig. 2.4 shows the average fitness value of the best chromosome at each generation from 20 simulations for all of the methods presented as well as the unaltered GA (i.e., no redundancy removal). From the figure, it is obvious that the TRRGA method provides significantly faster convergence than the traditional GA and CERRGA and is better than IERRGA by approximately 80 generations. Therefore, the results presented in the remaining sections correspond to those obtained by the TRRGA. 2.4.3 Particle Swarm Optimization PSO was developed by Eberhart and Kennedy [29] as an optimization procedure based on models of social behaviors such as bird flocking. In many ways, it bears a strong resemblance to the GA except that there are no genetic operators such as crossover and mutation. PSO handles the generation of a ?new population? in a much different way as will be explained. Similar to the GA, PSO has received much attention in a variety of 27 0 50 100 150 200 10?6 10?4 10?2 100 Generation Fitness Traditional GA All redundancies Elite redundancies 90% Elite redundancies PSO Figure 2.4: Comparison of convergence for the traditional GA and GA?s using three redundancy removal methods. The convergence rate of PSO is also shown. 28 fields including EM as previously mentioned. For a thorough review of previous research regarding PSO, the interested reader is referred to Hu [95]. The PSO algorithm utilizes a population of ?particles? which are given a position and velocity within the search space. Each of these particles, depending on the complex- ity of the algorithm, has a set of neighbors, so-called individuality and sociality factors which determine how much they are influenced by their neighbors, and an inertia factor which includes the effects of previous iterations. Each particle usually is also given a ?memory? of the best solution it has encountered which can also be made to affect the particle?s propagation through the solution space. The following paragraphs detail the PSO algorithm implemented in this work. As with most Evolutionary Algorithms, PSO begins by initializing a population of particles with a random location and velocity (see Fig. 2.6). In this implementation, individuality and sociality factors and neighbors are also initialized. A common choice of individuality and sociality factors is 2 [29]. A five-particle neighborhood is established so that each particle is influenced by the two particles to its left and right (i.e., particle 3 has neighbors 1, 2, 4, and 5 while particle 4 has neighbors 2, 3, 5, and 6) as shown in Fig. 2.5. The inertia factor is initialized to 1 such that the particle is also influenced heavily by the previous iteration?s information. Each particle?s fitness is then tested by evaluating the fitness function (forward solution). If this fitness is better than the best previous fitness of that particle, the current location is stored as the best. Also, if any fitness is encountered which is below a fitness threshold, then the algorithm ends because a sufficient solution has been found. 29 X31 X32 X33 X34 X35 X36 X37 X38 X39X31X30X31X31 X31X32 X31X33 X31X34 X31X35 X31X36 X31X37 X31X38 Figure 2.5: Illustration of the neighborhood structure implemented in PSO. Each particle is connected to its two neighbors to the left and right. 30 After every particle?s fitness has been evaluated, the velocity and, hence, next lo- cation can be calculated. Below, a pseudocode is presented to explain this calculation. This process is repeated for each particle. %Return the neighbor who has had the best fitness n = GetNeighborWithBestFitnessSoFar() %Return the best location of that neighbor n_Best = GetBestLocationOfParticle(n) %Multiply the sociality factor by a random number [0,1] Reduced_S = S_Factor*rand1 %Multiply reduced sociality factor by difference between %neighbor?s best location and particle?s current location n_Change = Reduced_S*(n_Best - p_Current) %Multiply individuality factor by a random number [0,1] Reduced_I = I_Factor*rand2 %Multiply reduced individuality factor by difference between %particle?s best location and current location p_Change = Reduced_I*(p_Best - p_Current) %Get change in location for the particle by summing up all %change contributions including the inertia factor multiplied %by the current velocity Change = inertia*p_Velocity + n_Change + p_Change %Check to see if the change is greater than what is allowed. %This helps limit changes that are too small or too large %so that the particles better cover the search space p_Velocity = CheckChange(Change, Max_Change, Min_Change) %Update next location by adding current location to the %calculated change in location p_Next = p_Current + p_Velocity Once this process is complete, the particles? locations are updated for the next iteration, and the procedure of fitness evaluation begins again. At each iteration, the inertia factor is reduced by a uniform amount until it reaches 0.2. By reducing the inertia at each iteration, the particles are slowly forced to converge. It was found that this method resulted in a faster convergence rate than simply leaving the factor set to 1. The PSO algorithm presented here was tested using the same one layer computer generated data as presented in the previous section. Again, 20 runs were performed 31 and the average convergence rate calculated. Fig. 2.4 shows the convergence curve for the PSO as compared to the GA?s presented. While the PSO returns equally accurate answers as compared to the GA, it requires many more iterations to do so than the TRRGA method. Therefore, a complete analysis of the cases considered in the next sections was not conducted using PSO. For brevity, the cases that were analyzed using PSO are not presented since their results are the same as what was found using the GA. 2.5 Numerical Results Computer generated S-parameter data sets for 1-, 2-, and 3-layer cases over X band frequencies were initially used to test the accuracy of each extraction scheme. The respective complex permittivities and thicknesses are shown in Table 2.1 (single layer permittivity extraction used layer 1; 2-layer extraction used layers 1 and 2; 3-layer extraction used all layers ordered accordingly). The initial complex permittivity guess for all layers was set to ??initial = 5?j0.4, while the search space for both algorithms was limited to 1 < Re(??c) < 10 and 0 < Im(??c) < 0.8. Table 2.1: Complex Permittivities and Thicknesses for the 1-, 2-, and 3-Layer Computer Generated Data Sets Layer ??n dn (mm) 1 7?j0.01 1 2 3?j0.02 10 3 2?j0.1 2 The SQP algorithm terminated for error function or directional derivative values ? 10?16 or after 5000 iterations. The GA?s population size was set to a value of 100 and terminated after 200, 1000, and 5000 generations for the 1-, 2-, and 3-layer cases, 32 X49X6EX69X74X69X61X6CX69X7AX65X20 X50X6FX70X75X6CX61X74X69X6FX6E X45X76X61X6CX75X61X74X65X20 X46X69X74X6EX65X73X73 X43X61X6CX63X75X6CX61X74X65X20X4EX65X77X20 X4CX6FX63X61X74X69X6FX6EX73 X55X70X64X61X74X65X20 X50X61X72X74X69X63X6CX65X73 X52X65X64X75X63X65X20 X49X6EX65X72X74X69X61 X45X78X69X74X20 X43X6FX6EX64X69X74X69X6FX6EX3F X45X78X69X74 X4EX6F X59X65X73 X49X6EX64X69X76X69X64X75X61X6CX69X74X79 X53X6FX63X69X61X6CX69X74X79 X4EX65X69X67X68X62X6FX72X73 X49X6EX65X72X74X69X61 Figure 2.6: Flowchart illustrating the PSO procedure. 33 respectively. The GA utilized an 80% crossover rate and 10% mutation rate (these values are within the range of nominal values given by [38]). In Table 2.2, the results for the 1-, 2-, and 3-layer cases using each algorithm (single or multiobjective) and error function are shown. An ?X? in the table indicates that the algorithm extracted incorrect complex permittivity values for the indicated sample (i.e., the algorithm became trapped in a local minimum). For the multiobjective cases, each term in the error functions (2.7)-(2.9) becomes an element of an error function vector. For instance, in a multiobjective format, Eqn. (2.8) would then be given by ???Err = ? ?? ?? ?? ?? ?? ?? ?? ? radicalBigg 1 N Nsummationtext i=1 (|S11,f,i|?|S11,m,i|)2 radicalBigg 1 N Nsummationtext i=1 (|S21,f,i|?|S21,m,i|)2 radicalBigg 1 N Nsummationtext i=1 (|S12,f,i|?|S12,m,i|)2 radicalBigg 1 N Nsummationtext i=1 (|S22,f,i|?|S22,m,i|)2 ? ?? ?? ?? ?? ?? ?? ?? ? . (2.13) The vectorization of the error functions was expected to increase the sensitivity of the SQP algorithm. The results shown in the table indicate that all methods performed extremely well for the single and double layer complex permittivity extraction. However, permittivity extraction for the 3-layer sample was unsuccessful using SQP with error functions (2.8) and (2.9). This is due to local minimum trapping indicating a poor initial estimate (guess) of the complex permittivity. Alternatively, the GA was able to successfully extract complex permittivity values using all three error functions independent of the number of layers. It should also be noted that the error for both the real and imaginary 34 Table 2.2: Results using ideal S-parameters Error One Layer Two Layers Three LayersAlgorithm Objective Function Error Time (s) Error Time (s) Error Time (s) (2.7) 1.837e-9 59.1 1.729e-8 4.3 2.431e-8 6.8 Single (2.8) 2.702e-8 1.9 1.780e-7 65.6 X X (2.9) 1.959e-11 50.1 1.396e-7 65.6 X X (2.7) 1.712e-13 2.1 2.012e-9 72.7 6.245e-9 83.3 Multiple (2.8) 6.397e-9 54.4 8.538e-8 4.4 X X SQP (2.9) 2.183e-9 54.2 1.656e-7 3.1 X X (2.7) 3.561e-7 68.8 2.636e-4 969 2.079e-4 6467.2 GA Single (2.8) 3.106e-6 103.1 2.003e-3 1136.4 1.151e-2 9683.3 (2.9) 3.147e-6 124.7 7.266e-3 2059 1.058e-2 11715 35 parts of the complex permittivity is less than O(10?7) for the SQP using (2.7) and less than O(10?3) for the GA using any of (2.7)-(2.9) for the 1-, 2-, and 3-layer cases. 2.6 Measurements Results 2.6.1 Single Layer Measurements For measurements of single layer materials, samples of Bakelite (d = 3.277 mm), Ce- ramic (d = 2.845 mm), Garlock Rubber (d = 1.702 mm), and Nano Material (d = 3.099 mm) were used. S-parameter data sets were generated using the extracted permittivity value from each algorithm and then compared to the measured S-parameter values. Ad- ditionally, the Agilent 85071 Materials Measurement Software was used to extract the single layer permittivities for comparison to the values returned by SQP and the GA. Fig. 2.7 shows a comparison of |S11| and |S21| for the measured S-parameters of the Bakelite sample and those generated from the extracted permittivities returned by the GA using error function (2.9). Fig. 2.8 shows similar results for the Garlock sample. The excellent agreement shown in both figures is also observed for all other materials? S-parameter comparisons using all algorithm/error function combinations. 36 8.5 9 9.5 10 10.5 11 11.5 120.5 0.55 0.6 0.65 0.7 0.75 0.8 Frequency (GHz) |S| |S11,f| |S11,m| |S21,m| |S 21,f| Figure 2.7: Magnitude of S11 and S21 for single layer Bakelite sample. The calculated S-parameters are generated using the GA with error function (2.9). 37 8.5 9 9.5 10 10.5 11 11.5 120.5 0.6 0.7 0.8 0.9 1 Frequency (GHz) |S| |S11,f| |S11,m| |S21,m| |S21,f| Figure 2.8: Magnitude of S11 and S21 for single layer Garlock Rubber sample. The calculated S-parameters are generated using the GA with error function (2.9). 38 Tables 2.3-2.6 show the results of the extracted complex permittivities for each material using all optimization algorithms and error functions. Since error functions (2.8) and (2.9) are calculated in decibels, a direct comparison of the error function values is not necessarily indicative of the accuracy of the solution. Rather, summations of the RMS errors between the magnitudes and phases of the formulated and measured S-parameters (see Columns 6-7 in Tables 2.3-2.6) are better qualifiers of the accuracy of each solution and expressed by summationtext ij radicalBig |Sij,f|2 ?|Sij,m|2 summationtext ij radicalBig (negationslash Sij,f)2 ?(negationslash Sij,m)2 . (2.14) The RMS errors between the measured S-parameters and those generated using the Agilent 85071 Materials Measurement Software are also listed for comparison. It can be concluded that all of the optimization techniques were effective in minimizing the RMS errors of the S-parameters (both phase and magnitude). However, the GA using either (2.8) or (2.9) was found to consistently produce the lowest magnitude error of the S-parameters (column VI) and maintained a phase error generally no worse than the other extraction methods. 39 Table 2.3: Extracted Permittivity & Error for Bakelite ErrorAlgorithm Objective Function ? ? ??? summationtext|Sij| err summationtextnegationslash S ijerr (2.7) 3.7244 0.2237 0.02875 0.19899 Single (2.8) 3.7907 0.2528 0.01339 0.20948 (2.9) 3.7909 0.2530 0.01339 0.20952 (2.7) 3.5124 0.3464 0.09824 0.29814 Multiple (2.8) 3.7896 0.2656 0.01434 0.20685 SQP (2.9) 3.8018 0.2577 0.01364 0.21206 (2.7) 3.7244 0.2238 0.02876 0.19898 GA Single (2.8) 3.7907 0.2528 0.01339 0.20948 (2.9) 3.7910 0.2530 0.01339 0.20954 85071 X X 3.6032 0.2347 0.06947 0.24004 40 Table 2.4: Extracted Permittivity & Error for Ceramic ErrorAlgorithm Objective Function ? ? ??? summationtext|Sij| err summationtextnegationslash S ijerr (2.7) 1.1422 0.0000 0.03263 0.15330 Single (2.8) 1.1819 0.0004 0.00415 0.19111 (2.9) 1.1819 0.0006 0.00415 0.19151 (2.7) 1.1344 0.0000 0.03888 0.14662 Multiple (2.8) 1.1782 0.0518 0.04234 0.70040 SQP (2.9) 1.1823 0.0050 0.00670 0.21894 (2.7) 1.1423 0.0000 0.03261 0.15333 GA Single (2.8) 1.1818 0.0001 0.00417 0.18989 (2.9) 1.1820 0.0010 0.00422 0.19331 85071 X X 1.1484 0.0043 0.02802 0.18255 41 Table 2.5: Extracted Permittivity & Error for Garlock Rubber ErrorAlgorithm Objective Function ? ? ??? summationtext|Sij| err summationtextnegationslash S ijerr (2.7) 7.7850 0.0409 0.04200 0.13297 Single (2.8) 7.9837 0.1642 0.01669 0.17187 (2.9) 8.1358 0.1304 0.02075 0.20980 (2.7) 7.6184 0.0869 0.06706 0.08846 Multiple (2.8) 8.0340 0.1604 0.01425 0.18396 SQP (2.9) 8.0347 0.1602 0.01424 0.18415 (2.7) 7.6184 0.0869 0.06707 0.08846 GA Single (2.8) 8.0340 0.1604 0.01425 0.18396 (2.9) 8.0348 0.1602 0.01424 0.18417 85071 X X 7.5320 0.1610 0.08267 0.06454 42 Table 2.6: Extracted Permittivity & Error for Nano Material ErrorAlgorithm Objective Function ? ? ??? summationtext|Sij| err summationtextnegationslash S ijerr (2.7) 2.5165 0.0000 0.05783 0.08785 Single (2.8) 2.5989 0.0119 0.01058 0.18065 (2.9) 2.5989 0.0119 0.01058 0.18061 (2.7) 2.4916 0.0000 0.07457 0.06210 Multiple (2.8) 2.5842 0.0000 0.01671 0.15915 SQP (2.9) 2.5837 0.0064 0.01530 0.16201 (2.7) 2.5165 0.0000 0.05779 0.08791 GA Single (2.8) 2.5989 0.0119 0.01057 0.18069 (2.9) 2.5989 0.0119 0.01057 0.18069 85071 X X 2.6872 0.0260 0.05989 0.27670 43 2.6.2 Two Layer Measurements The 2-layer structures used for the S-parameter measurements were Garlock Rub- ber/Bakelite, Garlock Rubber/Ceramic, and Garlock Rubber/Nano Material with Gar- lock Rubber used as the first layer for all cases. As discussed in Section 2.6.1, the three error functions (2.7)- (2.9) were minimized using SQP and the GA. No change was made to any paramater of the SQP algorithim, whereas the GA was allowed 1000 generations before termination. The extracted permittivities were used to generate S-parameters data sets for RMS error computation and were also contrasted to the Agilent 85071 re- sults of the previous section. It should be mentioned that the Agilent 85071 software can only determine the permittivity of single materials or provide bulk permittivity estimates for composite structures. Tables 2.7-2.9 show the extracted permittivities and S-parameter RMS errors for the three 2-layer samples. Using 2-layer S-parameter data sets generated from the single layer Agilent 85071 extracted permittivities, the RMS errors between these data sets and the measured S-parameters were calculated and are also included in the tables. For the Garlock Rubber/Bakelite sample, all algorithm/error function combinations achieved RMS magnitude errors O(10?2) and phase errors O(10?2 ? 10?1). The magnitude of the extracted permittivities for each layer showed excellent agreement with the extracted single layer values. However, the imaginary part of the extracted permittivity for Garlock Rubber showed an increased value when compared to the values given in Table 2.5. For the Garlock Rubber/Ceramic and Garlock Rubber/Nano Material samples (see Tables 2.8 and 2.9), the SQP algorithm was only effective at estimating the complex permittivity values using error function (2.7). Using error functions (2.8) and (2.9), 44 Table 2.7: Extracted Permittivity & Error for Garlock/Bakelite Error Garlock Garlock Layer 2 Layer 2Algorithm Objective Function ?? ??? ?? ??? summationtext|S ij|err summationtextnegationslash S ijerr (2.7) 7.9050 0.2207 4.0271 0.3274 0.0400 0.0647 Single (2.8) 7.5613 0.3626 3.9900 0.2563 0.0366 0.1997 (2.9) 7.8488 0.3343 3.9881 0.2583 0.0387 0.1036 (2.7) 7.8211 0.3080 4.0086 0.2908 0.0379 0.0991 Multiple (2.8) 7.9793 0.2813 3.9827 0.2980 0.0410 0.0699 SQP (2.9) 7.9793 0.2813 3.9827 0.2980 0.0410 0.0699 (2.7) 7.9050 0.2207 4.0271 0.3274 0.0400 0.0647 GA Single (2.8) 7.5620 0.3625 3.9900 0.2563 0.0366 0.1994 (2.9) 7.8523 0.3336 3.9883 0.2585 0.0387 0.1027 85071 X X 7.5320 0.1610 3.6032 0.2347 0.1136 0.3289 45 Table 2.8: Extracted Permittivity & Error for Garlock/Ceramic Error Garlock Garlock Layer 2 Layer 2Algorithm Objective Function ?? ??? ?? ??? summationtext|S ij|err summationtextnegationslash S ijerr (2.7) 7.9478 0.1169 1.2254 0.0373 0.0268 0.0601 Single (2.8) 1.0000 0.0712 6.1163 0.0963 0.0687 6.3055 (2.9) 1.0000 0.1167 6.1037 0.0656 0.0695 6.3189 (2.7) 7.8983 0.0405 1.2596 0.0649 0.0378 0.0686 Multiple (2.8) 1.0000 0.2616 6.1208 0.0066 0.0732 6.3536 SQP (2.9) 1.0000 0.2616 6.1208 0.0066 0.0732 6.3536 (2.7) 7.9479 0.1169 1.2254 0.0373 0.0268 0.0601 GA Single (2.8) 8.0688 0.1916 1.1930 0.0000 0.0144 0.0639 (2.9) 8.0866 0.1927 1.2167 0.0000 0.0146 0.0654 85071 X X 7.5320 0.1610 1.1484 0.0043 0.0803 0.1267 46 Table 2.9: Extracted Permittivity & Error for Garlock/Nano Material Error Garlock Garlock Layer 2 Layer 2Algorithm Objective Function ?? ??? ?? ??? summationtext|S ij|err summationtextnegationslash S ijerr (2.7) 7.6909 0.1382 2.7840 0.0585 0.0307 0.0490 Single (2.8) 2.5193 0.1860 5.4421 0.0110 0.0178 6.5273 (2.9) 2.5193 0.1978 5.4555 0.0020 0.0178 6.5223 (2.7) 7.6700 0.1435 2.7751 0.0380 0.0308 0.0552 Multiple (2.8) 2.4387 0.1856 5.3484 0.0152 0.0217 6.5841 SQP (2.9) 2.4729 0.1596 5.4659 0.0147 0.0204 6.5315 (2.7) 7.6909 0.1382 2.7841 0.0585 0.0307 0.0490 GA Single (2.8) 7.7636 0.2344 2.6881 0.0142 0.0132 0.0942 (2.9) 7.7544 0.2510 2.6866 0.0068 0.0136 0.0972 85071 X X 7.5320 0.1610 2.6872 0.0260 0.0294 0.1103 47 SQP was unable to correctly estimate each layer?s complex permittivity as evidenced by the unacceptable large phase errors. It is apparent from the discrepancies in the estimated permittivity values and the large S-parameter phase errors that this was a result of local minimum trapping. Unlike SQP, the GA was successful in extracting the complex permittivities using all error functions. The real and imaginary parts of the permittivities showed excellent agreement with the extracted single layer values, and the GA achieved RMS magnitude errors O(10?2) and phase errors O(10?2 ? 10?1). 2.6.3 Three Layer Measurements The S-parameters for a 3-layer composite structure (Nano Material/Garlock Rub- ber/Garlock Rubber) were measured. In accordance with the previous sections, Ta- ble 2.10 shows the extracted permittivities and S-parameter RMS errors for the compos- ite structure. SQP was effective in determining the complex permittivites of the sample for all but the multiobjective form of (2.7). As in the previous cases, the GA succesfully estimated the permittivities using all error functions. SQP and the GA returned RMS magnitude errors O(10?2) and phase errors O(10?2 ? 10?1), respectively. It was already determined that the algorithms and error functions would return complex permittivity values nearly indistinguishable from the actual values for all com- puter generated cases (SQP using (2.8) and (2.9) excepted). Therefore, the possible reasons for the errors sited previously warrant an explanation. The most likely sources of error for the single layer extractions stem from inaccuracies associated with instru- mentation and sample thickness measurements. In addition to these errors, compound structures may suffer from the presence of small air gaps between the layers as well as 48 Table 2.10: Extracted Permittivity & Error for Nano Material/Garlock/Garlock Error Layer 1 Layer 1 Layer 2 Layer 2 Layer 3 Layer 3Algorithm Objective Function ?? ??? ?? ??? ?? ??? summationtext|S ij|err summationtextnegationslash S ijerr (2.7) 2.8234 0.0615 8.3424 0.2413 8.1992 0.2086 .03571 .06789 Single (2.8) 2.9228 0.0064 8.2085 0.3069 7.3469 0.3813 .01132 .28668 (2.9) 2.9255 0.0133 8.1900 0.3014 7.3365 0.3561 .01178 .28771 (2.7) 2.7634 0.0371 8.3823 0.3812 8.3505 0.3224 .05696 .06735 Multiple (2.8) 4.1458 0.0082 4.7201 0.2681 2.2123 0.8000 .23606 3.6167 SQP (2.9) 4.6797 0.0000 2.4478 0.6366 2.8608 0.0000 .14824 4.5302 (2.7) 2.8238 0.0625 8.3417 0.2402 8.1981 0.2051 .03568 .06775 GA Single (2.8) 2.9159 0.0125 8.1326 0.3000 7.4216 0.3630 .01254 .26070 (2.9) 2.9259 0.0136 8.1716 0.3000 7.3279 0.3552 .01185 .29188 85071 X X 2.6872 0.0260 7.5320 0.1610 7.5320 0.1610 .10697 .63053 49 sample misalignment. Also, results obtained using (2.7) would be adversely affected by any uncertainty in the phase planes? positions. 2.7 Summary In this chapter, the performance of complex permittivity extraction methods based on SQP and the GA was contrasted using S-parameter measurements for 1-, 2-, and 3-layer samples. Three different error function definitions were also used to quantify the performance of each algorithm in terms of the amount of S-parameter information (magnitude only or magnitude and phase) available for the inversion process. Computer generated S-parameter data was initially used to determine the attainable accuracy of each algorithm/error function combination. The results of this portion of the study clearly indicated that the extracted permittivity from the single and multilayer cases was nearly identical to that used to generate the data when the GA was used. This was also found to be true for SQP in all but two cases possibly due to local minima trapping. This demonstrates that the GA is extremely accurate and would therefore be limited only by the precision of the S-parameter measurements, whereas the performance of SQP would also likely be limited by the accuracy of the initial guess. The algorithm/error functions were used to extract the complex permittivity from single layer S-parameter measurements, and it was evident that all of the optimization techniques were highly effective at minimizing the RMS error(s) of the S-parameters (both phase and magnitude). However, the GA using either (2.8) or (2.9) was found to consistently produce the lowest magnitude error of the S-parameters and hence the best estimate of the complex permittivity. 50 When the same techniques were used for complex permittivity extraction from mul- tilayer composite structures, the GA, for all cases considered, provided approximately the same level of accuracy as that observed for the single layer cases. SQP, however, failed to obtain accurate results for several of the cases considered. In summary, the GA appeared to be the more robust algorithm in terms of its ability to always achieve a low S-parameter RMS error and accurately obtain each layer?s complex permittivity. 51 Chapter 3 Complex Constitutive Parameter Extraction 3.1 Overview This chapter extends the procedure presented in Chapter 2 to allow extraction of both complex permittivity and permeability for each layer in a multilayer sample. In order to determine the accuracy of the method, the extracted constitutive parameters (CP) are used to generate S-parameters for comparison with measured S-parameters. Also, the single layer, complex CP extraction technique developed by Wolfson and Went- worth [96,97] is employed to provide further verification of the correct operation of the code. X band waveguide S-parameter measurements of three materials (Teflon, F40, and F125) were obtained using an HP-8510C Vector Network Analyzer for the frequency range of 8.2 - 10 GHz. Teflon was used to ensure the method obtained accurate results for nonmagnetic and low loss materials. Further, two radar-absorbing materials (RAM), F40 and F125 [98], were used to demonstrate extraction of complex permittivity and permeability. Multilayer samples were constructed by placing the single layer samples adjacent to one another in different combinations. 52 3.2 Extraction Techniques 3.2.1 Wolfson-Wentworth Method The Wolfson-Wentworth method consists of placing a sample in a section of rectan- gular waveguide and measuring the S11 parameter with the waveguide terminated by two offset shorts of differing length [96,97]. The S-parameter measurements are required to calculate the input impedance of the sample for each case. The input impedances along with transmission line equations are then used to extract the complex permittivity and permeability of the sample over the frequency range of interest. This approach requires the material sample thickness to be less than one-half wavelength for the TE10 mode inside the sample in order to avoid exciting higher order modes. However, a low loss sample must be thick enough to provide significant reflections. If both of these conditions cannot be met, spurious data for the extracted values of permittivity and permeability can result at high frequencies. A detailed description of this method is given in [97]. 3.2.2 CP Extraction Method Modifications As discussed previously, the method of Chapter 2 is a 2-port technique requiring a full set of S-parameters to extract complex permittivity for each layer of an n-layer sample (see Fig. 2.1). Here, however, the method is modified to also account for mag- netic materials. Appropriate modification of the forward solution only requires that the propagation constant and wave impedance of (2.1) be expressed as ?i = radicalbigg (?/a)2 ??2?0?0 parenleftBig ??r,i ?j???r,i parenrightBigparenleftBig ??r,i ?j???r,i parenrightBig (3.1) 53 Zi = j??0 parenleftBig ??r,i ?j???r,i parenrightBig ?i (3.2) where a is the maximum cross-sectional dimension of the waveguide. Calculation of the S-parameters then follows the same procedure outlined in Section 2.2. Once again, an error function analysis was carried out due to the fact that the extraction both both ?? and ?? results in a higher dimensional search space. Therefore, the conclusions previously drawn may be invalid for this situation. Initially, Err2 (2.8) was considered and determined to be ill-suited for accurately estimating the CP since S-parameter phase information becomes critical to the extraction process?s ability to obtain low RMS errors. Therefore, two error functions involving both magnitude and phase information were tested and are given by Err = 1N Nsummationdisplay i=1 parenleftBigvextendsinglevextendsingle vextendsingle[S]f,i vextendsinglevextendsingle vextendsingle? vextendsinglevextendsingle vextendsingle[S]m,i vextendsinglevextendsingle vextendsingle parenrightBig2 + 1N Nsummationdisplay i=1 parenleftBig negationslash [S]f,i ?negationslash [S]m,i parenrightBig2 (3.3) Err = 1N Nsummationdisplay i=1 parenleftBig Re parenleftBig [S]f,i parenrightBig ?Re parenleftBig [S]m,i parenrightBigparenrightBig2 + 1N Nsummationdisplay i=1 parenleftBig Im parenleftBig [S]f,i parenrightBig ?Im parenleftBig [S]m,i parenrightBigparenrightBig2 (3.4) where [S]f,i and [S]m,i are the formulated and measured S-parameters at frequency point i, respectively, and N is the number of frequencies. After a number of studies, (3.4) was found to give a lower overall RMS error for both magnitude and phase of the formulated 54 and measured S-parameters. This may be due to a number of factors such as numer- ical precision of the calculations and unequal weighting of the magnitude and phase information in (3.3). The Genetic Algorithm (GA) of Chapter 2 was determined to be an extremely robust method for determining accurate values for complex permittivities from S-parameter measurements. Sequential Quadratic Programming (SQP), a local optimization tech- nique, was only able to accurately determine permittivity values in some cases due to its severe dependence on the algorithm?s initial starting point (??initial). However, SQP was found to be 50 to 1700 times faster than the GA depending on the number of layers (higher speed-up for larger number of layers). Also, SQP was shown to more accurately obtain the value of the global minimum than the GA when the initial starting point was in the vicinity of the minimum. In this section, a modified SQP algorithm is developed to exploit SQP?s speed and accuracy while eliminating the issue of local minima trapping. The novel multi-point SQP (MPSQP) presented here relies on the generation of P randomly1 distributed initial guesses (??initial, ??initial) to reach the global minimum in the bounded solution space. The SQP algorithm is performed on each initial guess resulting in P solutions. If a sufficient number of points are taken, the MPSQP will accurately and quickly determine the global minimum by taking the solution with minimum error of all the returned solutions. To ensure that the minimum has been accurately determined, the MPSQP algorithm is repeated but bounded by a reduced search space centered about the previously determined solution. In all cases considered, this added step returned results identical to the results of the initial set of iterations. A study comparing the performance 1Uniform probability density function. 55 of the GA and MPSQP for the extraction of single layer material parameters is presented in Section 3.4.1. 3.3 Measurements The test ports of an HP-8510C Vector Network Analyzer were connected to WR90 (X band) waveguide via coax-to-waveguide adapters. The calibration procedure follows Wolfson?s outline of recommended instrument settings for improved measurement accu- racy [99]. Higher accuracy was obtained when operating the analyzer in step sweep mode with 128X averaging. A 2 ms dwell time was used to account for the analyzer settling time and propagation delay due to coaxial cable length. Finally, the TE10 mode cutoff frequency is entered as the waveguide delay (6.557 GHz for the WR90 waveguide used in this work). For the Wolfson-Wentworth method, the measurement reference plane is established at the open end of the waveguide of sufficient length to allow unwanted modes to attenu- ate before reaching the measurement reference plane. The reference plane was defined at the end of this waveguide using an offset short calibration procedure [100]. The procedure uses a short at the reference plane and offset shorts of lengths approximately ?/8 and 3?/8, where ? is chosen to give maximum phase separation for the offset shorts across the band [100]. For WR90 rectangular waveguide the optimum offset short lengths become .483 cm and 1.455 cm. Following calibration, S11 is measured for the reference plane terminated by both offset short loads. A short section of waveguide is then attached to the reference plane, with the material sample placed in the far end of the guide. S11 is then measured for the sample terminated by both offset short loads. These values 56 along with S11 measured for the offset shorts are inserted into the routine described by Wolfson [97] to extract ?? and ??. A through-response measurement for the method presented here requires two sec- tions of waveguide for attenuation of unwanted modes. The measurement reference planes are established at the ends of these two waveguides using a full 2-port calibra- tion that employs a short, the pair of offset shorts described in the Wolfson-Wentworth approach, and a direct connection. Following calibration, the material sample is placed into a short section of waveguide which is then inserted between the two reference planes. A complete set of S-parameters are then measured and the results used in the method presented here. 3.4 Results 3.4.1 Single Layer Measurements For the single layer materials, samples of Teflon (d = 4.8 mm), F40 (d = 1.524 mm), and F125 (d = 3.3 mm) were used. S-parameter data sets were calculated from the values extracted using the GA with the appropriately modified forward solver, MPSQP, and the Wolfson-Wentworth method. These data sets were in turn compared to the measured S-parameters. In addition to S-parameter comparisons, the extracted CP from each algorithm were directly compared. It was initially assumed for the GA and MPSQP that the CP were constant over the frequency band of interest, whereas the Wolfson- Wentworth method makes no such assumptions. The results presented show that this assumption was valid, and, therefore, it is used throughout this chapter. For materials having frequency dependent CP, the algorithm can be easily modified. 57 For all single layer optimizations, the GA utilized a population size of 100, crossover rate of 80%, and mutation rate of 10%. The redundancy removal scheme of the TRRGA was also employed to accelerate convergence. Since a single layer problem has a four dimensional search space, the GA was allowed to iterate for 1000 generations which corresponds to the parameter utilized in Section 2.6.2 for 4-D problems. For MPSQP, the number of initial guesses (P) was set to 100. The real parts of both the relative permittivity and permeability were restricted to the range of .1 to 25, while the imaginary parts were restricted from 0 to 8. For the Teflon sample, the extracted constitutive parameters are nearly identical for the GA and MPSQP as shown in Table 3.1. Good agreement is also shown between these two methods and the results of the Wolfson-Wentworth method. The constitutive parameter values listed for the Wolfson-Wentworth method are the average values over the frequency band. The RMS errors between the MPSQP-generated and measured S-parameter magnitudes and phases are listed in Table 3.2 and show excellent agree- ment. The GA results were nearly identical to the MPSQP results and are not shown. Analogous results were also obtained for the two RAM samples as shown in the tables. In all cases, MPSQP returned slightly lower RMS S-parameter errors (generated versus measured) than the Wolfson-Wentworth method and the GA. Figs. 3.1 and 3.2 show comparisons of the extracted CP from MPSQP and the Wolfson-Wentworth method over the entire frequency range. Figs. 3.3 and 3.4 show excellent agreement between the MPSQP-generated S-parameters and the measured data. Table 3.3 illustrates the improvement in performance when employing MPSQP in- stead of the GA. Overall, the MPSQP showed an average speed-up of 21.1 (speed-up 58 taken as the ratio of GA runtime to MPSQP runtime). Additionally, the MPSQP was able, in all cases, to achieve slightly lower error function values which, in turn, directly corresponded to lower S-parameter RMS errors. It is also worth noting that in Chapter 2 the speed-up of SQP over the GA scaled superlinearly with increasing number of layers. Although not directly presented here, this was also observed with MPSQP and the GA for the multilayer cases presented in the following sections. Therefore, it appears that the MPSQP is an improved method for obtaining the material parameter values not only in terms of its speed but also its accuracy in determining the global minimum. Table 3.1: Single Layer Results Method Teflon F40 F125 GA ?? 2.05 12.7 6.97 MPSQP ?? 2.03 12.7 6.98 Wolfson ?? 2.07 12.9 7.95 GA ??? 0 0.187 0 MPSQP ??? 0 0.203 0 Wolfson ??? 0.027 0.673 0.693 GA ?? 0.958 1.55 0.440 MPSQP ?? 0.965 1.55 0.434 Wolfson ?? 0.952 1.63 0.490 GA ??? 4.52e-3 1.24 0.516 MPSQP ??? 4.63e-3 1.24 0.515 Wolfson ??? 0.017 1.06 0.474 59 Table 3.2: MPSQP Error for Single Layer Sample Teflon F40 F125 RMS Error |S11| 1.47e-2 6.60e-3 1.25e-2 RMS Error |S12| 8.17e-3 5.14e-3 7.86e-3 RMS Error |S21| 8.20e-3 3.72e-3 1.43e-2 RMS Error |S22| 1.29e-2 4.35e-3 7.04e-3 RMS Error negationslash S11 2.85e-2 2.60e-2 2.08e-2 RMS Error negationslash S12 1.01e-2 2.05e-2 4.71e-2 RMS Error negationslash S21 1.50e-2 1.41e-2 3.97e-2 RMS Error negationslash S22 4.82e-2 2.77e-2 2.15e-2 Table 3.3: Performance Comparison of GA and MPSQP Error Time (sec)Material GA MPSQP GA MPSQP Teflon 1.53e-3 1.47e-3 810.2 48.6 F40 8.93e-4 8.93e-4 782.6 32.1 F125 1.35e-3 1.34e-3 790.9 35.6 60 8.5 9 9.5 10?2 0 2 4 6 8 10 12 Frequency (GHz) ??, ??? MPSQP Wolfson Method MPSQP ?? Wolfson ?? MPSQP ???Wolfson ??? Figure 3.1: Comparison of MPSQP and Wolfson-Wentworth method extracted complex permittivity values for F125 sample. 61 8.5 9 9.5 100 0.2 0.4 0.6 0.8 1 Frequency (GHz) ??, ??? MPSQP Wolfson Method MPSQP ?? Wolfson ??MPSQP ??? Wolfson ??? Figure 3.2: Comparison of MPSQP and Wolfson-Wentworth method extracted complex permeability values for F125 sample. 62 8.5 9 9.5 100 0.2 0.4 0.6 0.8 1 Frequency (GHz) |S 11 |, |S 21 | MPSQP Wolfson Method Measurements |S11| |S21| Figure 3.3: Measured and generated S-parameter magnitudes for the F125 sample. 63 8.5 9 9.5 10?200 ?150 ?100 ?50 0 50 100 150 Frequency (GHz) ?S 11 , ? S 21 MPSQP Wolfson Method Measurements ?S11 ?S21 Figure 3.4: Measured and generated S-parameter phases for the F125 sample. 64 3.4.2 Two Layer Measurements ThesamplesusedfortwolayerS-parametermeasurementswereF125/Teflon, Teflon/F40, and F125/F40. To ensure an accurate solution, P was set to 500 initial points. The in- crease in P was necessary since the search space for a two layer problem has 8 dimensions. The upper and lower bounds remained the same as those of the single layer case. Table 3.4 shows the extracted parameter values as well as the error function value and RMS errors between the generated and measured S-parameters. Overall, excellent agreement is shown by the low error values for all material combinations. Additionally, the extracted parameter values agree well with the values obtained from the single layer optimization. Figs. 3.5 and 3.6 show the magnitude and phase of the generated and measured S11 and S21 for the F125/Teflon sample. As with the data from Table 3.4, the agreement between the extracted and measured data is well within the tolerance of instrumentation error. 65 Table 3.4: Results for Two Layer Samples Sample F125/Teflon Teflon/F40 F125/F40 Layer F125 Teflon Teflon F40 F125 F40 ?? 6.82 2.05 2.36 11.8 6.2 14.0 ??? 0 0 0 0.091 0 0 ?? 0.444 0.997 1.08 1.06 0.402 0.492 ??? 0.546 0 9.72e-3 1.23 1.78 1.46 Error Function Value 2.29e-3 3.02e-3 2.73e-3 RMS Error |S11| 1.85e-2 7.30e-3 7.99e-3 RMS Error |S12| 1.39e-2 8.21e-3 1.04e-2 RMS Error |S21| 1.17e-2 1.03e-2 1.04e-2 RMS Error |S22| 2.02e-2 2.57e-2 1.85e-2 RMS Error negationslash S11 2.88e-2 2.90e-2 4.96e-2 RMS Error negationslash S12 4.91e-2 2.82e-2 6.12e-2 RMS Error negationslash S21 4.40e-2 4.25e-2 6.83e-2 RMS Error negationslash S22 2.64e-2 5.87e-2 5.99e-2 66 8.5 9 9.5 100 0.2 0.4 0.6 0.8 1 Frequency (GHz) |S 11 |, |S 21 | MPSQP Measurements |S11| |S21| Figure 3.5: Measured and generated S-parameter magnitudes for the F125/Teflon sam- ple. 67 8.5 9 9.5 10?250 ?200 ?150 ?100 ?50 0 50 100 150 Frequency (GHz) ?S 11 , ? S 21 MPSQP Measurements ?S11 ?S21 Figure 3.6: Measured and generated S-parameter phases for the F125/Teflon sample. 68 3.4.3 Three Layer Measurements A three layer F125/F40/Teflon sample?s S-parameters were measured and the com- plex constitutive parameters extracted using MPSQP. The upper and lower bounds re- mained unchanged while P was increased to 1000 due to the 12 dimensional search space. Table 3.5 shows the extracted constitutive parameters for each layer as well as the error function value and RMS errors between the measured and generated S-parameters. As in the prior cases, the algorithm was able to successfully match the generated and mea- sured S-parameters as evidenced by the very low RMS error values. Inspection of the extracted parameter values shows that each of the materials was correctly classified and agreement with the single and two layer parameter values was very good. Table 3.5: Results for F125/F40/Teflon Sample F125 F40 Teflon ?? 5.98 12.1 2.04 ??? 0 0 0.107 ?? 0.411 1.90 1.08 ??? 0.495 1.49 0 Error Function Value 7.61e-4 RMS Error |S11| 7.03e-3 RMS Error |S12| 8.77e-3 RMS Error |S21| 8.48e-3 RMS Error |S22| 1.10e-2 RMS Error negationslash S11 1.22e-2 RMS Error negationslash S12 5.78e-2 RMS Error negationslash S21 4.79e-2 RMS Error negationslash S22 2.59e-2 69 3.5 Summary In this chapter, a method was presented for accurately extracting the complex per- mittivity and permeability from each individual layer in a multilayer sample using S- parameter waveguide measurements. The forward solution formulation presented in Chapter 2 was modified to account for magnetic materials. Also, a modified SQP al- gorithm (MPSQP) was presented which utilized a large number initial guess points to alleviate the possibility of local minima trapping, a problem with gradient-based opti- mization methods. Use of such an algorithm was beneficial since it provides significant computational gains over traditional global optimization methods such as the GA. Specif- ically, for the single layer cases presented here, the MPSQP showed an average speed-up of 21.1 over the GA as well as improved accuracy. S-parameter measurements were conducted on three material samples used to con- struct multilayer samples. The MPSQP was used successfully to extract the complex CP for each layer. These values were then compared with values extracted using the Wolfson- Wentworth method (single layer cases only). Also, S-parameters were generated using the extracted values and compared with the measured data. In all cases, results were found to be in excellent overall agreement with both the Wolfson-Wentworth method values and measured data. In summary, the MPSQP was found to be a computationally efficient and robust algorithm for extracting complex CP from multilayer materials. 70 Chapter 4 Waveguide Filters Using Conductor Doping of Dielectrics 4.1 Overview In the previous chapters, novel optimization algorithms were presented for the so- lution of inverse problems, namely, complex constitutive parameter extraction of mul- tilayered materials. In this and the remaining chapters, the topic of optimal design is addressed. The specific problem considered is that of optimizing the geometry of waveguide filters to realize a specified transmission response. Novel frequency-dependent reflectivity and transmittivity characteristics of dielec- tric slabs have recently been shown to be achievable by embedding randomly oriented thin conducting wires (inclusions) within a host slab [101]. A wide variety of electro- magnetic behaviors have been observed by adjusting the wire length and doping level (wire density). In addition to wires, thin conducting patches or small particles of var- ious geometries may also be used as the dopant and, in the case of wires and patches, either embedded or limited to the surface of the dielectric. The focus of this study is to investigate the utility of random doping as a method for obtaining consistent reflection and transmission frequency-dependent profiles. Additional consideration is given to the feasibility of fabricating such materials. 4.2 Numerical Modeling To determine the reflection and transmission characteristics of a doped dielectric slab, a numerical procedure is necessary. Many numerical methods including Finite 71 Element Method (FEM) (and hybrids) [101?109], Method of Moments (MoM) [110?113], and Finite Difference Time Domain (FDTD) [114?122] have been applied to similar problems with great success. However, edge-based FEM is most efficient for the problem at hand because of its ability to easily account for inhomogeneous materials as well as its suitability for closed region problems (i.e., waveguides, cavities, etc.) [123]. The computational domain of the FEM is created by placing the dielectric slab in- side a waveguide as shown in Fig. 4.1 and discretizing the geometry using a tetrahedral elements. In this manner, the reflection and transmission properties of the material can be obtained in a computationally efficient manner. A mode-matching (MM) scheme, which allows for the development of TE and TM modes, is used to truncate the com- putational domain and permits the faces of the waveguide to be placed very close to the slab [123,124]. Further reduction in computational intensity is achieved by utilizing a multi-point Asymptotic Waveform Evaluation (AWE) [123] to quickly determine the frequency response of the material over a wide band. The matrix solution of the FEM problem is used to calculate the reflection (S11) and transmission (S21) coefficients which define the filter?s response. To further simplify the geometry and mesh generation as well as reduce solution time, the following methodology is employed to generate the doped medium. In order to model thin wires or conducting triangular patches, assume that the locations of the wires or patch edges occupy the same locations as tetrahedral edges within the discretized slab (see Fig. 4.2). Those edges where wires or patches exist can simply be assigned tangential electric field values of zero. Since the unknowns in an edge-based FEM are the tangential electric fields along the tetrahedral edges, this is equivalent to removing the appropriate 72 x y z a b Doped Slab 2 1 3 ? ? ? Figure 4.1: Loaded waveguide geometry with doped slab. Boundary notations for FEM formulation are indicated. 73 rows and columns of the FEM matrix equation which results in a reduced set of linear equations. Additionally, by this method, several different cases may be considered us- ing the same matrix equation. All that is required is the appropriate removal of rows and columns of the matrix, and, hence, no matrix reassembly is required. This results in a computationally efficient procedure, especially beneficial to optimization problems requiring many simulations. Although the edge locations within a given mesh are fixed, modern unstructured mesh generators provide enough variety in edge orientation that the wires/patches may still be assumed to be randomly placed [101]. Also notice that, the wires are assumed to be infinitesimally thin elements and there is no provision to include wire radius in the current formulation. To better illustrate this process, consider the following example. To generate a doped slab with 100 wires, a random number generator is used to select 100 numbers from a list numbering 1 to Nedges, the number of edges within the discretized slab. This list maps to the corresponding edge numbers which are then removed from the matrix equation as described. Generating a list of edges corresponding to conducting patches is performed in the same manner with the additional requirement that the three connecting edges composing a tetrahedral face form one patch. Therefore, a doping level of 100 triangular patches would correspond to the removal of 300 rows and columns from the FEM matrix equation. 4.3 Finite Element Formulation The FEM presented here is formulated to solve for the S-parameters from a loaded waveguide (see Fig. 4.1). Inside the waveguide, electromagnetic wave propagation is governed by Maxwell?s equations. From these equations, the vector Helmholtz equation 74 X58 X59 X5A Figure 4.2: Discretized slab illustrating wires lying along tetrahedral edges. Notice some edges join to form long, erratically shaped wires. Wires may also be located on interior edges but are not shown here. 75 for the electric field can be derived and is given by ????1?? vectorE ??2?vectorE = 0 (4.1) where vectorE is the electric field inside the waveguide, ? is the radian frequency of the wave, and ? and ? are the permeability and permittivity, respectively, at any given point in the domain. Within the waveguide, the electric field is subject to the boundary conditions ?n? vectorE = 0 on ?1 (4.2) ?n??? vectorE +??n? parenleftBig ?n? vectorE parenrightBig = vectorU on ?2,3 (4.3) Eqn. (4.2) forces the tangential electric field to be zero along the waveguide walls. Eqn. (4.3) is referred to as a boundary condition of the third kind. Boundary con- ditions of this type are necessary to accurately truncate the computational domain of the problem, and a discussion of their parameters will be given in a later section. Phys- ically, (4.3) represents a relationship between the tangential electric and magnetic fields along the boundary. As with most numerical methods, the FEM procedure begins with a discretization of the geometry into elements. Although there are many elements available, tetrahedral elements, as shown in Fig. 4.3, are usually chosen due to their ability to conform to arbitrarily shaped structures. Inside each element, the electric field can be approximated by an expansion such as vectorEe = 6summationdisplay j=1 vectorNj (x,y,z)Ej (4.4) 76 where vectorEe is the electric field inside the eth element, vectorNj (x,y,z) is the vector basis function associated with the jth tetrahedral edge, and Ej is defined as the electric field value at the center of edge j projected along that edge?s direction. The vector basis functions used in this work are termed the zeroth-order vector basis functions and given by vectorNj = ?j (Lj1?Lj2 ?Lj2?Lj1) (4.5) where ?j is the length of edge j and Lj? is the linear interpolation function (volume coordinate) for the ?th node bounding edge j (see Fig. 4.3) [123]. The necessity of edge- based (vector) elements as opposed to more traditional node-based elements has been discussed extensively in the literature [6,103,125?155]. Following Galerkin?s procedure [123], the weighted residual of (4.1) is minimized according to Rei = integraldisplay ?e vectorNi ?parenleftBig????1?? vectorE ??2?vectorEparenrightBig d?e (4.6) where Rei is the weighted residual corresponding to edge i of element e; vectorNi, defined by (4.5), is the weighting function for edge i; and ?e represents the domain of element e. Using the first vector Green?s theorem, (4.6) can be cast in the weak form which is given by Rei = integraldisplay ?e ?? vectorNi ??? vectorE ??2??vectorNi ? vectorE d?e ? contintegraldisplay ?e parenleftBigvector Ni ??? vectorE parenrightBig ? ?n d?e (4.7) 77 1 2 3 4 I II III IV V VI Figure 4.3: Tetrahedral element showing edge and node numbering. 78 From a vector identity, the integrand of the surface integral in (4.7) can be rewritten in two different ways as shown by parenleftBigvector Ni ??? vectorE parenrightBig ? ?n = parenleftBig ?n? vectorNi parenrightBig ? parenleftBig ?? vectorE parenrightBig = ?vectorNi ? parenleftBig ?n??? vectorE parenrightBig (4.8) Since Galerkin?s method calls for the weighting and basis functions to be identical, vectorNi must satisfy the same boundary conditions as the electric field. Therefore, the second term in (4.8) shows that the integrand is zero for the portion of the surface integral over ?1. Using the third form in (4.8), the weight residual can be rewritten as Rei = integraldisplay ?e ?? vectorNi ??? vectorE ??2??vectorNi ? vectorE d?e + integraldisplay ?e2,3 vectorNi ?parenleftBig?n??? vectorEparenrightBig d?e2,3 (4.9) 4.3.1 Absorbing Boundary Condition Many domain truncation schemes exist as outlined in [123]. A simple first-order absorbing boundary condition (ABC) method is presented here followed by a more accu- rate Mode-Matching (MM) technique. Substituting (4.3) for the integrand of the surface integral in (4.9), the weighted residual is given by Rei = integraldisplay ?e ?? vectorNi ??? vectorE ??2??vectorNi ? vectorE d?e + integraldisplay ?e2,3 vectorNi ?parenleftBigvectorU ???n? ?n? vectorEparenrightBig d?e2,3 (4.10) The remaining task is to determine appropriate values for ? and vectorU. Assuming the surfaces ?2 and ?3 are sufficiently far from the discontinuity in the waveguide, it can be assumed that only the TE10 mode is propagating (no higher order or evanescent modes). 79 It is also assumed that a unit amplitude TE10 mode wave is incident on ?2. The total field on ?2 may then be written as vectorE = vectorEinc + vectorEref = sinparenleftBig?x a parenrightBig e?jkzz1 +Rsin parenleftBig?x a parenrightBig ejkzz1 (4.11) whereRisthereflectioncoefficient, kz isthepropagationconstantoftheemptywaveguide, and z1 is the position of ?2. Applying (4.11) to the first term of the left-hand side of (4.3), the following is determined: ?n??? vectorE = ??z ??? vectorE = ?jkz vectorEinc +jkz vectorEref = jkz vectorE ?2jkz vectorEinc (4.12) Equating this with (4.3), it is apparent that ? = jkz and vectorU = ?2jkz vectorEinc since the total field is assumed to be in the form of the TE10 mode. Along ?3, a similar analysis can be performed. The field along ?3 can be expressed as vectorE = vectorEtr = T sinparenleftBig?x a parenrightBig e?jkzz2 (4.13) where T is the transmission coefficient and z2 is the location of ?3. Applying the same procedure, it is easily determined that, along ?3, ? = jkz and vectorU = 0. Substituting these results into (4.10), the weighted residual becomes Rei = integraldisplay ?e ?? vectorNi ??? vectorE ??2??vectorNi ? vectorE d?e ? integraldisplay ?e2,3 vectorNi ???n? ?n? vectorE d?e2,3 + integraldisplay ?e2 vectorNi ? vectorU d?e2 (4.14) 80 which can be rewritten as Rei = integraldisplay ?e ?? vectorNi ??? vectorE ??2??vectorNi ? vectorE d?e + integraldisplay ?e2,3 ??n? vectorNi ? ?n? vectorE d?e2,3 ? integraldisplay ?e2 ?n? vectorNi ? vectorU ? ?n d?e2 (4.15) Substituting (4.4) into (4.15), the final form of the residual is obtained and given by Rei = 6summationdisplay j=1 Ej integraldisplay ?e ?? vectorNi ??? vectorNj ??2??vectorNi ? vectorNj d?e + integraldisplay ?e2,3 ??n? vectorNi ? ?n? vectorNj d?e2,3 ? integraldisplay ?e2 ?n? vectorNi ? vectorU ? ?n d?e2 (4.16) Following the matrix assembly procedure outlined in [123], the final matrix equation can be obtained and is written as ([S]?[T]+[Q]){E} = {b} (4.17) The entries in the local element matrices are, therefore, given by Seij = integraldisplay ?e ?? vectorNi ??? vectorNj d?e Teij = integraldisplay ?e ?2??vectorNi ? vectorNj d?e Qeij = integraldisplay ?e2,3 ??n? vectorNi ? ?n? vectorNj d?e2,3 bei = integraldisplay ?e2 ?n? vectorNi ? vectorU ? ?n d?e2 (4.18) 81 4.3.2 Mode-Matching Domain Truncation The above domain truncation method requires at least one wavelength separation of ?2 and ?3 from the discontinuity in the waveguide [123]. This method is sufficient when the frequency is low or when the waveguide load does not require a significant number of unknowns. For complex structures or high frequencies, the ABC method is somewhat inefficient. Therefore, an expansion of the unknown field at the boundary in terms of the orthogonal waveguide modes (eigenfunctions) presents a means by which the computational domain may be significantly reduced. For the MM technique [6], the electric field along ?2 can be expressed by vectorE = vectorEinc + ?summationdisplay m=0 ?summationdisplay n=0 amnvectorETEmne?mnz1 + ?summationdisplay m=1 ?summationdisplay n=1 bmn bracketleftBigvector ETMmn + ?zETMzmn bracketrightBig e?mnz1 (4.19) where amn and bmn are expansion coefficients and ?mn is the propagation constant of the specified mode in the empty waveguide. The terms vectorETEmn, vectorETMmn , and ETMzmn are the normalized waveguide modes and are given by vectorETEmn = ?? m?n nmn parenleftBign b cos parenleftBigm?x a parenrightBig sin parenleftBign?y b parenrightBig ?x? ma sin parenleftBigm?x a parenrightBig cos parenleftBign?y b parenrightBig ?y parenrightBig (4.20) vectorETMmn = 2 nmn parenleftBigm a cos parenleftBigm?x a parenrightBig sin parenleftBign?y b parenrightBig ?x+ nb sin parenleftBigm?x a parenrightBig cos parenleftBign?y b parenrightBig ?y parenrightBig (4.21) ETMzmn = 2k 2c ??mnnmn sin parenleftBigm?x a parenrightBig sin parenleftBign?y b parenrightBig (4.22) where ?? is Neumann?s number, nmn = radicalBig n2ab +m2 ba, and kc is the cutoff wavenumber of the specified mode. Note that (4.20-4.22) are orthogonal according to the following 82 relationships: bintegraldisplay 0 aintegraldisplay 0 vectorETEmn ? vectorETEm?n?dxdy = ?mm??nn? (4.23) bintegraldisplay 0 aintegraldisplay 0 vectorETMmn ? vectorETMm?n?dxdy = ?mm??nn? (4.24) bintegraldisplay 0 aintegraldisplay 0 vectorETEmn ? vectorETMm?n?dxdy = 0 (4.25) bintegraldisplay 0 aintegraldisplay 0 ETMzmnETMzm?n?dxdy = ?mm??nn? k 2c ?2mn (4.26) Therefore, the coefficients of (4.19) may be expressed according to the following deriva- tion. Firstly, the reflected electric field is written as vectorEref = ?summationdisplay m=0 ?summationdisplay n=0 amnvectorETEmne?mnz1 + ?summationdisplay m=1 ?summationdisplay n=1 bmn bracketleftBigvector ETMmn + ?zETMzmn bracketrightBig e?mnz1 (4.27) from which amn may be obtained by taking a vector dot product on both sides by vectorETEm?n?e??m?n?z1 and integrating over ?2 as is shown by integraldisplay ?2 vectorETEm?n?e??m?n?z1 ? vectorErefd?2 = bintegraldisplay 0 aintegraldisplay 0 vectorETEm?n?e??m?n?z1 ? vectorErefdxdy = ?summationdisplay m=0 ?summationdisplay n=0 bintegraldisplay 0 aintegraldisplay 0 vectorETEm?n?e??m?n?z1 ?amnvectorETEmne?mnz1dxdy + ?summationdisplay m=1 ?summationdisplay n=1 bintegraldisplay 0 aintegraldisplay 0 vectorETEm?n?e??m?n?z1 ?bmnbracketleftBigvectorETMmn + ?zETMzmnbracketrightBige?mnz1dxdy (4.28) 83 which results in amn = integraldisplay ?2 vectorETEmne??mnz1 ? vectorErefd?2 (4.29) Similarly, bmn is found to be bmn = integraldisplay ?2 vectorETMmn e??mnz1 ? vectorErefd?2 (4.30) Following a procedure similar to that of the previous section, (4.12) becomes ?n??? vectorE = ?n??? vectorEinc + ?summationdisplay m=0 ?summationdisplay n=0 amn?mnvectorETEmne?mnz1 ? ?summationdisplay m=1 ?summationdisplay n=1 bmn k 20 ?mn vectorETMmn e?mnz1 (4.31) which by substitution of (4.29) and (4.30) into (4.31) can be written in terms of the total field and incident field as ?n??? vectorE = vectorP + vectorU (4.32) where vectorP = ?summationdisplay m=0 ?summationdisplay n=0 ?mnvectorETEmn integraldisplay ?2 vectorETEmn ? vectorEd?2 ? ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn vectorETMmn integraldisplay ?2 vectorETMmn ? vectorEd?2 (4.33) vectorU = ?n??? vectorEinc ? ?summationdisplay m=0 ?summationdisplay n=0 ?mnvectorETEmn integraldisplay ?2 vectorETEmn ? vectorEincd?2 + ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn vectorETMmn integraldisplay ?2 vectorETMmn ? vectorEincd?2 (4.34) 84 Similarly, on ?3, (4.32) is used with vectorU = 0. Therefore, (4.9) can now be expressed by Rei = integraldisplay ? ?? vectorNi ??? vectorE ??2??vectorNi ? vectorE d? + integraldisplay ?2 vectorNi ?parenleftBigvectorP + vectorUparenrightBig d?2 + integraldisplay ?3 vectorNi ?parenleftBigvectorPparenrightBig d?3 (4.35) where vectorP and vectorU are evaluated on the appropriate surfaces. When expanded, (4.35) is given by Rei = integraldisplay ? ?? vectorNi ??? vectorE ??2??vectorNi ? vectorEd? + ?summationdisplay m=0 ?summationdisplay n=0 ?mn integraldisplay ?2,3 vectorNi ? vectorETEmnd?2,3 integraldisplay ?2,3 vectorETEmn ? vectorEd?2,3 ? ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn integraldisplay ?2,3 vectorNi ? vectorETMmn d?2,3 integraldisplay ?2,3 vectorETMmn ? vectorEd?2,3 + integraldisplay ?2 vectorNi ?parenleftBig?n??? vectorEincparenrightBigd?2 ? ?summationdisplay m=0 ?summationdisplay n=0 ?mn integraldisplay ?2 vectorNi ? vectorETEmnd?2 integraldisplay ?2 vectorETEmn ? vectorEincd?2 + ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn integraldisplay ?2,3 vectorNi ? vectorETMmn d?2,3 integraldisplay ?2,3 vectorETMmn ? vectorEincd?2,3 (4.36) 85 Substituting (4.4) yields Rei = 6summationdisplay j=1 Ej integraldisplay ? ?? vectorNi ??? vectorNj ??2??vectorNi ? vectorNjd? + ?summationdisplay m=0 ?summationdisplay n=0 ?mn integraldisplay ?2,3 vectorNi ? vectorETEmnd?2,3 integraldisplay ?2,3 vectorETEmn ? vectorNjd?2,3 ? ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn integraldisplay ?2,3 vectorNi ? vectorETMmn d?2,3 integraldisplay ?2,3 vectorETMmn ? vectorNjd?2,3 + integraldisplay ?2 vectorNi ?parenleftBig?n??? vectorEincparenrightBigd?2 ? ?summationdisplay m=0 ?summationdisplay n=0 ?mn integraldisplay ?2 vectorNi ? vectorETEmnd?2 integraldisplay ?2 vectorETEmn ? vectorEincd?2 + ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn integraldisplay ?2,3 vectorNi ? vectorETMmn d?2,3 integraldisplay ?2,3 vectorETMmn ? vectorEincd?2,3 (4.37) which can be written in the matrix form of (4.17). Now, however, the elements Qij and bi are expressed by Qij = ?summationdisplay m=0 ?summationdisplay n=0 ?mn integraldisplay ?2,3 vectorNi ? vectorETEmnd?2,3 integraldisplay ?2,3 vectorETEmn ? vectorNjd?2,3 ? ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn integraldisplay ?2,3 vectorNi ? vectorETMmn d?2,3 integraldisplay ?2,3 vectorETMmn ? vectorNjd?2,3 (4.38) 86 bi = ? integraldisplay ?2 vectorNi ?parenleftBig?n??? vectorEincparenrightBigd?2 + ?summationdisplay m=0 ?summationdisplay n=0 ?mn integraldisplay ?2 vectorNi ? vectorETEmnd?2 integraldisplay ?2 vectorETEmn ? vectorEincd?2 ? ?summationdisplay m=1 ?summationdisplay n=1 k20 ?mn integraldisplay ?2,3 vectorNi ? vectorETMmn d?2,3 integraldisplay ?2,3 vectorETMmn ? vectorEincd?2,3 (4.39) The resulting matrix equation can be solved using any single matrix inversion scheme. 4.3.3 Scattering parameters After a matrix solution is obtained, the S-parameters of the waveguide section must be obtained. The scattering coefficients along the waveguide cross-sections can be written as S11 = vectorEref vectorEinc vextendsinglevextendsingle vextendsinglevextendsingle vextendsinglez=z 1 = vectorE ? vectorEinc vectorEinc vextendsinglevextendsingle vextendsinglevextendsingle vextendsinglez=z 1 = vectorE vectorEinc vextendsinglevextendsingle vextendsinglevextendsingle vextendsinglez=z 1 ?1 (4.40) S21 = vectorEtrvextendsinglevextendsinglevextendsingle z=z2 vectorEincvextendsinglevextendsinglevextendsingle z=z1 = vectorEvextendsinglevextendsinglevextendsingle z=z2 vectorEincvextendsinglevextendsinglevextendsingle z=z1 (4.41) where the electric fields expressed are the average field over the waveguide cross-section. If only the TE10 mode is assumed to exist, this formulation of the S-parameters is sufficient. However, in the case of the MM formulation, many modes may exist at the boundary. Therefore, the S-parameter definition must be augmented to remove any modes other than the TE10 mode. The S-parameters for this situation are expressed by S11 = integraltext ?2 vectorETE10 ? vectorEd?2 integraltext ?2 vectorETE10 ? vectorEincd?2 ?1 (4.42) 87 S21 = integraltext ?3 vectorETE10 ? vectorEd?3 integraltext ?2 vectorETE10 ? vectorEincd?2 (4.43) 4.3.4 Asymptotic Waveform Evaluation The computation of filter?s response requires the solution of (4.17) over a broad range of frequencies. This can be very time-consuming since the matrix must be inverted at each discrete frequency point. When a filter has rapid variation over a narrow band, more points are needed to accurately sample the response which only compounds the matrix solution problem. In this section, a method for quickly evaluating (4.17) over a broad frequency range is presented. This technique is referred to as Asymptotic Waveform Evaluation [156]. In AWE, the matrix equation is expanded using a Taylor series and is then converted to a rational Pad?e function. The procedure is as follows. Consider a matrix equation of the following form Ax = b (4.44) In general, A, the square matrix, is frequency dependent as well as the unknown x and the source term b. The AWE procedure begins with the expansion of x into a Taylor series [156] as x(f) = Qsummationdisplay n=0 mn (f ?f0)n (4.45) 88 where the mn are unknown coefficients. A and b are then expanded into Taylor series as well. For simplicity, consider only a two term expansion where x(f) = m0 +m1 (f ?f0) A(f) = A(f0)+A?(f0)(f ?f0) b(f) = b(f0)+b?(f0)(f ?f0) (4.46) wheretheprimesindicatederivativeswithrespecttofrequency. Substituting thisinto(4.44) yields A(f0)m0 +A(f0)m1 (f ?f0)+A?(f0)m0 (f ?f0) +A?(f0)m1 (f ?f0)2 = b(f0)+b?(f0)(f ?f0) (4.47) Matching like powers of (f ?f0) then gives the set of equations A(f0)m0 = b(f0) A(f0)m1 +A?(f0)m0 = b?(f0) (4.48) from which the moments of x may be solved recursively as m0 = A?1 (f0)b0 m1 = A?1 (f0)bracketleftbigb?(f0)?A?(f0)m0bracketrightbig (4.49) More generally, the solution of an arbitrary number of moments is given as m0 = A?1 (f0)b0 mn = A?1 (f0) bracketleftBigg b(n) (f0) n! ? nsummationdisplay i=1 A(i) (f0)mi?1 i! bracketrightBigg (4.50) 89 where A(i) is the ith derivative of A and b(n) is the nth derivative of b. The Taylor series expansion has an inherently limited bandwidth. A prohibitive number of terms are necessary to achieve practical bandwidths [123]. Therefore, x can be represented as a well-behaved rational Pad?e function of the form x(f) = Lsummationtext i=0 ai (f ?f0)i 1+ Msummationtext j=1 gj (f ?f0)j (4.51) where L+M = Q and the best performance is obtained when M = L. Substituting (4.51) into (4.45) yields the equations for gj and ai after some manipulation. They are given by ? ?? ?? ?? ?? ?? ?? ? mL mL?1 mL?2 ??? mL?M+1 mL+1 mL mL?1 ??? mL?M+2 mL+2 mL+1 mL ??? mL?M+3 ... ... ... ... ... mL+M?1 mL+M?2 mL+M?3 ??? mL ? ?? ?? ?? ?? ?? ?? ? ? ??? ??? ??? ??? ? ??? ??? ??? ??? ? g1 g2 g3 ... gM ? ??? ??? ??? ??? ? ??? ??? ??? ??? ? = ? ? ??? ??? ??? ??? ? ??? ??? ??? ??? ? mL+1 mL+2 mL+3 ... mL+M ? ??? ??? ??? ??? ? ??? ??? ??? ??? ? (4.52) ai = isummationdisplay j=0 gjmi?j 0 ? i ? L (4.53) Theoretically, more and more terms can be added to (4.45) to achieve higher band- width. However, there are several practical issues to consider. First, AWE is memory intensive since the storage of the derivatives of A is required. Therefore, to maintain reasonable memory limits, the number of terms must usually be kept low. Additionally, solution of the higher order moments is difficult since these equations are ill-conditioned. 90 Therefore, it may be necessary to perform the AWE procedure at multiple points within the frequency band (A must be inverted each time). For this case, a multi-point AWE procedure is needed where a few frequencies are judiciously chosen to maximally utilize the bandwidth obtainable from a given level of expansion. In this work, an 8th order expansion is employed (L = M = 4) as well as the Complex Frequency Hopping [157] to automatically select the AWE frequencies. In other words, this method is a binary search where initially the minimum and maximum frequencies of interest are chosen as expansion points. After evaluation of these AWE solutions, the error between them is checked at a number of points (2 evenly spaced points in this work). If the error is below a certain tolerance, the solution is obtained. If not, the AWE is performed again at the center frequency of the band. Once again, the error is checked at a number of points between [fmin,fcenter] and [fcenter,fmax]. This procedure continues to divide the frequency band as necessary until the error is below the tolerance. It was found that for the problems considered here, around 5 AWE evaluations were required. 4.4 Method Verification To verify the accuracy of each formulation presented, the S-parameters for an X band waveguide loaded with a 1 cm thick slab (?r = 2.2) were calculated over the full X band range. After running several test cases, the ABC formulation was determined to be unsuitable for this study since the size of the computational domain necessary to achieve reasonable accuracy was prohibitively large (number of unknowns). Using the MM formulation, the slab was discretized as shown in Fig. 4.4. For the undoped slab, ?2 and ?3 were chosen to coincide with the faces of the dielectric. The resolution shown by 91 Figure 4.4: Waveguide mesh used for method verification. 92 9 10 11 120 0.1 0.2 0.3 0.4 0.5 0.6 0.7 |S11| Frequency (GHz) Discrete Sweep AWE Sweep Theory Figure 4.5: Comparison of FEM and theoretical |S11| for the loaded waveguide. 93 9 10 11 120.8 0.85 0.9 0.95 1 |S21| Frequency (GHz) Discrete Sweep AWE Sweep Theory Figure 4.6: Comparison of FEM and theoretical |S21| for the loaded waveguide. 94 the mesh is far greater than what is actually necessary but is used anyway for illustrative purposes. The matrix solution (9458 unknowns) was carried out using UMFPACK [158], an efficient direct solver. Figs. 4.5 and 4.6 show the magnitude of S11 and S21, respec- tively, as obtained from theory and the numerical procedure. Two numerical solutions are shown: one which employs the AWE and one which is a traditional discrete sweep. In both cases, the results obtained are indistinguishable from the theoretical result. The discrete sweep consumed 835.08 seconds of CPU time whereas the AWE sweep needed only 84.22 seconds. Therefore, for optimizations presented later, the AWE will be used in all work presented in the remainder of this dissertation. 4.5 Results and Discussion All results discussed in this section refer to a host dielectric material (?r = 2.2, ?r = 1, d = 1 mm) placed in an X band rectangular waveguide (dimensions a = 2.29 cm, b = 1.02 cm) as shown in Fig. 4.1. The longitudinal dimension of waveguide was set to 5 mm with the slab located at the center. The MM scheme was applied to the faces of the waveguide and convergence of the solution was tested for a variety of situations (e.g., undoped and doped slabs). In all cases, the solution was found to converge when the MM scheme allowed up to the first 40 TE modes (and corresponding degenerate TM modes). To verify the accuracy of the FEM, the model was compared to the analytical solu- tion for undoped slabs presented in Chapter 2 and found to be in excellent agreement. The FEM results for doped cases were compared to models created in HFSS, a com- mercially available FEM software tool [159]. For these comparisons, the dopants were 95 manually inserted into the HFSS model. In all cases, agreement with the FEM was found to be excellent as shown in Fig. 4.7. Presently, fabrication techniques do not exist that would allow control of the place- ment and orientation of dopants embedded within a host medium. Therefore, for ran- domly doped materials to be useful they must exhibit nearly identical frequency responses (S-parameters) for a given set of doping conditions (e.g., doping level, feature size, etc.). To determine the consistency of the frequency response when using wire dopants, 20 different random wire distributions were selected for each doping level simulated and their results compared. The discretized geometry detailed previously contained a total of 1921 edges (average edge length 1.25 mm). Doping levels of 50, 100, 150, and 200 wires were each simulated (80 total simulations). The results of these simulations clearly showed that randomly doped dielectric slabs can, for some cases, exhibit notch filter behavior. However, the results lacked any signifi- cant consistency, with number of notches, notch frequency, and level of attenuation being unpredictable. Fig. 4.8 shows two profiles generated using 150 wires and illustrates the extreme variation in the obtained responses. Upon further investigation, the cause of the inconsistency was determined to be the wire selection process. During this procedure, no mechanism was in place to prevent wires from connecting. Consequently, in many cases, random selection of wire locations resulted in the presence of long, erratically shaped wires an example of which is illustrated in Fig. 4.2. Therefore, a second set of simulations was conducted with the algorithm appropriately modified to prohibit wire connection. For this condition, only negligible notch behavior was developed in most cases. How- ever, if a notch did appear, there was noticeable consistency in notch frequency (fnotch 96 9 10 11 12?50 ?40 ?30 ?20 ?10 0 Frequency (GHz) S 21 (dB) HFSS FEM Figure 4.7: Comparison of magnitude of S21 for triangular patch doping using present method and HFSS. 97 = 12.3 GHz) among responses for a given doping level as well as a correlation between doping level and notch frequency. Furthermore, there was no consistency in the level of attenuation for a given doping level as shown in Fig. 4.9. Therefore, since neither the allowance or prevention of wire connections would be controllable during fabrication of such a material, it must be concluded that random doping of a dielectric slab with short wires does not result in a material with any predictable and useful properties. Even though randomly embedded wires in the dielectric showed no promise of pro- ducing a useful material, random wire placement on the dielectric surface has shown some utility. Fabrication of such a material is also relatively easy using printed circuit and microelectronic fabrication techniques. To easily characterize the behavior of surface- doped dielectrics, valid wire locations were limited to vertical (y-directed) edges on the front face of the slab with no wire connections permitted. The limitation of vertical wires on the slab face was motivated by the fact that this particular orientation of wires would have the most dramatic effect on the wave propagation. Doping levels of 10, 20, 35, 50, 65, and 150 wires (selected from Nsurface = 380 total vertical edges) were each simulated 20 times to examine the variation in the frequency responses. In all cases, the reflection and transmission coefficients were found to be nearly identical to the undoped slab?s response. Therefore, for this particular frequency range, 1 mm wires did not result in any form of substantial frequency-dependent behavior. In order to investigate the response?s dependence on wire length, 2 mm wires were generated by connecting two vertically adjacent edges. Again, doping levels of 10, 20, 35, and 50 2 mm wires were each simulated 20 times to determine the utility of such materials. In this case, all 10 wires cases produced responses nearly indistinguishable 98 9 10 11 12?40 ?35 ?30 ?25 ?20 ?15 ?10 ?5 0 Frequency (GHz) S 21 (dB) Figure 4.8: A wide variation of S21 responses is shown for the cases of randomly embed- ded wires. 99 9 10 11 12?20 ?15 ?10 ?5 0 Frequency (GHz) S 21 (dB) Figure 4.9: Illustration of the wide range of notch behaviors observed even when wire connections within the dielectric are not allowed. 100 from the undoped slab?s response. For the 20, 35, and 50 wire cases, a few of the responses exhibited a notch filter behavior as shown in Fig. 4.10. A series of 3 mm wire simulations were also performed which included 20 simulations for each doping level of 10, 20, and 35 wires. Results for the 10 and 20 wire cases generally showed a notch filter response in the 10-12 GHz range (see Fig. 4.11). In the 35 wire simulations, however, the response typically showed the presence of two notches, one between 9-10 GHz and another between 11-12 GHz (see Fig. 4.11). In the final analysis, we can conclude that, for the cases studied, short wires do not yield especially useful reflection and transmission properties for the doped dielectrics. Although, their utility may be questioned, dielectrics with wires only on the surface are certainly realizable using current printed circuit technology. Therefore, the consistency of the random wire configurations? responses no longer presents the same fabrication problem observed in the embedded wire case. However, these conclusions may not be valid for finite wire radius doping. Noting the above behavioral trends seen in the thin wire cases, the response of thin triangular patches (?1 mm edge length) placed on the dielectric surface was also investigated. For these simulations, no restrictions were placed on the possibility of patches connecting to form larger conducting shapes since the fabrication is relatively simple for any arbitrary surface pattern. Twenty simulations were performed for each of the 10, 20, 30, 40, 50, and 75 patch doping levels (120 simulations total). For the case of 10 patches, no significant deviations from the undoped response were observed. At the higher doping levels, a large variation in the notch frequency, bandwidth, level of attenuation, and number of notches was found. Fig. 4.12 shows several of the possible 101 9 10 11 12?30 ?25 ?20 ?15 ?10 ?5 0 Frequency (GHz) S 11 , S 21 (dB) |S21| |S11| Figure 4.10: Magnitude of S11 and S21 for typical case of a dielectric slab doped with 50 2 mm wires. 102 9 10 11 12?50 ?40 ?30 ?20 ?10 0 Frequency (GHz) S 21 (dB) 10 wires 20 wires 35 wires Figure 4.11: Comparison of magnitude of S21 for typical cases of a dielectric slab doped with 10, 20, and 35 3 mm wires. 103 9 10 11 12?60 ?50 ?40 ?30 ?20 ?10 0 Frequency (GHz) S 21 (dB) 75 patches 30 patches 20 patches Figure 4.12: Comparison of magnitude of S21 for cases of a dielectric slab doped with 20, 30, and 75 triangular patches. 104 types of responses realizable using the higher doping levels. As expected, these results indicate that surface patterning based on triangular patches allows a wider variety of notch transmission profiles than thin wires. Therefore, should a certain type of notch behavior be desired, it is necessary only to find the right combination of patches to realize the required response. This problem is addressed in the next chapter. 4.6 Summary An investigation into the possible use of randomly doped dielectrics in order to obtain consistent and useful frequency response profiles has been presented. Using the FEM, the types of responses and feasibility of fabrication of the doped dielectric slab is discussed. The dopants considered include short wires embedded in the dielectric and wires and triangular patches distributed on the surface of the dielectric. Generally, it was found that random doping produces materials with random notch filter behavior. In the case of embedded wires, it is impractical to fabricate such materials due to the sensitivity of the frequency response to wire location and orientation. In other words, randomly doped slabs do not exhibit any predictable or repeatable filter behavior and, therefore, are unsuitable to practical use. However, printed circuit and microelectronic fabrication techniques would easily allow the creation of useful dielectric slabs with wires or patches printed on the surface. In the above studies, it was shown that wires less than 3 mm in length (for X band) are likely to give transmission responses which are too narrowband to be of any practical value. Alternatively, doping with thin triangular patches showed a wide range of notch filter profiles of which many show potential utility. 105 Chapter 5 Waveguide Filter Optimization Using Surface Patches 5.1 Overview ThischapterpresentsanoveloptimizationschemeutilizingahybridMode-Matching/Finite Element Method (MM/FEM) and Genetic Algorithm (GA) to evaluate the versatility of transverse waveguide filter designs. The novelty of the method lies in its completely autonomous framework. By specifying a desired transmission response in terms of scat- tering parameters, the GA optimizes the filter design using the MM/FEM discussed in Chapter 4 and [7,160] without any further input from the user. The power of this method is demonstrated through the optimization of a number of useful X band (8.2 - 12.4 GHz) waveguide filters. Additionally, practical concerns related to the fabrication and ele- mental connectivity within the waveguide are addressed and necessary modifications are discussed. 5.2 Numerical Procedure 5.2.1 Theory The filter geometry used in this work is composed of a dielectric slab located within an X band waveguide as represented in Fig. 5.1. Conducting patches printed on the front face of the slab will perturb the transmission and reflection response (S-parameters) of the slab possibly resulting in novel filter behavior within the frequency band of interest. To characterize the response of such a geometry, a FEM forward solution technique is 106 x y z a b Filter 2 1 3 ? ? ? Figure 5.1: Rectangular waveguide with optimized filter. 107 employed. The details of three-dimensional edge-based FEM formulations are discussed in Chapter 4. A Mode-Matching (MM) scheme is employed on ?2 and ?3 which effectively truncates the FEM domain [6,7,160]. 5.2.2 Optimization Technique The front face of the dielectric slab located within the waveguide was discretized into a regular grid as shown in Fig. 5.2. During the iterative process, sets of pixels are selected to be conductors, and the responses of the resulting structures are determined. Over the course of several iterations, a set of pixels is found which yields the optimal response. Given the discrete nature of the gridding, the GA is an advisable choice for carrying out the optimization. An overview of the GA have been presented in Chapter 2. Fig. 5.3 shows a flowchart of the GA implemented in this chapter. The pixel material is chosen as the optimization parameter such that the front face of the slab is encoded into a binary string where ?1? represents a perfectly electrically conducting (PEC) patch and ?0? represents the dielectric (no conductor). For all the cases presented in Section 5.3, the slab surface is divided into a 20 x 10 grid (x vs. y) resulting in 200 pixels. To reduce the search space, 4-fold symmetry is assumed as illustrated in Figs. 5.1 and 5.2. Therefore, the search space is reduced to a 50 element binary string. Initially, the typical elitism, crossover, and mutation operators are employed in the GA. At each generation, the two best solutions are chosen as elite children and are directly copied into the population of the next generation. Binary tournament selection is used to choose members to be evaluated for crossover [38]. Single-bit and double-bit 108 mutation is implemented as outlined in Chapter 2. The crossover and mutation rates are chosen to be 80% and 10%, respectively [38]. Once again the Total Redundancy Removal GA (TRRGA) scheme is employed to more aggressively search the solution space. The population size is set as 20, and the algorithm terminates after 500 generations if a fitness tolerance (0.01) has not been met. In the examples presented here, |S21| is the desired filter characteristic, and the error function used to evaluate each population member is given by Err = 1N Nsummationdisplay i=1 parenleftBigvextendsinglevextendsingle SFEM21 (fi)vextendsinglevextendsingle? vextendsinglevextendsingle vextendsingleSDesired21 (fi) vextendsinglevextendsingle vextendsingle parenrightBig2 (5.1) where N is the number of frequency points used and vextendsinglevextendsingleSFEM21 (fi)vextendsinglevextendsingle and vextendsinglevextendsingleSDesired21 (fi)vextendsinglevextendsingle are the calculated and desired S21 parameters evaulated at frequency fi, respectively. Once the best population member has a fitness value below a certain value (0.15 in this work), the redundancy removal scheme becomes inefficient. This tolerance level represents an approximate point at which a ?fine-tuning? scheme is better-suited to further reduce the error than aggressive interrogation of the entire solution space. It was found that switching to a high mutation rate (80%) as shown in Fig. 5.3 yielded better performance since the search was then centered around the member with the best fitness. 5.2.3 Practical Considerations Ohira, et al. [72] cite a practical issue associated with FSS optimization. In their work, a MoM code is employed to optimize an FSS element shape. They show that when metal patches are joined only at the corners, the adopted MoM scheme does not allow 109 dx dy b a Connected metallic patches 4-fold symmetry enforced for optimization Figure 5.2: Gridded substrate face showing metallization for fabrication and FEM sim- ulation. 110 Initialize Population Evaluate Fitness Population Generation Elitism Selection Crossover Any redundancies in population? Eliminate extra population Fill population with new members Geometry Refinement Exit Yes No Exit Condition Best fitness below tolerance? Mutation (10%) Mutation (80%) No No Yes Yes Figure 5.3: Flowchart illustrating the GA procedure with geometry refinement, redun- dancy removal, and high mutation procedures. 111 current to flow between the patches. This represents a significant analysis difficulty in terms of assuring that the optimized design would be realizable when fabricated. To correct this deficiency, an automated geometry modification procedure known as the Geometry Refinement Method (GRM) is required. Fig. 5.4 shows two corner-connected patches and two mesh triangles at the patch junction. The traditional FEM vector basis functions suffer from a weakness similar to that described in [72]. The electric field inside the triangular patches on either side of the metallized patches are represented by the set of basis functions. However, when two sides of the triangle have no tangential electric field, only the basis function associated with the outer edge remains to describe the field in this corner region (in the case of zeroth-order basis functions). By definition, this function is zero-valued (both tangential and normal components) at the corner where the PEC patches meet. Therefore, there can be no electrical connection between the PEC patches (no surface charge is present). Similar arguments hold even if the corner junction is discretized in a different fashion such as when only one triangle edge rests along the patch or higher order functions are used. This issue is addressed by the GRM which is part of the GA as shown in Fig. 5.3. After the generation of a new population member, the binary string is scanned for corner-connected patches. When this condition is encountered, one of the two empty pixels is filled (random decision) so as to eliminate the corner connection. The process is repeated until no corner connections are present on the substrate grid. In addition to resolving the problem presented by the basis functions, this method also eliminates any 112 fabrication uncertainty as to whether or not a physical connection was actually made in the measured sample. Similarly, for the algorithm presented here, no PEC patches are allowed along the outer ring of the grid (Fig. 5.4). This ensures that the edge patches are electrically isolated from the waveguide walls. Thus, this requirement further reduces the search space to a 36 bit binary string. 5.2.4 Forward Solution The fitness evaluation of any population member can be carried out via a forward solution using the MM/FEM. This method was found to be an efficient means of deter- mining the S-parameters of the loaded substrate. For additional method verification, an optimization was performed for a certain bandpass filter characteristic and then com- pared to the same geometry simulated in HFSS [159]. Measurements of the fabricated filter were also conducted to ensure that both computational techniques returned accu- rate data. Fig. 5.5 shows the desired response compared to the optimized FEM response as well as the results of HFSS and measurements. As indicated, the FEM scheme led to errors in the computation of the fields. A probable explanation for these errors is the mesh resolution. For the FEM to be a computationally efficient method, the mesh generation was done before optimization so that one mesh, and therefore one matrix could be used throughout the optimization. Also, each pixel was discretized as two tri- angles in the mesh generation scheme. However, in the vicinity of patch corners, the fields vary rapidly and a linear basis functions cannot realize this type of behavior. Ef- forts were made to increase the resolution of the mesh, but the mesh generator available 113 X50X61X74X63X68X65X73X20X6EX6FX74 X61X6CX6CX6FX77X65X64X20X6FX6E X73X75X62X73X74X72X61X74X65X20X65X64X67X65 X43X6FX72X6EX65X72 X63X6FX6EX6EX65X63X74X69X6FX6EX73 X6EX6FX74X20X61X6CX6CX6FX77X65X64 Figure 5.4: The lower left quarter of the waveguide filter is shown illustrating the rules governing the geometry refinement process. The inset illustrates the problem with corner-connected patches. 114 was not sophisticated enough to do this efficiently. The matrix resulting from this high resolution mesh was too computationally intensive to invert and maintain a reasonable optimization time. HFSS, however, returned results that were in excellent agreement with the mea- surements. Therefore, HFSS was integrated into the optimization scheme as the forward solver. The accuracy of HFSS can be attributed to two things: (1) HFSS uses a higher order basis function as opposed to the zeroth-order basis functions used in Chapter 4, and (2) HFSS utilizes an adaptive solver which refines the mesh as necessary in the vicinities of rapid field variation. HFSS employs a MM scheme slightly different from that presented here. Through a hybridization procedure, the field within the waveguide structure and the unknown waveguide eigenfunction amplitudes required for the impedance match at ?2 and ?3 are solved simultaneously. The forward solution begins with a 2.5D eigenvalue solution at each of the ports. Now the field variation for each mode is known at the ports and the unknown amplitudes can be solved for during the FEM solution. Therefore, this has the benefit of directly providing the magnitude and phase of the S-parameters (TE10 mode) at the waveguide ports. 5.3 Results In this section, four design examples are presented to demonstrate the utility of the proposed optimization technique. For all examples, the filter patterns are optimized on a substrate (?r = 2.2, ?r = 1, thickness = 0.7874 mm) placed in a section of WR90 waveguide as represented by Fig. 5.1. The model utilized in HFSS for the optimization 115 9 10 11 120 0.2 0.4 0.6 0.8 1 Frequency (GHz) |S 21 | Ideal FEM Measurement HFSS Figure 5.5: Comparison of ideal data, MM/FEM, HFSS, and measurements. The MM/FEM matches well with the ideal (requested) response. However, the geometry generated during the process yields a different response in actuality as given by HFSS and measurements. 116 procedure is shown in Fig. 5.6. The geometry is generated prior to the optimization, and during the procedure, the algorithm simply selects which patches are to be conductors. As previously mentioned, the GA is allowed to optimize the response until the RMS error between the design goal and simulated responses is less than 0.01 or the algorithm reaches 500 generations. Measurements of fabricated samples are carried out according to the method presented in Section 3.3. 5.3.1 Notch Filter As an initial test, a notch filter was specified with a center frequency of 10.4 GHz and 3 dB bandwidth of 1.4 GHz. For the error function, only five points were evaluated: the center frequency (|S21| = 0), two 3 dB points (|S21| = .5), and the ends of the X band (|S21| = 1). After optimization, the 5-point RMS error between the simulated response and the ideal response was 0.0933. The fabricated sample is shown as pattern A in Fig. 5.7. Measurements showed excellent agreement with the simulation (RMS error = 0.0667) as well as the design goal. All three |S21| responses for this pattern are shown in Fig. 5.8. 5.3.2 Bandpass Filter A bandpass filter was optimized with center frequency of 10.3 GHz and 3 dB band- width of 1.6 GHz. In this case, 51 points were evaluated for the error function. The ideal response was specified as |S21| = exp parenleftBigg ? parenleftbigf ?10.3?109parenrightbig2 1018 parenrightBigg (5.2) 117 Figure 5.6: HFSS model used for optimization 118 Figure 5.7: Fabricated samples produced using the patterns generated by the optimiza- tion procedure. A dime is shown for size reference. 119 X39 X31X30 X31X31 X31X32 X2DX35X30 X2DX34X30 X2DX33X30 X2DX32X30 X2DX31X30 X30 X46X72X65X71X75X65X6EX63X79X20X28X47X48X7AX29 X7CX53 X32X31 X64X42 X7C X49X64X65X61X6C X53X69X6DX75X6CX61X74X69X6FX6E X4DX65X61X73X75X72X65X6DX65X6EX74X73 Figure 5.8: Comparison of S21 responses of ideal, simulated, and measured notch filter. 120 where f is the frequency. The resulting RMS error between the ideal and simulated responses was 0.0482. The filter pattern is noted as B in Fig. 5.7. Measured results were found to be in excellent agreement with the ideal and simulated results as shown in Fig. 5.9. The RMS error between the measurements and simulation was calculated to be 0.0252. 5.3.3 Low Pass Filter The next example optimized was a low pass filter with 3 dB frequency of 11.1 GHz. The ideal response was specified as shown in Fig. 5.10. Optimization yielded pattern C shown in Fig. 5.7 for which the RMS error was 0.0685. Again, agreement was quite good between the ideal, simulated, and measured responses (see Fig. 5.10) with the measurement error being 0.0252. 5.3.4 High Pass Filter Finally, a high pass filter with 9.1 GHz 3 dB frequency was desired as shown in Fig. 5.11. The optimization yielded pattern D in Fig. 5.7 which gave an RMS error of 0.0368. Fig. 5.11 shows the ideal, simulated, and measured responses, all of which show excellent agreement with the exception of a small error near 12 GHz. This small error resulted in a slightly higher measurement error of 0.134 although, clearly, overall agreement is excellent. 5.4 Summary A general optimization scheme for designing transverse rectangular waveguide filters has been presented and shown to yield excellent results for a wide variety of useful filter 121 X39 X31X30 X31X31 X31X32X2DX36X30 X2DX35X30 X2DX34X30 X2DX33X30 X2DX32X30 X2DX31X30 X30 X46X72X65X71X75X65X6EX63X79X20X28X47X48X7AX29 X7CX53 X32X31 X20X64X42 X7C X49X64X65X61X6C X53X69X6DX75X6CX61X74X69X6FX6E X4DX65X61X73X75X72X65X6DX65X6EX74X73 Figure 5.9: Comparison of S21 responses of ideal, simulated, and measured bandpass filter. 122 X39 X31X30 X31X31 X31X32X2DX33X35 X2DX33X30 X2DX32X35 X2DX32X30 X2DX31X35 X2DX31X30 X2DX35 X30 X46X72X65X71X75X65X6EX63X79X20X28X47X48X7AX29 X7CX53 X32X31 X20X64X42 X7C X49X64X65X61X6C X53X69X6DX75X6CX61X74X69X6FX6E X4DX65X61X73X75X72X65X6DX65X6EX74X73 Figure 5.10: Comparison of S21 responses of ideal, simulated, and measured low pass filter. 123 X39 X31X30 X31X31 X31X32X2DX35X30 X2DX34X30 X2DX33X30 X2DX32X30 X2DX31X30 X30 X46X72X65X71X75X65X6EX63X79X20X28X47X48X7AX29 X7CX53 X32X31 X20X64X42 X7C X49X64X65X61X6C X53X69X6DX75X6CX61X74X69X6FX6E X4DX65X61X73X75X72X65X6DX65X6EX74X73 Figure 5.11: Comparison of S21 responses of ideal, simulated, and measured high pass filter. 124 responses. Additionally, this technique can be applied to waveguides of arbitrary cross section while completely alleviating the burden of analytical design techniques. This method employs a modified Genetic Algorithm to aggressively search a large discontinu- ous solution space until a low RMS error has been achieved. At this point the fine-tuning mechanism becomes active. Experiments included optimization of notch, bandpass, low pass, and high pass filters, all of which showed extremely good agreement between ideal, simulated, and measured |S21| responses. 125 Chapter 6 Conclusions With the wireless communications industry and military demand for smaller, higher performance devices, traditional electromagnetic design methods are no longer suitable. A hybridization of computational tools and optimization algorithms has become neces- sary for the design challenges faced by today?s engineers. In response, this dissertation presented several novel methods utilizing fast and efficient optimization techniques that may lead to fully autonomous engineering design and inversion codes in the future. In Chapter 2, a thorough characterization of optimization techniques and error functions was conducted with regard to the extraction of complex permittivities from multilayered dielectric slabs. Sequential Quadratic Programming (SQP), the Genetic Algorithm (GA), and Particle Swarm Optimization (PSO) were applied to ideal and measured S-parameter data, and the complex permittivities of each individual layer were estimated. SQP was found to be a very fast method, but often suffered from local minima trapping. GA and PSO, however, were able to accurately extract the permittivities for multilayer materials (1, 2, and 3 layer materials discussed). After application of the Total Redundancy Removal scheme to the GA, the convergence rate was improved significantly, thereby establishing GA as the preferred method. Of the three error functions considered, the function which utilized only magnitude information achieved superior results in many cases. However, as the number of layers increased, the inclusion of phase information within the error function was determined to be necessary. 126 In Chapter 3, the method of Chapter 2 was extended to include the extraction of complex permeability. To the author?s knowledge, no other efforts to date have been successful in extracting the full set of complex constitutive parameters from individual layers of multilayered slabs. Since this problem results in a much higher dimensionality search space, a novel Multi-Point SQP algorithm was developed and presented to exploit the efficiency and accuracy of SQP while improving its ability to avoid local minima. Overall, the algorithm showed better accuracy than the GA for single layer cases and was over 20 times faster in finding the global minimum. In Chapter 4, a novel methodology for realizing waveguide filters was presented. In this method, small, thin conducting wires (or other conducting shapes) are used to randomly dope a dielectric slab. With proper selection of wire length and density, it was shown that novel filter behaviors can result. However, when practical constraints on the doping methods were applied, the resulting doped slabs showed either no useful filter behavior or erratic behavior not controllable by modern fabrication processes. When conducting patches were used and restricted to the surface of the slab, a variety of novel notch filter characteristics were observed. Using this approach, fabrication difficulties can be avoided since standard printed circuit technology can easily be used to manufacture the filters. To realize specific filter responses, optimization algorithms could be employed to correctly locate the placement of the conducting patches. In Chapter 5, a method is presented to automate the design of waveguide filters using surface patches on dielectric slabs loaded into the guide. This technique utilizes a combination of a FEM forward solver (HFSS) in the context of a GA. The GA presented 127 here was modified with a fine-tuning mechanism which becomes operational once the er- ror has dropped below a certain threshold. Also, some practical fabrication concerns are addressed by employing the Geometry Refinement Method to eliminate corner-connected patches. The novelty and utility of this method were demonstrated by allowing the al- gorithm to optimize patch configurations for a number of filter responses. In all cases, the resulting structure showed excellent agreement with the desired response and was further verified by fabrication and measurements. 6.1 Future Work 6.1.1 Constitutive Parameter Extraction A number of suggestions are given in this section regarding possible avenues for future research regarding the constitutive parameter extraction technique. As previously addressed, the availability of limited information over a band of frequencies determines the number of layers that can be accurately handled. With the addition of more layers and uncertainties in layer thicknesses or number of layers, significantly more information must be provided to limit the number of local minima within the search space. This problem is not only affected by the dimensionality of the search space but by the accuracy of the measured data available. To increase the amount of information available to the extraction algorithm, a free space measurement technique [52] is suggested in which information over both frequency and angle of incidence can be provided. To realize such a scheme, modification of the forward solver is all that is required. With the current availability of numerical electromagnetics codes to handle free space scattering, it is likely 128 that this task could be implemented using FEM, Method of Moments, Finite Different Time Domain, hybrid techniques, etc. 6.1.2 Waveguide Filters Extensions to the waveguide filter design method of Chapter 5 are numerous. The work can easily be extended to account for arbitrary waveguide cross sections since the HFSS geometry is the only parameter requiring modification. Also, the method may be extended to handle multiple layers or to also optimize dielectric thickness. These modifications are trivial but would greatly increase the search space and, therefore, computational intensity of the optimization. Also, it is not recommended, although theoretically possible, to optimize for the dielectric constant of the slab since these materials may not be readily available. 6.1.3 Extension To Antenna Design The methodology utilized throughout this work may also be extended for the design of antennas. Wire antenna shape can be optimized using the techniques presented here to obtain novel radiation patterns or bandwidth characteristics. Similarly, the waveguide filter design method can directly be applied to the automation of patch antenna design. Simply changing the forward solver to an efficient Method of Moments code or modifying the FEM code, the same optimization algorithm can be applied to design novel printed circuit antennas. In conclusion, the present work opens up new areas of research for accurate design, simulation, and characterization of problems in a wide variety of fields. The inversion 129 techniques presented have enormous significance to an array of topics including well- logging, mineral location, unexploded ordnance discrimination, finance, etc. Also, the hybridization of design methodologies and optimization techniques such as those pre- sented here will ultimately lead to fully autonomous engineering tools thereby rapidly accelerating the pace of technological development. 130 Bibliography [1] M. Gen and R. Cheng, Genetic Algorithms and Engineering Optimization. New York: Wiley-Interscience, 1999. [2] C. Bennett, ?Time domain inverse scattering,? IEEE Transactions on Antennas and Propagation, vol. AP-29, no. 2, pp. 213?219, Mar 1981. [3] K. Preis, O. Biro, M. Friedrich, A. Gottvald, and C. Magele, ?Comparison of different optimization strategies in the design of electromagnetic devices,? IEEE Transactions on Magnetics, vol. 27, no. 5, pp. 4154?4157, Sep 1991. [4] J. M. Johnson, Y. Rahmat-Samii, and M. Manteghi, ?Genetic algorithms and method of moments (GA/MOM) for the design of integrated antennas,? IEEE Transactions on Antennas and Propagation, vol. 47, no. 10, pp. 1606?1614, Oct 1999. [5] F. J. Villegas, T. Cwik, Y. Rahmat-Samii, and M. Manteghi, ?A parallel elec- tromagnetic genetic-algorithm optimization (EGO) application for patch antenna design,? IEEE Transactions on Antennas and Propagation, vol. 52, no. 9, pp. 2424?2435, Sep 2004. [6] K. Ise, K. Inoue, and M. Koshiba, ?Three-dimensional finite element solution of dielectric scattering obstacles in a rectangular waveguide,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-38, no. 9, pp. 1352?1359, Sept 1990. [7] Z. Cendes and J. Lee, ?The transfinite element method for modeling MMIC de- vices,? IEEE Transactions on Microwave Theory and Techniques, vol. 36, no. 12, pp. 1639?1649, Dec 1988. [8] T. G. Kolda, R. M. Lewis, and V. Torczon, ?Optimization by direct search: New perspectives on some classical and modern methods,? SIAM Review, vol. 45, no. 3, pp. 385?482, 2003. [9] R. M. Lewis, V. Torczon, and M. W. Trosset, ?Direct search methods: Then and now,? Journal of Computational and Applied Mathematics, vol. 124, pp. 191?207, 2000. [10] W. H. Swann, Direct search methods. London: Academic Press, 1972, pp. 13?28. [11] G. Berman, ?Minimization by successive approximation,? SIAM Journal of Nu- merical Analysis, vol. 3, pp. 123?133, 1966. 131 [12] ??, ?Lattice approximations to the minima of functions of several variables,? Journal of the ACM, vol. 16, pp. 286?294, 1969. [13] E. Polak, Computational Methods in Optimization: A Unified Approach. New York: Academic Press, 1971. [14] V. Torczon, ?On the convergence of the multidirectional search algorithm,? SIAM Journal on Optimization, vol. 1, p. 123145, 1991. [15] P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization. San Diego, CA, USA: Academic Press, 1981. [16] J. J. More, ?The levenberg-marquardt algorithm: Implementation and theory,? in Numerical Analysis, ser. Lecture Notes in Mathematics 630, G. A. Watson, Ed. Berlin, Germany: Springer-Verlag, 1977, pp. 105?116. [17] M. C. Biggs, ?Constrained minimization using recursive quadratic programming,? in Towards Global Optimization, L. C. W. Dixon and G. P. Szergo, Eds. Amster- dam, Netherlands: North-Holland, 1975, pp. 341?349. [18] S. P. Han, ?A globally convergent method for nonlinear programming,? Journal of Optimization Theory and Applications, vol. 22, no. 3, pp. 297?309, July 1977. [19] W. Hock and K. Schittowski, ?A comparative performance evaluation of 27 non- linear programming codes,? Computing, vol. 30, no. 4, pp. 335?354, Apr 1983. [20] K. Schittkowski, ?Nlqpl: A fortran-subroutine solving constrained nonlinear programming problems,? Annals of Operations Research, vol. 5, pp. 485?500, 1985/1986. [21] R. Fletcher, Practical Methods of Optimization, 2nd ed. New York, USA: John Wiley and Sons, 2000. [22] A. E. Eiben and J. E. Smith, Introduction to Evolutionary Computing. Berlin: Springer-Verlag, 2003. [23] G. L. Nemhauser and L. A. Wolsey, Integer and Combinatorial Optimization. New York: John Wiley & Sons, 1988. [24] R. W. Becker and G. V. Lago, ?A global optimization algorithm,? in Proceedings of the 8th Allerton Conference on Circuits and Systems Theory, 1970, pp. 3?12. [25] J. H. Holland, ?Outline for a logical theory of adaptive systems,? Journal of the Association for Computing Machinery, vol. 3, p. 297314, 1962. [26] J. Mockus, ?Application of bayesian approach to numerical methods of global and stochastic optimization,? Journal of Global Optimization, vol. 4, no. 4, pp. 347?356, 1994. 132 [27] G. L. Bilbro and W. E. Snyder, ?Optimization of functions with many minima,? IEEE Transactions on Systems, Man, and Cybernetics, vol. 21, no. 4, pp. 840?849, 1991. [28] D. E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learn- ing. New York, USA: Addison-Wesley, 1989. [29] R. C. Eberhart and J. Kennedy, ?A new optimizer using particle swarm theory,? in Proceedings of the Sixth International Symposium on Micromachine and Human Science, 1995, pp. 39?43. [30] R. Storn and K. Price, ?Differential evolution - a simple and efficient adaptive scheme for global optimization over continuous spaces,? ICSI, Technical Report TR-95-012, 1995. [31] A. S. Fraser, ?Simulation of genetic systems by automatic digital computers,? Australian Journal of Biological Sciences, vol. 10, pp. 484?491, 1957. [32] N. A. Barricelli, ?Symbiogenetic evolution pocesses realized by artificial methods,? Methodos, vol. 9, no. 35-36, pp. 143?184, 1957. [33] G. E. P. Box, ?A method for increasing industrial productivity,? Journal of the Royal Statistical Society C, vol. 6, no. 2, pp. 81?101, 1957. [34] Y. Rahmat-Samii and E. Michielssen, Eds., Electromagnetic Optimization by Ge- netic Algorithms. New York, USA: Wiley-Interscience, 1999. [35] A. Boag, E. Michielssen, and R. Mittra, ?Design of electrically loaded wire anten- nas using genetic algorithms,? IEEE Transactions on Antennas and Propagation, vol. 44, no. 5, pp. 687?695, May 1996. [36] R. L. Haupt, ?Thinned arrays using genetic algorithms,? IEEE Transactions on Antennas and Propagation, vol. 42, no. 7, pp. 993 ? 999, July 1994. [37] E. Michielssen, J.-M. Sajer, S. Ranjithan, and R. Mittra, ?Design of lightweight, broad-band microwave absorbers using genetic algorithms,? IEEE Transactions on Microwave Theory and Techniques, vol. 41, no. 6, pp. 1024 ? 1031, Jun 1993. [38] J. M. Johnson and Y. Rahmat-Samii, ?Genetic algorithms in engineering electro- magnetics,? IEEE Antennas and Propagation Magazine, vol. 39, no. 4, pp. 7?21, Aug 1997. [39] R. L. Haupt, ?An introduction to genetic algorithms for electromagnetics,? IEEE Antennas and Propagation Magazine, vol. 37, no. 2, pp. 7?15, Apr 1995. 133 [40] Y. Shi and R. Eberhart, ?A modified particle swarm optimizer,? in Evolutionary Computation Proceedings, IEEE World Congress on Computational Intelligence, May 1998, pp. 69 ? 73. [41] J. Kennedy and R. Mendes, ?Neighborhood topologies in fully-informed and best- of-neighborhood particle swarms,? in Proceedings of the 2003 IEEE International Workshop on Soft Computing in Industrial Applications 2003 (SMCia/03), 2003, pp. 45?50. [42] J. Robinson and y. Rahmat-Samii, ?Particle swarm optimization in electromagnet- ics,? IEEE Transactions on Antennas and Propagation, vol. 52, no. 2, pp. 397 ? 407, Feb 2004. [43] D. W. Boeringer and D. H. Werner, ?Particle swarm optimization versus genetic algorithms for phased array synthesis,? IEEE Transactions on Antennas and Prop- agation, vol. 52, no. 3, pp. 771 ? 779, Mar 2004. [44] J. Vesterstrom and R. Thomsen, ?Particle swarm optimization versus genetic al- gorithms for phased array synthesis,? vol. 52, no. 3, Mar 2004, pp. 771 ? 779. [45] D. G. Kurup, M. Himdi, and A. Rydberg, ?Synthesis of uniform amplitude un- equally spaced antenna arrays using the differential evolution algorithm,? IEEE Transactions on Antennas and Propagation, vol. 51, no. 9, pp. 2210?2217, Sep 2003. [46] A. Qing, ?Electromagnetic inverse scattering of multiple two-dimensional perfectly conducting objects by the differential evolution strategy,? IEEE Transactions on Antennas and Propagation, vol. 51, no. 6, pp. 1251?1262, Jun 2003. [47] K. Gupta and P. S. Hall, Analysis and Design of Integrated Circuit-Antenna Mod- ules. New York, USA: John Wiley and Sons, 1999. [48] S. Chakravarty and R. Mittra, ?Application of the micro-genetic algorithm to the design of spatial filters with frequency-selective surfaces embedded in dielectric media,? IEEE Transactions on Electromagnetic Compatibility, vol. 44, no. 2, pp. 338?346, May 2002. [49] S. Chakravarthy, R. Mittra, and N. Williams, ?Application of the micro-genetic algorithm to the design of spatial filters with frequency-selective surfaces embedded in dielectric media,? IEEE Transactions on Antennas and Propagation, vol. 50, no. 3, pp. 284?296, Mar 2002. [50] M. D. Deshpande and K. Dudley, ?Estimation of complex permittivity of compos- ite multilayer material at microwave frequency using waveguide measurements,? NASA Langley Research Center, Hampton, VA, USA, NASA Technical Memoran- dum 212398, 2003. 134 [51] C.-W. Chang, K.-M. Chen, and J. Qian, ?Nondestructive determination of electro- magnetic parameters of dielectric materialsat x-band frequenciesusing a waveguide probe system,? IEEE Transactions on Instrumentation and Measurement, vol. 46, no. 5, pp. 1084?1092, Oct 1997. [52] T. Zwick, J. Haala, and W. Wiesbeck, ?A genetic algorithm for the evaluation of material parameters of compound multilayered structures,? IEEE Transactions on Microwave Theory and Techniques, vol. 50, no. 4, pp. 1180?1187, Apr 2002. [53] M. D. Janezic and J. Baker-Jarvis, ?Full-wave analysis of a split-cylinder resonator for nondestructive permittivity measurements,? IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 10, pp. 2014?2020, Oct 1999. [54] P. Queffelec and P. Gelin, ?New method for determining the permeability tensor of magnetized ferrites in a wide frequency range,? IEEE Transactions on Instru- mentation and Measurement, vol. 48, no. 8, pp. 519?522, Aug 2000. [55] C. J. Reddy, M. D. Deshpande, and G. A. Hanidu, ?Learning tool for estima- tion of complex permittivity of dielectric material,? in Proceedings of Frontiers in Education Conference, Nov 1996, pp. 1195 ?1196. [56] M. D. Janezic and J. A. Jargon, ?Complex permittivity determination from propa- gation constant measurements,? IEEE Microwave and Guided Wave Letters, vol. 9, no. 2, pp. 76?78, Feb 1999. [57] B. A. Sanadiki and M. Mostafavi, ?Inversion of inhomogeneous continuously vary- ing dielectric profiles using open-ended waveguides,? IEEE Transactions on An- tennas and Propagation, vol. 39, no. 2, pp. 158?163, Feb 1991. [58] Y. Tajima and Y. Sawayama, ?Design and analysis of a waveguide-sandwich mi- crowave filter,? IEEE Transactions on Microwave Theory and Techniques, vol. 22, no. 9, pp. 839?841, Sep 1974. [59] R. Vahldieck and W. J. R. Hoefer, ?Finline and metal insert filters with improved passband separation and increased stopband attenuation,? IEEE Transactions on Microwave Theory and Techniques, vol. 33, no. 12, pp. 1333?1339, Dec 1985. [60] M. Guglielmi, ?Simple CAD procedure for microwave filters and multiplexers,? IEEE Transactions on Microwave Theory and Techniques, vol. 42, no. 7, p. 13471352, Jul 1994. [61] R. J. Langley, ?A dual-frequency band waveguide using FSS,? IEEE Microwave and Guided Wave Letters, vol. 3, no. 1, pp. 9?10, Jan 1993. 135 [62] F. Arndt, R. Beyer, J. M. Reiter, T. Sieverding, and T. Wolf, ?Automated design of waveguide components using hybrid mode-matching/numerical EM building- blocks in optimization-oriented CAD frameworks - State-of-the-Art and Recent Advances,? IEEE Transactions on Microwave Theory and Techniques, vol. 45, no. 5, pp. 747?760, May 1997. [63] T. Shen, H. T. Hsu, K. A. Zaki, A. E. Atia, and T. G. Dolan, ?Full-wave design of canonical waveguide filters by optimization,? IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 2, pp. 504?511, Feb 2003. [64] D. Budamir and G. Goussetis, ?Design of asymmetrical RF and microwave band- pass filters by computer optimization,? IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 4, pp. 1174?1178, Apr 2003. [65] F. Arndt, J. Brandt, V. Catina, J. Ritter, I. Rullhusen, J. Dauelsberg, U. Hilge- fort, and W. Wessel, ?Fast CAD and optimization of waveguide components and aperture antennas by hybrid mm/fe/mom/fd methods - State-of-the-Art and Re- cent Advances,? IEEE Transactions on Microwave Theory and Techniques, vol. 52, no. 1, pp. 292?305, Jan 2004. [66] S. Bila, D. Baillargeat, M. Aubourg, S. Verdeyme, F. Seyfert, L. Barachart, C. Boi- chon, F. Thevenon, J. Puech, C. Zanchi, L. Lapierre, and J. Sombrin, ?Finite- element modeling for the design optimization of microwave filters,? IEEE Trans- actions on Magnetics, vol. 40, no. 2, pp. 1472?1475, Mar 2004. [67] J. V. M. Ros, P. S. Pacheco, H. E. Gonzalez, V. E. B. Esbert, C. B. Martin, M. T. Calduch, S. C. Borras, and B. G. Martinez, ?Fast automated design of waveguide filters using aggressive space mapping with a new segmentation strategy and a hybrid optimization algorithm,? IEEE Transactions on Microwave Theory and Techniques, vol. 53, no. 4, pp. 1130?1142, Apr 2005. [68] D. S. Lockyer and J. C. Vardaxoglou, ?Reconfigurable FSS response from two layers of slotted dipole arrays,? Electronics Letters, vol. 32, no. 6, pp. 512?513, Mar 1996. [69] R. D. Seager, J. C. Vardaxoglou, and D. S. Lockyer, ?Close coupled resonant aper- ture inserts for waveguide filtering applications,? IEEE Microwave and Wireless Components Letters, vol. 11, no. 3, pp. 112?114, Mar 2001. [70] A. Monorchio, G. Manara, U. Serra, G. Marola, and E. Pagana, ?Design of waveguide filters by using genetically optimized frequency selective surfaces,? IEEE Microwave and Wireless Components Letters, vol. 15, no. 6, pp. 407?409, June 2005. 136 [71] G. Manara, A. Monorchio, and R. Mittra, ?Frequency selective surface design based on genetic algorithm,? Electronics Letters, vol. 35, no. 17, pp. 1400?1401, Aug 1999. [72] M. Ohira, H. Deguchi, M. Tsuji, and H. Shigesawa, ?Multiband single-layer frequency selective surface designed by combination of genetic algorithm and geometry-refinement technique,? IEEE Transactions on Antennas and Propaga- tion, vol. 52, no. 11, pp. 2925? 2931, Nov 2004. [73] E. A. Parker, A. D. Chuprin, J. C. Batchelor, and S. B. Savia, ?Ga optimisation of crossed dipole fss array geometry,? Electronics Letters, vol. 37, no. 16, pp. 996?997, Aug 2001. [74] S. Chakravarty and R. Mittra, ?Design of a frequency selective surface (fss) with very low cross-polarization discrimination via the parallel micro-genetic algorithm (pmga),? IEEE Transactions on Antennas and Propagation, vol. 51, no. 7, pp. 1664? 1668, July 2003. [75] M. Bozzi, G. Manara, A. Monorchio, and L. Perregrini, ?Automatic design of inductive fsss using the genetic algorithm and the mom/bi-rme analysis,? Antennas and Wireless Propagation Letters, vol. 1, pp. 91? 93, 2002. [76] S. Chakravarty and R. Mittra, ?Application of the micro-genetic algorithm to the design of spatial filters with frequency-selective surfaces embedded in dielectric media,? IEEE Transactions on Electromagnetic Compatibility, vol. 44, no. 2, pp. 338?346, May 2002. [77] S. Chakravarty, R. Mittra, and N. R. Williams, ?Application of a microgenetic algorithm (mga) to the design of broadband microwave absorbers using multiple frequency selective surface screens buried in dielectrics,? IEEE Transactions on Antennas and Propagation, vol. 50, no. 3, pp. 284?296, Mar 2002. [78] Z. Li, P. Y. Papalambros, and J. L. Volakis, ?Frequency selective surface design by integrating optimisation algorithms with fast full wave numerical methods,? IEE Proceedings on Microwaves, Antennas and Propagation, vol. 149, no. 3, pp. 175?180, Jun 2002. [79] Agilent 85071E Materials Measurement Software: Technical Overview, Agilent, Palo Alto, CA, USA, 2003. [80] Product note on Hewlett-Packard 85071 Materials Measurement Software ?The Measurement of Both Permittivity and Permeability of Solid Materials?, Hewlett- Packard, Palo Alto, CA, USA. [81] D. Pozar, Microwave Engineering, 3rd ed. New York, USA: Wiley, 2003. 137 [82] R. Ludwig and P. Bretchko, RF Circuit Design: Theory and Applications. Upper Saddle River, NJ, USA: Prentice-Hall, 2000. [83] M. J. D. Powell, A fast algorithm for nonlinearly constrained optimization calcu- lations, G. A. Watson, Ed. Berlin: Springer-Verlag, 1977. [84] ??, ?Algorithms for nonlinear constraints that use langrangian functions,? Math- ematical Programming, vol. 14, no. 1, pp. 224?248, 1978. [85] ??, The convergence of variable metric methods for nonlinearly constrainted op- timization calculations, O. Mangasarian, R. Meyer, and S. Robinson, Eds. New York: Academic Press, 1978. [86] U. M. Garcia-Palomares and O. L. Mangasarian, ?Superlinearly convergent quasi- Newton,? Mathematical Programming, vol. 11, pp. 1?13. [87] S. P. Han, ?Superlinearly convergent variable metric algorithms for general non- linear programming problems,? Mathematical Programming, vol. 11, pp. 263?282, 1976. [88] MATLAB Documentation, Mathworks, Inc., 2005. [Online]. Available: http://www.mathworks.com [89] P. E. Gill, W. Murray, and M. H. Wright, Numerical Linear Algebra and Opti- mization. New York: Addison Wesley, 1991, vol. 1. [90] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright, ?Procedures for opti- mization problems with a mixture of bounds and general linear constraints,? ACM Transactions on Mathematical Software, vol. 10, pp. 282?298, 1984. [91] J. T. Alander, ?Indexed bibliography of genetic algorithms in biosciences,? Univer- sity of Vaasa, Department of Information Technology and Production Economics, Report 94-1-BIO, 1995. [92] ??, ?Indexed bibliography of genetic algorithms and neural networks,? Univer- sity of Vaasa, Department of Information Technology and Production Economics, Report 94-1-NN, 1995. [93] ??, ?Indexed bibliography of genetic algorithms in economics,? University of Vaasa, Department of Information Technology and Production Economics, Report 94-1-ECO, 1995. [94] ??, ?A bibliography collection of evolutionary algorithms,? University of Vaasa, Department of Information Technology and Production Economics, Report 94-1-*, 1997. 138 [95] X. Hu, ?PSO bibliography.? [Online]. Available: http://www.swarmintelligence.org/bibliography.php?sort=2 [96] B. Wolfson and S. Wentworth, ?Complex permittivity and permeability measure- ment using rectangular waveguide,? Microwave and Optical Technology Letters, vol. 27, no. 3, pp. 180?182, Nov 2000. [97] ??, ?Complex permittivity and permeability measurement at elevated temper- atures using rectangular waveguide,? Microwave and Optical Technology Letters, vol. 38, no. 6, pp. 449?553, Sep 2003. [98] E. . Cuming. [Online]. Available: http://www.eccosorb.com [99] B. Wolfson, ?Complex permittivity and permeability extraction at elevated tem- peratures using rectangular waveguide,? Master?s thesis, Auburn University, Al- abama, 2001. [100] Introduction to Basic Measurements Using the HP 8510, Hewlett-Packard, Palo Alto, CA, USA, 1984. [101] V. M. Deshpande and M. D. Deshpande, ?Study of electromagnetic wave propa- gation through dielectric slab doped randomly with thin metallic wires using finite element method,? IEEE Microwave and Wireless Components Letters. [102] M. Koshiba and M. Suzuki, ?Finite-element analysis of h-plane waveguide junction with arbitrarily shaped ferrite post,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-34, no. 1, pp. 103?109, Jan 1986. [103] K. Ise, K. Inoue, and M. Koshiba, ?Three-dimensional finite-element method with edge elements for electromagnetic waveguide discontinuities,? IEEE Transactions on Microwave Theory and Techniques, vol. 39, no. 8, pp. 1289?1295, Aug 1991. [104] J. M. Jin and S. W. Lee, ?Hybrid finite element analysis of scattering and radiation by a class of waveguide-fed structures,? Microwave and Optical Technology Letters, vol. 7, no. 17, pp. 798?803, 1994. [105] K. Hirayama, S. Alam, Y. Hayashi, and M. Koshiba, ?Vector finite element method with mixed-interpolation-type triangular-prism element for waveguide discontinu- ities,? IEEE Transactions on Microwave Theory and Techniques, vol. 42, no. 12, pp. 2311?2316, Dec 1994. [106] M. D. Deshpande and C. J. Reddy, ?Application of fem to estimate complex per- mittivity of dielectric material at microwave frequency using waveguide measure- ments,? NASA Langley Research Center, Hampton, OH, USA, NASA Contractors Report CR-198203, 1995. 139 [107] H. B. Lee and T. Itoh, ?A systematic optimum design of waveguide-to-microstrip transition,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT- 45, no. 5, pp. 803?809, May 1997. [108] J. Rubio, J. Arroyo, and J. Zapata, ?Analysis of passive microwave circuits by using a hybrid 2-d and 3-d finite-element mode-matching method,? IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 9, pp. 1746?1749, Sept 1999. [109] K. L. Wu and H. Wang, ?A rigorous modal analysis of h-plane waveguide t-junction loaded with a partial-height post for wide-band applications,? IEEE Transactions on Microwave Theory and Techniques, vol. 49, no. 5, pp. 893?901, May 2001. [110] C.-W. Chang, K.-M. Chen, and J. Qian, ?Nondestructive measurements of complex tensor permittivity of anisotropic materials using a waveguide probe system,? IEEE Transactions on Microwave Theory and Techniques, vol. 44, no. 7, pp. 1081?1090, Jul 1996. [111] A. I. Khalil and M. B. Steer, ?A generalized scattering matrix method using the method of moments for electromagnetic analysis of multilayered structures in waveguide,? IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 11, pp. 2151?2157, Nov 1999. [112] A. I. Khalil, A. B. Yakovlev, and M. B. Steer, ?Efficient method-of-moments formu- lation for the modeling of planar conductive layers in a shielded guided-wave struc- ture,? IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 9, pp. 1730?1736, Sep 1999. [113] R. Bunger and F. Arndt, ?Moment-method analysis of arbitrary 3-d metallic n-port waveguide structures,? IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 4, pp. 531?537, Apr 2000. [114] F. Moglie, T. Rozzi, P. Marcozzi, and A. Schiavoni, ?A new termination condition for the application of fdtd techniques to discontinuity problems in close homoge- neous waveguide,? IEEE Microwave and Guided Wave Letters, vol. 2, no. 12, pp. 475?477, Dec 1992. [115] J. V. Hese and D. D. Zutter, ?Modeling of discontinuities in general coaxial waveguide structures by the fdtd-method,? IEEE Transactions on Microwave The- ory and Techniques, vol. 40, no. 3, pp. 547?556, Mar 1992. [116] J. F. Lee, R. Palandech, and R. Mittra, ?Modeling three-dimensional discontinu- ities in waveguides using nonorthogonal fdtd algorithm,? IEEE Transactions on Microwave Theory and Techniques, vol. 40, no. 2, pp. 346?352, Feb 1992. 140 [117] F. Moglie, T. Rozzi, and P. Marcozzi, ?Wideband matching of waveguide discon- tinuities by fdtd methods,? IEEE Transactions on Microwave Theory and Tech- niques, vol. 42, no. 11, pp. 2093?2098, Nov 1994. [118] J. C. Olivier and D. A. McNamara, ?Analysis of multiport discontinuities in waveguide using a pulsed fdtd approach,? IEEE Transactions on Microwave The- ory and Techniques, vol. 42, no. 12, pp. 2229?2238, Dec 1994. [119] L. A. Vielva, J. A. Pereda, A. Vegas, and A. Prieto, ?Simulating 3d waveguide discontinuities using a combination of prony?s method and fdtd with improved absorbing boundary conditions,? IEE Proceedings on Microwaves, Antennas and Propagation, vol. 141, no. 2, pp. 127?132, Apr 1994. [120] A. P. Zhao and A. V.Raisanen, ?Application of a simple and efficient source exci- tation technique to the fdtd analysis of waveguide and microstrip circuits,? IEEE Transactions on Microwave Theory and Techniques, vol. 44, no. 9, pp. 1535?1539, Sep 1996. [121] A. Kreczkowski, T. Rutkowski, and M. Mrozowski, ?Fast modal abc?s in the hybrid pee-fdtd analysis of waveguide discontinuities,? IEEE Microwave and Guided Wave Letters, vol. 9, no. 5, pp. 186?188, May 1999. [122] S. Wang and F. L. Teixeira, ?An equivalent electric field source for wideband fdtd simulations of waveguide discontinuities,? IEEE Microwave and Wireless Compo- nents Letters, vol. 13, no. 1, pp. 27?29, Jan 2003. [123] J. M. Jin, The Finite Element Method in Electromagnetics, 2nd ed. New York, USA: Wiley-Interscience, 2002. [124] R. F. Harrington, Time-Harmonic Electromagnetic Fields. New York, USA: McGraw-Hill, 1961. [125] J. Meixner, ?The behavior of electromagnetic fields at edges,? IEEE Transactions on Antennas and Propagation, vol. AP-20, no. 7, pp. 442?446, July 1972. [126] J. B. Andersen and V. V. Solodukhov, ?Field behavior near a dielectric wedge,? IEEE Transactions on Antennas and Propagation, vol. AP-26, no. 7, pp. 598?602, July 1978. [127] J. V. Bladel, ?Field singularities at metal-dielectric wedges,? IEEE Transactions on Antennas and Propagation, vol. AP-33, no. 4, pp. 450?455, Apr 1985. [128] R. D. Smedt and J. V. Bladel, ?Field singularities at the tip of a metallic cone of arbitrary cross-section,? IEEE Transactions on Antennas and Propagation, vol. AP-34, no. 7, pp. 865?870, July 1986. 141 [129] C. W. Crowley, P. P. Silvester, and H. Hurwitz, ?Covariant projection elements for 3d vector field problems,? IEEE Transactions on Magnetics, vol. MAG-24, pp. 397?400, 1988. [130] K. D. Paulsen and D. R. Lynch, ?Elimination of vector parasites in finite element maxwell solutions,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-39, no. 3, pp. 395?404, Mar 1991. [131] H. Whitney, Geometric Integration Theory. Princeton, NJ: Princeton University Press, 1957. [132] J. C. Nedelec, ?Mixed finite elements in r3,? Numerical Methods, vol. 35, pp. 315? 341, 1980. [133] A. Bossavit and J. C. Verite, ?A mixed fem-biem method to solve 3-d eddy current problems,? IEEE Transactions on Magnetics, vol. MAG-18, no. 3, pp. 431?435, Mar 1982. [134] M. Hano, ?Finite-element analysis of dielectric-loaded waveguides,? IEEE Trans- actions on Microwave Theory and Techniques, vol. MTT-32, no. 10, pp. 1275?1279, Oct 1984. [135] G. Mur and A. T. de Hoop, ?A finite-element method for computing three- dimensional electromagnetic fields in inhomogeneous media,? IEEE Transactions on Magnetics, vol. MAG-21, no. 11, pp. 2188?2191, Nov 1985. [136] J. S. van Welij, ?Calculation of eddy current in terms of h on hexahedra,? IEEE Transactions on Magnetics, vol. MAG-21, no. 11, pp. 2239?2241, Nov 1985. [137] M. L. Barton and Z. J. Cendes, ?New vector finite elements for three-dimensional magnetic field computation,? Journal of Applied Physics, vol. 61, no. 8, pp. 3919? 3921, Apr 1987. [138] M. Hano, ?Vector finite-element solution of anisotropic waveguides using novel triangular elements,? Electron. Commun. Jpn. Pt. 2, vol. 71, no. 8, pp. 71?80, 1988. [139] ??, ?Finite-element solution of three-dimensional resonator problems: Novel rec- tangular parallelepiped elements,? Electron. Commun. Jpn. Pt. 2, vol. 71, no. 7, pp. 27?34, 1988. [140] A. Bossavit, ?A rationale for ?edge-elements? in 3-d field computations,? IEEE Transactions on Magnetics, vol. MAG-24, no. 1, pp. 74?79, Jan 1988. [141] D. R. Tanner and A. F. Peterson, ?Vector expansion functions for the numerical solution of maxwell?s equations,? Microwave and Optical Technology Letters, vol. 2, no. 2, pp. 331?334, 1989. 142 [142] A. Bossavit and I. Mayergoyz, ?Edge-elements for scattering problems,? IEEE Transactions on Magnetics, vol. MAG-25, no. 7, pp. 2816?2821, July 1989. [143] A. Kameari, ?Calculation of transient 3d eddy current using edge-elements,? IEEE Transactions on Magnetics, vol. MAG-26, no. 3, pp. 466?469, Mar 1990. [144] X. Yuan, ?Three-dimensional electromagnetic scattering from inhomogeneous ob- jects by the hybrid moment and finite element method,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-38, no. 8, pp. 1053?1058, Aug 1990. [145] K. Sakiyama, H. Kotera, and A. Ahagon, ?3-d electromagnetic field mode analysis using finite element method by edge element,? IEEE Transactions on Magnetics, vol. MAG-26, no. 9, pp. 1759?1761, Sept 1990. [146] J. D?Angelo and I. D. Mayergoyz, ?Three dimensional rf scattering by the finite element method,? IEEE Transactions on Magnetics, vol. MAG-27, no. 9, pp. 3827? 3832, Sept 1991. [147] J. M. Jin and J. L. Volakis, ?Electromagnetic scattering by and transmission through a three-dimensional slot in a thick conducting plane,? IEEE Transactions on Antennas and Propagation, vol. AP-39, no. 4, pp. 543?550, Apr 1991. [148] J. S. Wang and N. Ida, ?Eigenvalue analysis in electromagnetic cavities using divergence free finite elements,? IEEE Transactions on Magnetics, vol. MAG-27, no. 9, pp. 3978?3981, Sept 1991. [149] Z. J. Cendes, ?Vector finite elements for electromagnetic field,? IEEE Transactions on Magnetics, vol. MAG-27, no. 9, pp. 3958?3966, Sept 1991. [150] J. F. Lee, D. K. Sun, and Z. J. Cendes, ?Tangential vector finite elements for electromagnetic field computation,? IEEE Transactions on Magnetics, vol. MAG- 27, no. 9, pp. 4032?4035, Sept 1991. [151] ??, ?Full-wave analysis of dielectric waveguides using tangential vector finite elements,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT- 39, no. 8, pp. 1262?1271, Aug 1991. [152] R. Miniowitz and J. P. Webb, ?Covariant-projection quadrilateral elements for the analysis of waveguides with sharp edges,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-39, no. 3, pp. 501?505, Mar 1991. [153] ??, ?Analysis of 3-d microwave resonators using covariant-projection elements,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-39, no. 11, pp. 1895?1899, Nov 1991. 143 [154] A. Chatterjee, J. M. Jin, and J. L. Volakis, ?Robust finite element solution for three-dimensional scattering,? Electronics Letters, vol. 28, no. 10, pp. 966?967, May 1992. [155] ??, ?Computation of cavity resonances using edge-based finite elements,? IEEE Transactions on Microwave Theory and Techniques, vol. MTT-40, no. 11, pp. 2106? 2108, Nov 1992. [156] L. T. Pillage and R. A. Rohrer, ?Asymptotic waveform evaluation for timing analy- sis,? IEEE Transactions on Computer-Aided Design, vol. CAD-9, no. 4, pp. 352? 366, April 1990. [157] M. A. Kolbehdari, M. Srinivasan, M. S. Nakhla, Q. J. Zhang, and R. Achar, ?Simul- taneous time and frequency domain solutions of EM problems using finite element and CFH techniques,? IEEE Transactions on Microwave Theory and Techniques, vol. 44, no. 9, pp. 1526?1534, Sep 1996. [158] T. A. Davis, ?Algorithm 832: UMFPACK - an unsymmetric-pattern multifrontal method with a column pre-ordering strategy,? ACM Transactions of Mathematical Software, vol. 30, no. 2, pp. 196?199, 2004. [159] HFSS Online Help Manual, Ansoft Corporation, Pittsburgh, PA, USA, 2004. [160] J. Bracken, D. Sun, and Z. Cendes, ?S-Domain methods for simultaneous time and frequency characterization of electromagnetic devices,? IEEE Transactions on Microwave Theory and Techniques, vol. 46, no. 9, pp. 1277?1290, Sep 1998. 144