TRANSIENT PORT EXTRACTION FOR OPTIMIZATION OF NON-LINEAR FEEDS TO BROADBAND ANTENNAS By Sean Frederick DePalma A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical and Computer Engineering โ€“ Master of Science 2023 ABSTRACT This work considers the computationally challenging problem of optimizing systems of broadband electromagnetic (EM) devices fed by strongly non-linear circuits. Traditionally, fine-tuning this system would require self-consistent fully coupled analysis of the EM-circuit network in the time domain. Optimization is then dependent on the cost of repeatedly evaluating the tightly-integrated broadband device and strongly non-linear connected circuit. A novel method for transient port parameter extraction enables a circumvention of this computational cost. It constructs a reduced-order EM-circuit representation in which the equivalent EM model is circuit-agnostic. This method of port extraction has been demonstrated to produce equivalent results to that of self-consistent analysis. The circuit- agnostic approach is amenable to any optimization framework given the EM system is invariant. An optimization scheme is developed for application to tightly-coupled broadband EM and non-linear circuit systems that is otherwise untenable with a self- consistent coupled method. This work demonstrates the use of a genetic algorithm to optimize such a system. The optimization scheme is applied to both linear and non-linear circuit systems for objectives including reflection coefficient minimization, reduction of distortion in a non- linear amplifier feeding a broadband antenna, and reverse beam-steering for an array with time-varying incident interference. Copyright by SEAN FREDERICK DEPALMA 2023 To my partner Sierra for her patience and continued support iv ACKNOWLEDGMENTS I would like to acknowledge those who guided me through this work and made this process possible. I extend my gratitude to my advisors and committee Dr. Shanker Balasubramaniam, Dr. Leo Kempel, and Dr. Nelson Sepulveda for their instruction and continued support. To my mentors Omkar Ramachandran and Zane Crawford for their guidance. To my colleagues Chad Moorman, Aditya Purandare, and the Michigan State University Electromagnetics Research Group for their advice and collaboration. To my family for having my back. To the Air Force Institute of Technology and AFROTC Detachment 380 for helping make this opportunity possible. v PREFACE In reference to IEEE copyrighted material which is used with permission in this thesis, the IEEE does not endorse any of Michigan State Universityโ€™s products or ser- vices. Internal or personal use of this material is permitted. If interested in reprint- ing/republishing IEEE copyrighted material for advertising or promotional purposes or for creating new collective works for resale or redistribution, please go to http: //www.ieee.org/publications_standards/publications/rights/rights_link.html to learn how to obtain a License from RightsLink. The views expressed in this article are those of the author and do not reflect the official policy or position of the United States Air Force, Department of Defense, or the U.S. Government. vi TABLE OF CONTENTS CHAPTER 1 INTRODUCTION AND BACKGROUND . . . . . . . . . . . . . . . . 1 CHAPTER 2 PORT PARAMETER EXTRACTION AND OPTIMIZATION OF LINEAR SYSTEMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 CHAPTER 3 PORT PARAMETER EXTRACTION AND OPTIMIZATION OF NON-LINEAR SYSTEMS . . . . . . . . . . . . . . . . . . . . . . . . . . 52 CHAPTER 4 CONCLUSION AND FUTURE WORK . . . . . . . . . . . . . . . . . . 63 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 vii CHAPTER 1 INTRODUCTION AND BACKGROUND 1.1 Introduction Simulation of full-wave electromagnetic (EM) systems attached to non-linear circuits is highly desirable for high-frequency radio frequency (RF) device design. As advances in RF design permit fabrication of components with rapidly decreasing feature size and higher operating frequencies along with tighter integration of components, characterization of coupling effects early in the design process is essential [1, 2]. Combining broadband EM devices with strongly non-linear circuits necessitates ana- lyzing the system in time domain. These tightly-coupled systems are traditionally solved self-consistently. Using a frequency domain approach such as harmonic balance is not feasible for these systems due to the broadband excitation and secondary as well as higher-order harmonics generated due to the inclusion of strongly non-linear elements in the circuit system. RC extraction is another method which may be considered, but it does not sufficiently include the effects of EM and circuit systems that are tightly coupled in their operating bands. A Schur complement approach can be used to extract port parameters for analysis of linear and weakly non-linear circuits and EM systems [3]. The Schur complement approach is circuit-agnostic and computationally efficient as the number of ports is smaller than the number of spatial degrees of freedom. While this is usually the case, this method is incompatible with strongly non-linear systems. A similar method is desirable for tightly-coupled transient analysis. Computing the numerical impulse response in time domain is known to have instabilities associated with deconvolution [4]. A novel approach for transient port parameter extraction in time domain has been developed which obtains results that are equivalent to the solution of a self-consistent method [5, 6]. In transient port parameter extraction, the EM system is circuit-agnostic. Given the 1 EM domain remains unchanged, the port response remains unchanged and can be reused within circuit solving methods while the circuit components and configurations are varied. This system is amenable to optimization by modification of the circuit system feeding a fixed EM device. Self-consistent solution methods for a tightly-coupled broadband EM system fed by non-linear circuits involves an iterative process to converge the non-linear components to a desired tolerance. The iterative process involves multiple inversions of the fully coupled EM-circuit matrix for each time-step. Optimization of these systems becomes prohibitively expensive and is unrealizable when using a fully coupled method. With port extraction, this becomes tenable as the matrix being inverted includes just the non-linear circuit system and port response rather than the EM system (e.g., a much lower number of degrees-of- freedom then the full-wave EM portion). Characterizing the broadband EM device in terms of a port response with port extraction in time domain allows for investigation into optimization of these systems. The key contributions of this thesis are as follows: (1) the development and incorporation of an optimization scheme for tightly-coupled EM and non-linear circuit systems is detailed, (2) the scheme is applied to both linear and non-linear EM-circuit systems, and (3) optimization using port extraction is utilized for a case of adaptive inference cancellation to demonstrate application to mitigation of EM interference. 1.2 Problem Statement A system consists of a radiating object contained within a computational domain ฮฉEM โˆˆ R3 , with ๐‘๐‘ ports each connecting to a lumped circuit system. The circuit system consists of tunable components producing coupling currents JCKT (r, ๐‘ก) across the ports to excite the EM system. Voltage sources within the circuit system are causal and band-limited between some ๐‘“๐‘š๐‘–๐‘› and ๐‘“๐‘š๐‘Ž๐‘ฅ . This analysis enables self-consistent solving for the electric and magnetic fields as well as voltage and currents within each component of the circuit 2 system. Consider the definition of a vector of design variables ๐“” which resides entirely within the tunable circuit system, and some non-linear function ๐‘“ (๐“”) dependent on the EM and circuit systems. It is important to distinguish that the design space of this problem excludes the EM system itself, which shall remain invariant. Using a problem-specific constraint ๐‘‘, we define a cost function ๐ฝ(๐“”) to minimize: ๐ฝ(๐“”) = |๐‘‘ โˆ’ ๐‘“ (๐“”)| (1.1) Our goal is to utilize this cost function to optimize behavior of EM-circuit networks containing non-linear circuit components through exclusively modifying the circuit system. 1.3 Mixed Finite Element Method The mixed finite element method (MFEM) is chosen for this work. MFEM is based upon finite element approximations for a system of partial differential equations, where a continuous domain is represented using a mesh containing a set of discrete points. Using a method such as finite-difference time-domain (FDTD) requires using a step size limited by the size of the mesh cells. For high-frequency simulations, this step size becomes very short due to the minimum physical dimension. The advantage of utilizing MFEM and the flexibility in usable time-step size will later contribute to the feasibility of this work. Working with a large number of steps will become untenable, and using MFEM prevents the time-step size from becoming an impediment. The finite element method (FEM) was originally utilized outside of electrical engineering applications [7]. The first uses of it for electromagnetic devices included investigations into waveguides and cavities in the 1960s [8, 9]. Within FEM, a domain is divided into non-overlapping cells to create a mesh. This mesh is defined by nodes, edges, faces, and cells. To recreate some quantity within the domain, such as electric or magnetic fields, the quantity must be discretized. This discretization is dependent on the number of nodes, edges, faces, or cells, as well as a a basis set. Vector quantities can be approximated using 3 equation (1.2a) for some vector basis set s๐‘– (r), and scalar quantities with equation (1.2b) for some scalar basis set ๐‘”๐‘– (r). ร•๐‘ X(๐‘ก, r) โ‰ˆ ๐‘ฅ ๐‘– (๐‘ก)s๐‘– (r) (1.2a) ๐‘–=1 ร•๐‘ ๐‘Œ(๐‘ก, r) โ‰ˆ ๐‘ฆ ๐‘– (๐‘ก)๐‘”๐‘– (r) (1.2b) ๐‘–=1 This discretization generates a sparse system of equations to work with based on the defined mesh. For application to electromagnetic problems, relevant equations and relations must first be defined. The electric field E(r, ๐‘ก) and magnetic flux density B(r, ๐‘ก) obey Maxwellโ€™s equations. ๐œ•B(r, ๐‘ก) โˆ‡ ร— E(r, ๐‘ก) = โˆ’ (1.3a) ๐œ•๐‘ก ๐œ•D(r, ๐‘ก) โˆ‡ ร— H(r, ๐‘ก) = + J(r, ๐‘ก) (1.3b) ๐œ•๐‘ก โˆ‡ ยท B(r, ๐‘ก) = 0 (1.3c) โˆ‡ ยท D(r, ๐‘ก) = ๐œŒ(r, ๐‘ก) (1.3d) The electric field and magnetic flux density are related through constitutive parameters to the electric flux density and magnetic field. D(r, ๐‘ก) = ๐œ€E(r, ๐‘ก) (1.4a) B(r, ๐‘ก) = ๐œ‡H(r, ๐‘ก) (1.4b) The permittivity and permeability are ๐œ€ = ๐œ€0 ๐œ€๐‘Ÿ and ๐œ‡ = ๐œ‡0 ๐œ‡๐‘Ÿ respectively, where ๐œ‡0 and โˆš ๐œ€0 are defined as 1/ ๐œ‡0 ๐œ€0 = ๐‘ for ๐‘ being the speed of light. The relatively permittivity ๐œ€๐‘Ÿ and relative permeability ๐œ‡๐‘Ÿ are assumed equal to one in this work unless otherwise noted for dielectric substrates. The electric and magnetic fields obey the boundary conditions between two distinct regions. 4 ๐‘›ห† ร— (E1 (r, ๐‘ก) โˆ’ E2 (r, ๐‘ก)) = 0 (1.5a) ๐‘›ห† ร— (H1 (r, ๐‘ก) โˆ’ H2 (r, ๐‘ก)) = J๐‘  (r, ๐‘ก) (1.5b) ๐‘›ห† ยท (D1 (r, ๐‘ก) โˆ’ D2 (r, ๐‘ก)) = ๐œŒ ๐‘  (r, ๐‘ก) (1.5c) ๐‘›ห† ยท (B1 (r, ๐‘ก) โˆ’ B2 (r, ๐‘ก)) = 0 (1.5d) The continuity equation must likewise be satisfied. ๐œ•๐‘ก ๐œŒ(๐‘ก, r) = โˆ’โˆ‡ ยท J(๐‘ก, r) (1.6) Next, a FEM scheme is needed which can self-consistently solve for the fields within the domain ฮฉEM . Initially, applications of FEM to do this used scalar basis functions defined about the nodes of the mesh. Extending this to the vector wave equation results in spurious modes due to not properly satisfying Gaussโ€™ Law and boundary conditions [10]. A variety of approaches were developed to remedy this [11, 12, 13]. One method incorporates vector basis functions [14, 15]. The basis sets were proposed in [16] and satisfy tangential continuity at element boundaries. This satisfies boundary conditions for fields in equation (1.5a) without violating the boundary conditions for fluxes in equation (1.5d). This system enables accurate results with the vector wave equation and Maxwellโ€™s equations, opening investigation into many EM applications [17, 18, 19, 20, 21]. By introducing an additional basis function with normal continuity across cell faces to solve equations (1.3a) and (1.3b), the mixed finite element method is created [22, 23, 24, 25]. This method is used throughout this work. 1.4 Modified Nodal Analysis MNA Formulation The circuit system is solved using the SPICE-like Modified Nodal Analysis (MNA) method [26]. The circuit is represented in a nodal matrix by applying Kirchoffโ€™s current 5 law to non-reference nodes and Kirchoffโ€™s voltage law to independent loops. The sole reference node is typically selected as ground. When applied to the corresponding nodes and loops in the circuit, this generates a set of equations describing the system: YXCKT (๐‘ก) = fCKT (๐‘ก) + fCKT ๐‘›๐‘™ (XCKT , ๐‘ก) (1.7a) CKT CKT ยฉV (๐‘ก)ยช X (๐‘ก) = ยญ ยฎ (1.7b) ICKT (๐‘ก) ยซ ยฌ where Y is the admittance matrix, and VCKT (๐‘ก) and ICKT (๐‘ก) refer to the nodal voltages and branch currents respectively. These are contained within the forcing vector XCKT (๐‘ก), and the excitations due to linear and non-linear components are represented by fCKT (๐‘ก) and fCKT ๐‘›๐‘™ (XCKT , ๐‘ก). If there are no non-linear components, the non-linear contribution to the right-hand side is zero and XCKT (๐‘ก) can be directly solved for at each time-step by inverting the admittance matrix Y. If non-linear components are present, an iterative method must be used to solve for XCKT (๐‘ก) at each time-step. The MNA system can be solved using a method like the Newton- Raphson method, as done in [27]. A non-linear solver tolerance is set before simulation and an initial guess is made for XCKT (๐‘ก). The admittance matrix Y is inverted and the dot product taken with the right-hand side vectors. The voltage across the non-linear components is taken from the result, and the norm is calculated between this and the corresponding voltage in the initial guess. The solution approximation is updated using an appropriate non-linear iterator and the iterative process is repeated until the norm between the approximation and the solve is below the set tolerance. Within both a linear and non-linear MNA solve, Y is inverted. When non-linear components are included and the solve is repeated within an iterative method until the set tolerance is met, this inversion of Y is significant. It may take thousands of iterations, and subsequently thousands of matrix inversions, to converge the non-linear values to the set tolerance. A process reliant on many inversions of a large matrix Y within each time-step 6 can quickly become computationally infeasible. This point is critical for later discussion of EM-circuit simulation. MNA Stamps Applying the KVL and KCL conditions results in each lumped component in the circuit adding to Y in the form of regular matrix โ€™stampsโ€™. This name refers to the standardized formatting approach to adding in, or โ€™stampingโ€™, the contributions and relations for each component into the nodal matrix. These can be derived for the specific time basis functions used to evolve the system in time. Thus, if ๐‘๐‘ is the number of components in the system, the admittance matrix Y can be written as ร•๐‘๐ถ Y= ๐‘€๐‘– (1.8) ๐‘–=1 where ๐‘€ ๐‘– refers to the aforementioned stamp for the ๐‘–th lumped element. Constructing Y requires applying Kirchoffโ€™s voltage law to independent loops. For each component, the first and second rows of the stamp are set by the chosen current reference direction. The last column of the first and second rows correspond to the โ€™positiveโ€™ and โ€™negativeโ€™ nodes based on the reference direction taken. The last row is comprised of admittance information, which changes based on the component. Determining stamps for circuit components with differential voltage-current relations, such as capacitors and inductors, requires discretization of the differential relation. It is necessary to start with the basic differential equations representing these components. ๐‘‘๐‘‰๐ถ (๐‘ก) ๐ผ๐ถ (๐‘ก) = ๐ถ (1.9a) ๐‘‘๐‘ก ๐‘‘๐ผ๐ฟ (๐‘ก) ๐‘‰๐ฟ (๐‘ก) = ๐ฟ (1.9b) ๐‘‘๐‘ก A third-order backwards differentiation formula is used to represent the differential voltage-current relation. In this formula, โ„Ž is the step size, ๐‘› represents the step, ๐‘ก ๐‘› is the 7 current time-step, and ๐‘ฆ๐‘› is the voltage or current value at time-step ๐‘ก ๐‘› [28]. 6 6 18 9 2 โ„Ž๐‘ฆ โ€ฒ(๐‘ก ๐‘›+3 ) = โ„Ž ๐‘“ (๐‘ก ๐‘›+3 , ๐‘ฆ(๐‘ก ๐‘›+3 )) = ๐‘ฆ๐‘›+3 โˆ’ ๐‘ฆ๐‘›+2 + ๐‘ฆ๐‘›+1 โˆ’ ๐‘ฆ๐‘› (1.10) 11 11 11 11 11 The resulting representation is included within the respective element stamps below. The stamps corresponding to the components used in the rest of this work are specified next. Consider ๐‘€๐‘‰ , ๐‘€ ๐ฝ , and ๐‘€๐‘… denote the stamps for a voltage source, current source, and resistor respectively. ๏ฃฏ0 0 1๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๐‘€๐‘‰ = ๏ฃฏ0 (1.11a) ๏ฃฏ ๏ฃบ 0 โˆ’1๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ1 ๏ฃฐ โˆ’1 0 ๏ฃบ๏ฃป ๏ฃฏ0 0 1๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๐‘€ ๐ฝ = ๏ฃฏ0 (1.11b) ๏ฃฏ ๏ฃบ 0 โˆ’1๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ0 0 1 ๏ฃบ๏ฃป ๏ฃฐ ๏ฃฏ0 0 1๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๐‘€๐‘… = ๏ฃฏ 0 (1.11c) ๏ฃฏ ๏ฃบ 0 โˆ’1๏ฃบ ๏ฃฏ ๏ฃบ โˆ’1 ๏ฃฏ1 ๏ฃบ ๏ฃฐ๐‘… ๏ฃฏ ๐‘… โˆ’1๏ฃบ๏ฃป The respective right-hand side contributions for these components follow. The resistor does not contain any independent voltage or current information, so its right-hand side vector is null. ๏ฃฏ 0 ๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๐‘“๐‘‰๐ถ๐พ๐‘‡ (๐‘ก ๐‘› ) = ๏ฃฏ 0 ๏ฃบ (1.12a) ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ๐‘‰(๐‘ก ๐‘› )๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฐ ๏ฃป ๏ฃฏ 0 ๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๐ถ๐พ๐‘‡ ๐‘“๐ผ (๐‘ก ๐‘› ) = ๏ฃฏ 0 ๏ฃบ (1.12b) ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ๐ผ(๐‘ก ๐‘› )๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฐ ๏ฃป 8 ๏ฃฏ0๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๐‘“๐‘…๐ถ๐พ๐‘‡ (๐‘ก ๐‘› ) = ๏ฃฏ0๏ฃบ (1.12c) ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ0๏ฃบ ๏ฃฐ ๏ฃป Applying this formula to the reactive component relations above, the following stamps are produced for the capacitor and inductor. The factor ๐‘‘๐‘ก comes as the time-step size. ๏ฃฎ ๏ฃน ๏ฃฏ0 0 1 ๏ฃบ ๏ฃฏ ๏ฃบ ๐‘€๐ถ = ๏ฃฏ๏ฃฏ0 0 ๏ฃฏ ๏ฃบ โˆ’1 ๏ฃบ๏ฃบ (1.13a) ๏ฃฏ1 โˆ’1 โˆ’6 ๐‘‘๐‘ก ๏ฃบ ๏ฃฏ ๏ฃบ 11 ๐ถ ๏ฃป ๏ฃฏ ๏ฃบ ๏ฃฐ ๏ฃฎ ๏ฃน ๏ฃฏ 0 0 1 ๏ฃบ๏ฃบ ๏ฃฏ ๐‘€๐ฟ = ๏ฃฏ๏ฃฏ 0 ๏ฃฏ ๏ฃบ 0 โˆ’1๏ฃบ๏ฃบ (1.13b) ๏ฃฏ โˆ’11 ๐‘‘๐‘ก 11 ๐‘‘๐‘ก ๏ฃบ ๏ฃฏ ๏ฃฏ 6 ๐ฟ 1 ๏ฃบ๏ฃบ ๏ฃฐ 6 ๐ฟ ๏ฃป The corresponding right-hand side contributions are as follows. 0 ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๐ถ๐พ๐‘‡ ๐‘“๐ถ (๐‘ก ๐‘› ) = ๏ฃฏ (1.14a) ๏ฃฏ ๏ฃบ 0 ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๐‘‰๐ถ (๐‘ก ๐‘›โˆ’1 ) โˆ’ 9 ๐‘‰๐ถ (๐‘ก ๐‘›โˆ’2 ) + 2 ๐‘‰๐ถ (๐‘ก ๐‘›โˆ’3 )๏ฃบ ๏ฃฏ 18 ๏ฃบ ๏ฃฐ 11 11 11 ๏ฃป 0 ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๐‘“๐ฟ๐ถ๐พ๐‘‡ (๐‘ก ๐‘› ) = ๏ฃฏ (1.14b) ๏ฃฏ ๏ฃบ 0 ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๐ผ๐ฟ (๐‘ก ๐‘›โˆ’1 ) โˆ’ 9 ๐ผ๐ฟ (๐‘ก ๐‘›โˆ’2 ) + 2 ๐ผ๐ฟ (๐‘ก ๐‘›โˆ’3 )๏ฃบ ๏ฃฏ 18 ๏ฃบ ๏ฃฐ 11 11 11 ๏ฃป Next, we develop an ideal operational amplifier stamp for some of the system optimization problems studied in this work. The ideal operational amplifier follows the constraints of zero input potential difference, and zero input current into both positive and negative terminals [29]. Forming these constraints into the MNA system creates the following ideal operational amplifier stamp. 9 ๏ฃฏ0 0 0 0๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ0 0 0 0๏ฃบ ๏ฃฏ ๏ฃบ ๐‘€๐‘‚ = ๏ฃฏ๏ฃฏ ๏ฃบ ๏ฃบ (1.15) ๏ฃฏ0 0 0 1๏ฃบ๏ฃบ ๏ฃฏ ๏ฃฏ ๏ฃบ ๏ฃฏโˆ’1 1 0 0๏ฃบ๏ฃป ๏ฃฐ The component does not contain any independent voltage or current relations beyond the constraints, leading to a null right-hand side vector. ๏ฃฏ0๏ฃบ ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ0๏ฃบ ๏ฃฏ ๏ฃบ ๐ถ๐พ๐‘‡ ๐‘“๐‘‚ (๐‘ก ๐‘› ) = ๏ฃฏ๏ฃฏ ๏ฃบ๏ฃบ (1.16) ๏ฃฏ0๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ0๏ฃบ ๏ฃฐ ๏ฃป Lastly, a diode is established for use in non-linear circuit systems. The diode behavior is modeled using the Shockley equation and its first derivative with respect to the diode voltage. In the equation, ๐ผ๐‘† represents the diode saturation current, ๐‘‰๐ท represents the diode voltage, ๐œ‚ represents the ideality factor, and ๐‘‰๐‘‡ is the thermal voltage. In this notation, the current is set equal to a placeholder ๐›ผ ๐ท and the derivative is likewise set to ๐›ฝ ๐ท . ๐‘‰๐ท (๐‘ก ๐‘› ) ( ) ๐›ผ ๐ท = ๐ผ(๐‘ก ๐‘› ) = ๐ผ๐‘† ๐‘’ ๐œ‚๐‘‰๐‘‡ (1.17a) ๐‘‰๐ท (๐‘ก ๐‘› ) ๐‘‘๐ผ(๐‘ก ๐‘› ) ๐ผ๐‘† ( ๐œ‚๐‘‰๐‘‡ ) ๐›ฝ๐ท = = ๐‘’ (1.17b) ๐‘‘๐‘‰๐ท ๐œ‚๐‘‰๐‘‡ Using these equations, the diode stamp is then created. Fundamentally, the diode behaves like a current-dependant resistor with resistance 1/๐›ฝ ๐ท attached in parallel to a current source of amplitude ๐›ผ ๐ท . This yields a stamp as follows ๏ฃฏ ๐›ฝ ๐ท โˆ’๐›ฝ ๐ท ๏ฃบ ๏ฃฎ ๏ฃน ๐‘€๐ท = ๏ฃฏ๏ฃฏ ๏ฃบ (1.18) ๏ฃฏโˆ’๐›ฝ ๐ท ๐›ฝ ๐ท ๏ฃบ ๏ฃบ ๏ฃฐ ๏ฃป 10 The subsequent right-hand side contribution for the diode can be written as ๏ฃฏโˆ’๐›ผ ๐ท + ๐›ฝ ๐ท ๐‘‰๐ท (๐‘ก ๐‘› )๏ฃบ ๏ฃฎ ๏ฃน ๐‘“๐ท๐ถ๐พ๐‘‡ (๐‘ก ๐‘› ) = ๏ฃฏ๏ฃฏ ๏ฃบ. (1.19) ๏ฃฏ ๐›ผ ๐ท โˆ’ ๐›ฝ ๐ท ๐‘‰๐ท (๐‘ก ๐‘› ) ๏ฃบ ๏ฃบ ๏ฃฐ ๏ฃป 1.5 Transient Port Parameter Extraction A method for finding self-consistent solutions to a MFEM EM-circuit system in time domain allows for representing the EM-circuit system with a reduced-order model [5, 6]. Port extraction in time domain finds port responses which are equivalent to the numerical impulse response at each port within the EM system. The response is incorporated into the MNA solution method in place of the EM system. Given ฮฉ๐ธ๐‘€ is unchanged, the port response remains unchanged and can be reused within the MNA process while circuit components and configuration are modified. The solution produced by port extraction is equivalent to that of a self-consistent, fully coupled solution. A Maxwell solver is formulated with port extraction in time domain. The solver is used for the EM domain ฮฉ๐ธ๐‘€ . This domain is bounded by a truncating surface ๐œ•ฮฉ๐ธ๐‘€ . Boundary conditions must be satisfied within the domain. This is done by further defining surfaces ฮ“๐ท and ฮ“๐‘ as a subset of ๐œ•ฮฉ๐ธ๐‘€ to represent PEC and PMC regions of the boundary respectively. ๐‘›ห† ร— E(r, ๐‘ก) = 0 for r โˆˆ ฮ“๐ท (1.20a) ๐‘›ห† ร— H(r, ๐‘ก) = 0 for r โˆˆ ฮฉ๐‘ (1.20b) Within this work, ฮฉ๐ธ๐‘€ is truncated using a perfectly matched layers (PML) region [30, 31]. The exterior of the PML is defined as a PEC surface. The implementation of PML regions follows a stretched coordinate system in which the convolutions resulting from the system are directly evaluated. The stretched coordinate system is defined using the 11 transformation ๐‘  ๐‘  ยฉ ๐‘ฆ๐‘  ๐‘ฅ ๐‘ง 0 0 ยช ยญ ยฎ ฮ›(๐œ”) = ยญ 0 ยญ ๐‘ ๐‘ฅ ๐‘ ๐‘ง ยฎ (1.21) ๐‘ ๐‘ฆ 0 ยฎ ยญ ยฎ ยญ ๐‘ ๐‘ฅ ๐‘ ๐‘ฆ ยฎ 0 0 ๐‘ ๐‘ง ยฌ ยซ ๐œŽ๐‘– where ๐‘  ๐‘– = 1 + is used to match the absorbing layers to free space, and ๐œŽ๐‘– govern field ๐‘—๐œ”๐œ–0 loss. This system alters Maxwellโ€™s equations in frequency domain: ฮ›(๐œ”)โˆ’1 ยท โˆ‡ ร— E(r, ๐œ”) = โˆ’ฮ›(๐œ”)โˆ’1 ยท ๐‘—๐œ”B(r, ๐œ”) (1.22) โˆ’1 โˆ’1 โˆ‡ ร— ๐œ‡ ฮ›(๐œ”) ยท B(r, ๐œ”) = J(r, ๐œ”) + ๐œ–(r)ฮ›(๐œ”) ยท ๐‘—๐œ”E(r, ๐œ”) These equations are inverse Fourier transformed to obtain a time-marching scheme, resulting in ๐œ•B(r, ๐‘ก) L2 (๐‘ก) โˆ— โˆ‡ ร— E(r, ๐‘ก) = โˆ’L2 (๐‘ก) โˆ— (1.23a) ๐œ•๐‘ก โˆ‡ ร— ๐œ‡โˆ’1 L2 (๐‘ก) โˆ— B(r, ๐‘ก) = J(r, ๐‘ก) + ๐œ– 0 L1 (๐‘ก) โˆ— E(r, ๐‘ก) (1.23b) where โˆ— represents temporal convolution and L1 (๐‘ก) = F โˆ’1 (๐‘—๐œ”ฮ›(๐œ”)) (1.24a)   L2 = F โˆ’1 ฮ›(๐œ”)โˆ’1 (1.24b) To discretize this system, equation (1.23a) is tested with a W2 (r) basis function and (1.23b) is tested with W1 (r). PEC and PMC surfaces obey the respective boundary conditions (1.5a) and (1.5b). The electric and magnetic fields likewise vary as (1.3a) and (1.3b) respectively. Variational Formulation Variational equations with appropriate function spaces are created to find solutions to the equations in (1.3). The spaces ๐ป(๐‘๐‘ข๐‘Ÿ๐‘™; ฮฉ) and ๐ป(๐‘‘๐‘–๐‘ฃ; ฮฉ) are suitable for E(r, ๐‘ก) and B(r, ๐‘ก) respectively. ๐ป(๐‘๐‘ข๐‘Ÿ๐‘™; ฮฉ) ensures the tangential components of the function are 12 continuous across cell faces, and ๐ป(๐‘‘๐‘–๐‘ฃ; ฮฉ) ensures normal continuity of the function across cell faces. ๐ป(๐‘๐‘ข๐‘Ÿ๐‘™; ฮฉ) = {u โˆˆ L2 (ฮฉ); โˆ‡ ร— u โˆˆ L2 (ฮฉ)} (1.25) ๐ป(๐‘‘๐‘–๐‘ฃ; ฮฉ) = {u โˆˆ L (ฮฉ); โˆ‡ ยท u โˆˆ L (ฮฉ)} 2 2 The variational equations are constructed by taking the inner product of (1.3a) with the Whitney edge basis function Bโˆ— โˆˆ ๐ป(๐‘‘๐‘–๐‘ฃ; ฮฉ) and (1.3b) with the Whitney face basis function Eโˆ— โˆˆ ๐ป(๐‘๐‘ข๐‘Ÿ๐‘™; ฮฉ) [22]. ๐œ•B โˆ— โˆซ โˆซ 1 ๐‘‘๐‘‰๐œ‡ (r)โˆ’1 ยทB =โˆ’ ๐‘‘๐‘‰ โˆ‡ ร— E ยท Bโˆ— ๐œ•๐‘ก ๐œ‡(r) ฮฉ ฮฉ (1.26) ๐œ•E โˆ— โˆซ โˆซ   โˆ’1 โˆ— ๐‘‘๐‘‰ ๐œ–(r) ยทE =โˆ’ ๐‘‘๐‘‰ โˆ‡ ร— ๐œ‡ (r)B โˆ’ J CKT ยทE ๐œ•๐‘ก ฮฉ ฮฉ The fields are incorporated into the MFEM scheme using Whitney edge and face basis functions defined on the edges and faces of the mesh to represent E(r, ๐‘ก) and B(r, ๐‘ก). With ๐‘๐‘’ representing number of edges and ๐‘ ๐‘“ the number of faces in the mesh, this takes form as ร•๐‘๐‘’ E(r, ๐‘ก) = ๐‘’ ๐‘– (๐‘ก)W1๐‘– (r) (1.27a) ๐‘–=1 ๐‘๐‘“ ร• B(r, ๐‘ก) = ๐‘ ๐‘– (๐‘ก)W2๐‘– (r) (1.27b) ๐‘–=1 The vector e(๐‘ก) can be defined as a vector of ๐‘’ ๐‘– (๐‘ก) from ๐‘’1 (๐‘ก) to ๐‘’ ๐‘๐‘’ (๐‘ก), and likewise b(๐‘ก) as a vector of ๐‘ ๐‘– (๐‘ก) from ๐‘ 1 (๐‘ก) to ๐‘ ๐‘ ๐‘“ (๐‘ก). Galerkinโ€™s method is used to formulate the semi-discrete Maxwell system from (1.26): ๐œ•๐‘ก M๐œ‡โˆ’1 b(๐‘ก) = โˆ’ M๐‘ e(๐‘ก), (1.28) ๐œ•๐‘ก M๐œ– e(๐‘ก) =M๐‘‡๐‘ b(๐‘ก) โˆ’ M๐ผ e(๐‘ก) + JCKT (๐‘ก) 13 where the matrices are defined as โˆซ [M๐œ– ]๐‘–,๐‘— = ๐‘‘๐‘‰ ๐œ–(r)W1๐‘– (r) ยท W1๐‘— (r) โˆซฮฉ 1 ๐‘‘๐‘‰ W2๐‘– (r) ยท W2๐‘— (r)   M๐œ‡โˆ’1 = ๐‘–,๐‘— ๐œ‡(r) โˆซฮฉ 1 [M๐‘ ]๐‘–,๐‘— = ๐‘‘๐‘‰ W2๐‘– (r) ยท โˆ‡ ร— W1๐‘— (r) (1.29) ฮฉ ๐œ‡(r) โˆซ [M๐ผ ]๐‘–,๐‘— = ๐‘‘๐‘‰๐œ‚(r)๐‘›ห† ร— W1๐‘– (r) ยท ๐‘›ห† ร— W1๐‘— (r) โˆซฮ“๐ผ jCKT (๐‘ก) ๐‘‘๐‘‰W1๐‘– (r) ยท JCKT (r, ๐‘ก)   ๐‘– = ฮฉ The Newmark-๐›ฝ time-marching scheme for solving second order differential equations is used with this system [32, 33, 34]. Fixing ๐›พ = 0.5 and ๐›ฝ = 0.25 within the process forms an unconditionally stable time-marching scheme [35], solving for e(๐‘ก) and b(๐‘ก) at uniform steps in time. A testing function ๐‘Š(๐‘ก) is defined as ๐‘ก ๐‘› โˆ’๐‘ก ๏ฃฑ ๐‘ก โˆˆ [๐‘ก ๐‘›โˆ’1 , ๐‘ก ๐‘› ] ๏ฃด ๏ฃด ๏ฃด ๏ฃด ๏ฃด ฮ”๐‘ก ๏ฃด ๏ฃฒ ๏ฃด ๐‘Š(๐‘ก) = ๐‘กโˆ’๐‘ก ๐‘› ๐‘ก โˆˆ [๐‘ก ๐‘› , ๐‘ก ๐‘›+1 ] (1.30) ๏ฃด ฮ”๐‘ก ๏ฃด ๏ฃด ๏ฃด ๏ฃด ๏ฃด0 ๏ฃด otherwise ๏ฃณ Representing the vectors e(๐‘ก) and b(๐‘ก) in terms of second order Lagrange polynomials and testing (1.28) with ๐‘Š(๐‘ก) results in the recurrence relation [0.5A1 + 0.25ฮ”๐‘ก A0 ] x๐‘›+1 + [0.5ฮ”๐‘ก A0 ] x๐‘› (1.31) + [โˆ’0.5A1 + 0.25ฮ”๐‘ก A0 ] x๐‘›โˆ’1 = 0.25ฮ”๐‘ก f๐‘›+1 + 0.5ฮ”๐‘ก f๐‘›+1 + 0.25ฮ”๐‘ก f๐‘›+1 14 With matrix entries ยฉM๐œ‡โˆ’1 0 ยช A1 = ยญ ยฎ ยซ 0 M๐œ– ยฌ ยฉ 0 M๐‘ ยช A0 = ยญ ยฎ (1.32) ๐‘‡ ยซโˆ’M๐‘ M๐ผ ยฌ 0 f๐‘› = ยญ ยฉ ยช ยฎ jCKT (๐‘ก ) ยซ ๐‘› ยฌ Incorporating the Circuit System with the Time-Marching Scheme Returning to the circuit system described by (1.7a) and (1.7b), the system is discretized by using fourth-order backward-looking Lagrange polynomials to represent nodal voltages V(๐‘ก) in (1.7b) as ๐‘– ร– ๐‘ก โˆ’ ๐‘ก๐‘— ๐ฟ ๐‘– (๐‘ก) = ๐‘ก๐‘– โˆ’ ๐‘ก ๐‘— ๐‘—=๐‘–โˆ’3 ๐‘—โ‰ ๐‘– (1.33) ร•๐‘๐‘ก V(๐‘ก) = xCKT ๐‘– ๐ฟ ๐‘– (๐‘ก) ๐‘–=4 Testing this representation with a delta function centered at the ๐‘ก ๐‘›+1 time-step yields a time-marching scheme for the circuit system.   11 3 Y1 + ฮ”๐‘ก Y0 xCKT ๐‘›+1 = 3Y1 xCKT ๐‘› โˆ’ Y1 xCKT ๐‘›โˆ’1 6 2 1 (1.34) + Y1 xCKT ๐‘›โˆ’2 + ฮ”๐‘ก f CKT (๐‘ก ๐‘›+1 ) 3   + ฮ”๐‘ก fCKT๐‘›๐‘™ ๐‘ก ๐‘›+1 , xCKT ๐‘›+1 Coupling the EM and Circuit Systems The EM and circuit systems are coupled by the fields and current densities in the EM system and the voltages and currents in each circuit system. Each EM-circuit connection 15 spans a coupling edge l ๐‘˜ . For the ๐‘˜ th coupling edge of the EM system connecting to the ๐‘— th circuit system, the voltage across the edge ๐‘‰๐‘—CKT (๐‘ก) at a time ๐‘ก ๐‘– can be represented with a coupling coefficient ๐ถ ๐‘— ๐‘˜ which relates the EM system fields to the circuit system voltages. โˆซ โŸจ๐›ฟ(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐‘‰๐‘—CKT (๐‘ก)โŸฉ = โŸจ๐›ฟ(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐‘’ ๐‘˜ (๐‘ก) lฬ‚ ๐‘˜ ยท W1๐‘˜ (r)๐‘‘rโŸฉ |l ๐‘˜ | (1.35) = โŸจ๐›ฟ(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐‘’ ๐‘˜ (๐‘ก)๐ถ ๐‘— ๐‘˜ โŸฉ ๐ผ CP ๐‘— (๐‘ก) likewise represents the magnitude of the current impressed upon the EM system over the coupling edge by a circuit system. For the ๐‘˜ th edge connecting to the ๐‘— th circuit system, the current density JCKT can similarly be represented with a coupling coefficient ๐ถ ๐‘˜ ๐‘— . This relates circuit system currents to the EM system current densities through testing with the function ๐‘Š(๐‘ก) from (1.30). โˆซ โŸจ๐‘Š(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐ฝ๐‘˜CKT (๐‘ก)โŸฉ = โŸจ๐‘Š(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐ผ CP ๐‘— (๐‘ก) lฬ‚ ๐‘˜ ยท W1๐‘˜ ๐‘‘rโŸฉ |l ๐‘˜ | (1.36) = โŸจ๐‘Š(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐ผ CP ๐‘— (๐‘ก)๐ถ ๐‘˜ ๐‘— โŸฉ The edges in the FEM system are assumed to be infinitesimally thin, meaning the current density seen by the EM system can be approximated with the current on the circuit port edge. The feed model used within the system is a current probe. Additionally, the choice of testing functions in (1.35) and (1.36) creates identical coupling coefficients. Applying the coupling coefficients to (1.34) and (1.28) results in the following matrix system. ๏ฃฎ EM ๏ฃฏM 0 0.25C๏ฃบ ๏ฃฏ xEM ๏ฃบ ๏ฃฏ bEM ๏ฃบ ๏ฃน ๏ฃฎ ๏ฃน ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ CKT ๏ฃบ ๏ฃฏ CKT ๏ฃบ ๏ฃฏ CKT ๏ฃบ ยท = (1.37) ๏ฃฏ ๏ฃฏ 0 M B ๏ฃบ ๏ฃฏx ๏ฃบ ๏ฃฏb ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ T T ๏ฃบ ๏ฃฏ CP ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ C ๏ฃฐ B 0 ๏ฃบ ๏ฃป ๏ฃฐ ๏ฃฏ I ๏ฃบ ๏ฃป ๏ฃฐ ๏ฃฏ 0 ๏ฃบ ๏ฃป Within this system, MEM and bEM refer to the left- and right-hand sides of (1.28) respectively, MCKT and bCKT refer to the left- and right-hand side of (1.34), ICP is a vector containing the coupling currents, and B is a logical matrix connecting couplings edge to the respective voltage nodes in the circuit system. 16 Temporal Port Extraction The objective of temporal port extraction is to find the response of the EM system due to an excitation applied to each coupling edge. The response can then be convolved with the current across the edge to solve for fields within ฮฉ๐ธ๐‘€ . Consider a circuit port ๐‘ž spanning the set of FEM edges ๐‘(๐‘ž). A spatial excitation ๐‘’ ๐‘ž,๐‘˜ is defined and used to further define a temporal excitation e๐‘ž (๐‘ก). ๏ฃฑ ๏ฃฒ 1 if ๐‘˜ โˆˆ ๐‘(๐‘ž) ๏ฃด ๏ฃด ๏ฃด ๐‘’ ๐‘ž,๐‘˜ = (1.38a) ๏ฃด ๏ฃด 0 otherwise ๏ฃด ๏ฃณ ๏ฃฑ ๏ฃฒ e๐‘ž if ๐‘ก = ๐›ฟ ๐‘ก ๏ฃด ๏ฃด ๏ฃด e๐‘ž (๐‘ก) = (1.38b) ๏ฃด ๏ฃด 0 otherwise ๏ฃด ๏ฃณ Using the temporal excitation as the forcing function f within (1.28), the recurrence relation (1.31) becomes ๐‘ž [0.5A1 + 0.25ฮ”๐‘ก A0 ] x๐‘›+1 ๐‘ž + [0.5ฮ”๐‘ก A0 ] x๐‘› ๐‘ž + [โˆ’0.5A1 + 0.25ฮ”๐‘ก A0 ] x๐‘›โˆ’1 (1.39) = 0.25ฮ”๐‘ก e๐‘ž (๐‘ก ๐‘›+1 ) + 0.5ฮ”๐‘ก e๐‘ž (๐‘ก ๐‘› ) + 0.25ฮ”๐‘ก e๐‘ž (๐‘ก ๐‘›โˆ’1 ) From the vector ๐‘ฅ ๐‘ž solved from (1.39), the response matrix G can be constructed. ๐บ ๐‘˜๐‘ž = ๐‘ฅ ๐‘˜,๐‘(๐‘ž) (1.40) This allows for representation of ๐‘‰๐‘—CKT (๐‘ก) in (1.35) in terms of the coupling coefficient ๐ถ ๐‘— ๐‘˜ and the convolution of the response G with the coupling current ๐ผ ๐ถ๐‘ƒ . ๐‘๐‘ ร• ๐‘‰๐‘—CKT (๐‘ก ๐‘– ) = โŸจ๐›ฟ(๐‘ก โˆ’ ๐‘ก ๐‘– ), ๐ถ ๐‘— ๐‘˜ ๐บ ๐‘˜๐‘ž (๐‘ก) โˆ— ๐ผ ๐‘žCP (๐‘ก)โŸฉ (1.41) ๐‘ž=1 17 By extracting the response at each port in this way, the EM-circuit network can be represented using a lower-order model through characterizing the EM system in terms of the responses. The constructed EM model is circuit-agnostic, allowing for modifications to circuit configuration after port extraction. This model is therefore amenable to circuit modification given the EM system and its subsequent port response are unmodified. 1.6 Optimization Using Port Extraction General Optimization Framework A general non-constrained optimization problem can be defined using input variables ๐‘ฅ ๐‘› , objectives ๐‘“๐‘š , and a convergence criterion such as minimization of the objectives [36]. The problem has ๐‘ input variables subject to lower bound ๐‘ฅ ๐‘›๐ฟ and upper bound ๐‘ฅ ๐‘ˆ ๐‘› , used to minimize ๐‘€ objectives ๐‘“๐‘š using a defined operation within the body of the optimization problem. The input variables are defined within a space ๐‘ฅ โˆˆ ฮฉ. minimizing ๐‘“๐‘š , ๐‘š = 1, ..., ๐‘€ ๐‘ฅ ๐‘›๐ฟ โ‰ค ๐‘ฅ ๐‘› โ‰ค ๐‘ฅ ๐‘ˆ ๐‘› , ๐‘› = 1, ..., ๐‘ ๐‘ฅ โˆˆ ฮฉ When maximization is desired in place of minimization, the negative of the minimum can be used as the objective. For ๐‘€ = 1, the problem is single-objective and produces a single minimum result. For ๐‘€ > 1, the problem becomes multi-objective and produces a set of viable solutions representing multiple minimums, requiring further processing to isolate a desired solution. In the optimization scheme developed for this thesis, the cost function ๐ฝ(๐“”) from (1.1) is used as the objective ๐‘“๐‘š , with the ๐‘-size vector of input variables ๐“” taking value as ๐‘ฅ ๐‘› and the non-linear function ๐‘“ (๐“”) found from the results of simulating the EM-circuit system. Then, ๐‘‘ takes value as a desired objective value. Solutions to this optimization scheme minimize distance between desired value ๐‘‘ and calculated system 18 response ๐‘“ (๐“”) to circuit inputs ๐“”. Throughout optimization, the EM system is unchanged. Combining Port Extraction with a Genetic Algorithm Optimization of a tightly-coupled broadband EM system fed by strongly non-linear circuits is prohibitively expensive for large EM meshes. Each step in time requires an iterative method to converge non-linear quantities within the circuit system to a desired tolerance, and each evaluation within this iterative method involves inverting the fully coupled EM-circuit matrix. Consider some EM device and surrounding radiation box discretized with tetrahedral elements having average edge length ๐œ†min /30, or thirty times smaller than the minimum wavelength analyzed in the system; the size of the mesh describing this discretized EM system can reach up to a million or more elements in size. Inverting the resulting EM matrix many times within an iterative method, and further repeating the iterative method for each time-step makes optimization by repeatedly changing and resolving the coupled EM-circuit system unrealizable. Temporal port extraction characterizes the EM system in terms of a response for each port in the system, allowing for representation of the EM system in terms of these responses rather than the discretized EM system matrix. This lower-order model finds solutions equal to a self-consistent method. Furthermore, this allows for optimization of the circuit system given the EM system and its port responses are completely unchanged. Optimization of broadband antennas with strongly non-linear circuit feeds then becomes possible using port extraction to optimize the total network by exclusively modifying circuit components and configuration. A comparison of this concept using self-consistent fully coupled and port extraction methods is compared in Figure 1.1. In this figure, the self-consistent method shown on the left illustrates the fully coupled matrix inversion repeated within the iterative method for non-linear components, which is then within the optimization process. Comparatively, port extraction involves EM matrix inversion only to find port responses, which is independent from the iterative method and circuit-modifying optimization 19 process. Figure 1.1: Comparison of non-linear circuit optimization tied to EM system using self- consistent (denoted coupled) and port extraction methods This scheme is amenable to any optimization method. For tightly-coupled systems of broadband EM devices fed with strongly non-linear circuits, this work will investigate an optimization scheme which can sample the input space between specified bounds, solve the system using port extraction, and rank produced solutions compared to a desired value. Successive iterations which repeat this process and bias future input selections towards those which produced solutions closest to the desired value, converging towards the desired outcome. Furthermore, choosing a non-gradient optimization method offers more flexibility for strongly non-linear system responses, as gradient methods may converge to 20 local minima within the response instead of the global minimum. The chosen method is the genetic algorithm. Pymoo is a multi-objective python optimization module which has a simple genetic algorithm [36]. This algorithm will be used throughout this thesis. The genetic algorithm selected begins with a sampled population from within the provided input bounds. The body of the algorithm is then solved. In this application, the circuit solver of port extraction and subsequent processing is used to calculate a desired system characteristic, such as reflection coefficient. The solutions are then ranked by fitness. For this simple genetic algorithm, this is taken as the solution yielding minimum cost function. This requires setting an objective of (1.1) for which solutions are represented using distance from the desired objective. A crossover method is then used to select fit solutions for mating. This produces offspring within the next generation which carry on successful traits from the previous generation. Next, a mutation method is applied to increase diversity within the population. Finally, the algorithm iterates and repeats this process until a termination condition is met. The settings modified within the genetic algorithm include specifying a population size, number of offspring per generation, total number of generations, method of sampling, crossover method with probability and coefficient ๐œ‚๐ถ , mutation method with probability and coefficient ๐œ‚ ๐‘€ , and termination method. Most of the settings remain unchanged throughout all applications in this thesis. The termination method selected is to terminate after simulating a total number of generations. The mutation method is set as polynomial mutation and the mutation probability is fixed at 0.8 [37]. The crossover method is set as simulated binary crossover [38]. The crossover coefficient ๐œ‚๐ถ is fixed at 15, and the crossover probability is fixed to the mutation probability. The method of sampling the inputs between provided bounds is set as random float sampling. This method generates a random floating point number between 0 and 1, and scales it between the inputs of the most-fit and least-fit solutions within each generation. 21 Mutation and crossover increase diversity in the post-sampled population before the next generation is simulated. Mutation coefficient, population, offspring, and generations vary case to case. 1.7 Goals and Outline This thesis presents a method for using temporal port extraction for optimization of tightly-coupled broadband EM devices fed by strongly non-linear circuit systems. The principle contributions of this work include developing a scheme which utilizes port extraction to create a circuit-agnostic EM system which is amenable to optimization by exclusively changing circuit components. This is demonstrated in using a genetic algorithm for single-objective and multiple-objective optimization of EM systems fed by non-linear circuits, and applying this scheme to a case of adaptive interference cancellation for signal to interference and noise ratio (SINR) maximization. The outline of the thesis is as follows. Chapter 2 presents methodology and examples for constructing these optimization cases with linear circuits. Validation of port extraction is made in comparison to a measured device and to self-consistent simulation. Optimization is performed for circuit-only configurations, single-objective EM-circuit systems, multi- objective EM-circuit systems, and a linear circuit approach to the SINR case. Chapter 3 presents non-linear adaptions to the work of chapter 2, again validating port extraction against self-consistent simulation as well as a non-linear circuit approach to the SINR case and reduction of distortion in a non-linear amplifying circuit feeding a broadband antenna. Chapter 4 presents conclusions and suggestions for future work. 22 CHAPTER 2 PORT PARAMETER EXTRACTION AND OPTIMIZATION OF LINEAR SYSTEMS 2.1 Introduction Port extraction in time domain produces a solution equivalent to that of self-consistent analysis. Given the EM domain is invariant, the port responses remain unchanged and can be reused with varying circuit configurations. The EM-circuit network can then be optimized exclusively by circuit modification as shown in the right-hand side of Figure 1.1. By representing the EM system in terms of port responses and reusing them with changing circuit configurations, optimization of these systems becomes realizable whereas it is untenable using self-consistent analysis. To validate the optimization scheme developed in 1.6, linear cases are first investigated and the results are compared to the analytical or known solution where applicable. In this chapter, the results of port extraction for analysis of linear systems are compared to that of a fabricated and measured monopole as well as self-consistent simulation of a broadband Vivaldi antenna. An optimization scheme is developed starting with linear circuit configurations and working towards a more complex case of adaptive interference cancellation using an antenna array fed by linear circuits. These cases build a scheme applicable to the optimization of non-linear EM-circuit systems in chapter 3. 2.2 Port Extraction Comparisons Methodology for optimization of linear EM-circuit networks using transient port extraction is developed to build towards non-linear systems later on. This establishes proof-of-concept for subsequent non-linear experiments while working with linear circuits where the response can typically be calculated and compared to the optimized result. A linear system in this sense refers to an EM system fed by one or more circuits 23 containing lumped components which have a linear relation between voltage and current, such as resistors, capacitors, and inductors. To begin testing linear systems, an EM object must first be constructed and discretized in a geometry and mesh generation software. After extracting port responses from the EM system and solving the subsequent circuit system as detailed in 1.5 and 1.4 respectively, the resulting voltages and currents can be used to determine system characteristics such as radiated power or reflection coefficients. The voltage sources used are causal and defined as modified Gaussian signals with 2 /2๐œŽ 2 ๐‘ฃ(๐‘ก) = cos(2๐œ‹ ๐‘“0 ๐‘ก)๐‘’ โˆ’(๐‘กโˆ’6๐œŽ) where ๐œŽ = 3/(2๐œ‹ ๐‘“bw )โˆ’1 , using a maximum frequency found as ๐‘“max = ๐‘“0 + ๐‘“bw . The time-step size is defined as ฮ”๐‘ก = (20 ๐‘“๐‘š๐‘Ž๐‘ฅ )โˆ’1 . These definitions for a voltage source and time-step size hold constant for the remainder of this work unless otherwise specified. Comparison to Measured Device A comparison is made between port extraction and a fabricated and measured device. Radiation impedance and radiated power are measured for a strip above a finite ground plane [1]. This is a notional monopole. The strip measures 16.51cm in length and 2.51cm in width, and sits 1cm above a rectangular ground plane measuring 29.21cm by 30.48cm. This device is depicted in Figure 2.1. 24 Figure 2.1: Geometry of strip monopole above finite ground plane (dimensions in cm) The radiation impedance and power radiated due to a 1mW source voltage are calculated for the measured device using: 1 ๐‘… ๐‘Ÿ ๐‘Ž๐‘‘ ๐‘ƒ๐‘Ÿ ๐‘Ž๐‘‘ = ๐‘‰2 (2.1) 2 ๐‘… + ๐‘‹2 2 ๐‘Ÿ ๐‘Ž๐‘‘ ๐‘Ÿ ๐‘Ž๐‘‘ A model of the strip monopole above a finite ground plane is constructed based on the dimensions in Figure 2.1. A feeding circuit is attached to an edge centered between the strip and ground plane. The circuit consists of a unit-amplitude voltage source defined as a modulated Gaussian wave with center frequency 1GHz and bandwidth 900MHz, in series with a 100ฮฉ resistor. The geometry is centered within a radiation box extending 2 ๐œ†max from all sides of the system. A PML region is defined, extending 13 ๐œ†max beyond the 1 outer faces of the radiation box [31]. The entire geometry is discretized with tetrahedral mesh elements (tets) having average edge length ๐œ†min /30. The resulting model consists of 1.98M mesh elements. 25 The response of the EM system at the single port is found and used with the feeding circuit to find simulated radiation impedance. A comparison of this to the measured device is shown in Figure 2.2. Equation (2.1) is used to ensure identical calculation of radiated power to the publication for which the measured result is referenced from. Using the impedance found with (2.1) for a 1mW voltage source, the power radiated by the system is calculated and compared to the measured data in Figure 2.3. The source of the measured data compared the power radiated to their own finite-difference time-domain (FD-TD) simulation [1]. Their FD-TD comparison is additionally shown against transient port extraction. While the results plotted show relative agreement, the discrepancy between simulated and measured data can in part be explained by quality of the constructed model. A similar comparison is made for port extraction and the measured data of this device in [6], which better matches the measured radiated power. This discrepancy can be remedied by employing finer discretization around the port in the antenna model, to better capture higher-order field behavior near the excitation. 26 Figure 2.2: Measured radiation impedance of strip monopole compared to simulated results using transient port extraction 27 Figure 2.3: Comparison of radiated power for measured and FD-TD data versus port extraction results Analysis of a Vivaldi Antenna A comparison of the accuracy between port extraction and self-consistent analysis is made using an exponentially-tapered Vivaldi antenna. The antenna consists of a RT/duroid 5880 substrate measuring 41mm in length, 37.5mm in width, and 1.575mm in height, with symmetric PEC regions on either face of the substrate. The PEC regions taper exponentially from one end of the substrate face to the opposite end, with a 3.5mm square balun at the base of the taper. The antenna is fed using a port on the outer edge of the substrate near the balun, spanning from the middle of the substrate height to one of the symmetric PEC faces. The feeding circuit consists of a voltage source with center frequency 10GHz and bandwidth 2GHz, in series with a 50ฮฉ resistor. The antenna is centered in a radiation box extending 12 ๐œ†max beyond the edges of the antenna. A PML region extends 13 ๐œ†max past the outside of the radiation box on all sides. The antenna is shown in Figure 2.4. 28 Figure 2.4: One face of the expontentially-tapered Vivaldi antenna used in this example The circuit configuration used for this analysis is shown in Figure 2.5. Figure 2.5: Circuit used for simulating Vivaldi antenna The antenna was discretized using tets with average edge length ๐œ†min /60 around the feed and to ๐œ†min /6 at the outer edges of the radiation box and PML region. This produced a total of 218K elements in the mesh. The system was simulated for 10K time-steps. The voltage and current across the port were measured after simulation with both port extraction and self-consistent methods. The normalized difference, or โ„’ 2 error, is used for comparison, calculated as ||๐‘ฅ1 โˆ’ ๐‘ฅ2 || 2 |๐‘ฅ| result = (2.2) ||๐‘ฅ2 || 2 Using (2.2), the โ„’ 2 error of voltage across the port edge of the methods is 6.4 ร— 10โˆ’12 , which is near the set solver tolerance of 10โˆ’12 . The port voltages are shown in Figure 2.6. After 29 applying a Fourier Transform to the time domain data, the resulting computed reflection coefficient found with the data from both methods is plotted in Figure 2.7. Figure 2.6: Comparison of port voltage of Vivaldi antenna for self-consistent (coupled) and port extraction methods 30 Figure 2.7: Reflection coefficient of a Vivaldi antenna for self-consistent (coupled) and port extracted solutions 2.3 Optimization of Linear Systems Circuit-Only Optimization The genetic algorithm described in section 1.6 is used to optimize phase of a linear circuit system [39]. Consider the generalized cost function (1.1), where a waveform has phase ๐‘“ (๐“”) and the desired phase is ๐‘‘. The cost function is minimized using linear circuit components to create a phase shift. An RC pair is selected to create the phase shift. This circuit layout is shown in Figure 2.8. 31 Figure 2.8: RC circuit used to create phase shift in circuit-only optimization This system will use single-objective optimization with a single component value as input. The value of the resistor is fixed to 50ฮฉ, and the capacitance is variable. ๐œ‹ A phase shift of is desired. The voltage source produces a sinusoid with frequency 3 1GHz. The time delay between maxima of the source and node voltage at the output can be used to measure the relative phase shift: ๐‘ก ๐‘š๐‘Ž๐‘ฅ1 โˆ’ ๐‘ก ๐‘š๐‘Ž๐‘ฅ2 ๐œƒ๐‘š๐‘’ ๐‘Ž๐‘ ๐‘ข๐‘Ÿ๐‘’ ๐‘‘ = โˆ— 2๐œ‹ (2.3) ๐‘‡ This can be checked against the analytically found result using the formula: ๐‘‹๐ถ ๐œƒ๐‘Ž๐‘›๐‘Ž๐‘™ ๐‘ฆ๐‘ก๐‘–๐‘ = ๐‘ก๐‘Ž๐‘› โˆ’1 ( ) (2.4) ๐‘… Using equation (2.4) with a resistance of 50ฮฉ and frequency of 1GHz, the resulting ๐œ‹ capacitance to find phase shift is approximately 1.84pF. 3 The capacitance input bounds are set as 0.1pF to 1nF. The cost function is represented using ๐“” as the single-input capacitance value, ๐‘“ (๐“”) as the calculated phase shift, and ๐‘‘ as ๐œ‹ . Therefore, minimizing (1.1) leads to a selection of capacitance which gives the output 3 ๐œ‹ waveform phase shift from the voltage source. 3 Given the phase is calculated in time domain using a discrete waveform, some amount of error is to be expected based on the time-step size. Simulating the sinusoid for approximately two periods with a time-step size of 9.1ps, the minimum phase resolution between two time-steps found using equation (2.4) is 5.7 ร— 10โˆ’2 radians (approximately 3.3โ—ฆ ). An error up to this minimum resolution may be expected even when the genetic algorithm converges. 32 This experiment is run with modified settings of the genetic algorithm being a population of 10, offspring per generation of 10, total generations set to 200, and mutation factor ๐œ‚ ๐‘€ set to 20. The optimized result is an input capacitance of 2.089pF and a resulting phase ๐œ‹ shift of 1.0282rad, which is 0.019rad (approximately 1.1โ—ฆ ) away from . Figure 2.9 displays 3 the source versus optimized output for this case, and Figure 2.10 shows the convergence of ๐œ‹ the output phase to . The optimization scheme converged near the analytic solution for 3 the given expected error. Figure 2.9: Waveforms from linear circuit-only optimization for phase 33 Figure 2.10: Convergence plot from linear circuit-only optimization for phase using larger time-step size Using a smaller step size and additionally increasing the total number of evaluations within the genetic algorithm will lower the difference between optimized and analytic results. This is tested by decreasing time-step size down to 0.091ps, corresponding to a minimum phase resolution of 5.7 ร— 10โˆ’4 radians. The number of generations is increased to 1000, and a terminating condition is set for cost function value below 10โˆ’3 . The scheme reached the terminating condition with the resulting convergence plot shown in Figure 2.11. 34 Figure 2.11: Convergence plot from linear circuit-only optimization for phase using finer time-step size Single-Objective Optimization of Linear EM-Circuit System Single-objective optimization paired with port extraction will be demonstrated using the Vivaldi antenna from section 2.2. The reflection coefficient shown in Figure 2.7 was calculated over 8GHz to 12GHz using a circuit with unit amplitude modified Gaussian voltage source and series 50ฮฉ resistance before the port nodes. Observing a wider band, the same circuit will be modified to calculate and plot the reflection coefficient from 1GHz to 21GHz. Let a desirable reflection coefficient threshold be established as -10dB. Consider an optimization scheme to widen the band for which the reflection coefficient satisfies the threshold. There are many possible approaches to do this. The approach selected here is to maximize the percentage of the band meeting the threshold. This objective only checks whether or not the value of the reflection coefficient is less than -10dB at each step 35 in frequency. The source impedance will be modified with an objective of increasing the percentage of the band with reflection coefficient value less than -10dB. The same model is used from section 2.2 but changing discretization around the feed area to account for ๐‘“max equal to 21GHz for a total of 316K mesh elements. Using an initial feeding circuit with impedance 50ฮฉ, the reflection coefficient over the 1GHz to 21GHz band is calculated and shown in Figure 2.12. For this initial circuit, 35.36% of the band has reflection coefficient less than -10dB. The circuit impedance will be modified to maximize this percentage. Figure 2.12: Initial Vivaldi antenna reflection coefficient over 1GHz to 21GHz The genetic algorithm settings modified were setting a population and offspring of 10 for 50 generations with ๐œ‚ ๐‘€ set to 20. The input components were a resistor and inductor in series to generate a complex impedance. The resulting feeding circuit is shown in Figure 2.13. 36 Figure 2.13: Circuit used for optimizing broadband Vivaldi antenna The component values form the input vector (๐“”). The percentage of the resulting reflection coefficient under -10dB over the band 1GHz to 21GHz represents the single objective ๐‘“ (๐“”). The desired objective value ๐‘‘ was set as 1, such that minimizing (1.1) converges the system to the optimal solution for this scheme. Input bounds were set as 0ฮฉ to 1000ฮฉ for the resistor and 1pH to 1nH for the inductor. The optimized result converged upon a source impedance of 75.67ฮฉ and 27.2pH, producing 44.06% of the reflection coefficient under -10dB compared to the 50ฮฉ unmodified circuit producing 35.36% over the 1GHz to 21GHz band. The convergence plot of this system is shown in Figure 2.14. The optimized coefficient is plotted against the unmodified data in Figure 2.15. To highlight the most impacted section of the band, Figure 2.16 displays the same plot for 8GHz to 16GHz. The optimized reflection coefficient increases compared to the unmodified data between roughly 9GHz and 12GHz, but remains under -10dB for this range; this behavior was possible since the objective allows for any part of the coefficient with value less than -10dB to count towards a fit solution. Within the 8GHz to 16GHz band of Figure 2.16, 47.37% of the original reflection coefficient satisfies the threshold for value less than -10dB and 79.52% of the optimized coefficient meets this threshold. 37 Figure 2.14: Convergence plot from linear broadband antenna reflection coefficient opti- mization 38 Figure 2.15: Comparison of reflection coefficient for unmodified and optimized broadband antenna feed 39 Figure 2.16: Highlighting impacted band in comparison of reflection coefficient for unmodified and optimized broadband antenna feed Multiple-Objective Optimization of Linear EM-Circuit System A log periodic dipole array is created and simulated to demonstrate multiple-objective optimization with port extraction on a linear EM-circuit system. The array consists of four elements with varying length and separation distance. The total lengths of the dipoles are 14.9cm, 11.9cm, 9.5cm, and 7.6cm respectively. The distances from the center of the first dipole with length 14.9cm to the centers of the remaining dipoles are likewise 0.9cm, 1.6cm, and 2.2cm respectively. Each element has a diameter of 0.4cm. A 1mm gap spans between each set of dipole flanges, across which a linear circuit system is connected. The EM system is shown in Figure 2.17, and the circuit layout used is shown in Figure 2.5. The initial circuit feeding each dipole consists of a modified Gaussian voltage source centered at 1.25GHz with bandwidth 750MHz, connected in series to a resistor with initial value 50ฮฉ. A radiation box is defined 12 ๐œ†max beyond each exterior face of the array. 40 A PML region is defined 13 ๐œ†max beyond the outside faces of the radiation box. Figure 2.17: Layout of the log periodic dipole array This system will be used to demonstrate optimization using multiple inputs and multiple objectives. The genetic algorithm objectives take value as the individual reflection coefficients, and the impedances within each feeding circuit serve as the inputs. Since multiple objectives are optimized, a solution set is produced rather than a single solution. This requires an additional method for determining the final solution from the set. The method chosen is the minimum sum of the four reflection coefficients. A model of the geometry is constructed. The system is discretized with tets having average edge length ๐œ†min /30 on the surfaces of the dipoles and ๐œ†min /6 within the radiation box and PML region, totalling 1.82M mesh elements. Prior to optimization, reference reflection coefficient plots, noted as unmodified, are found using these feeding circuits. They are shown in Figure 2.18. 41 Figure 2.18: Reflection coefficients for log periodic dipole array with unmodified sources The circuit impedance feeding each dipole is modified to optimize the reflection coefficient of each element. The genetic algorithm settings modified for this result were a set population of 80 with 45 offspring per generation and 75 total generations. The mutation coefficient ๐œ‚ ๐‘€ was set as 20. To minimize the cost function, the input vector ๐“” represents the value of the resistance within each circuit, ๐‘“ (๐“”) is set as a vector of the calculated reflection coefficients, and ๐‘‘ is equal to zero. The minimum reflection coefficient sum of the result set was found. This result is compared to the reference reflection coefficients and is shown in Figure 2.19. In this figure, unmodified and optimized parts are grouped in labelling using 1 โ‰ค ๐‘– โ‰ค 4 to represent S11, S22, S33, and S44. The optimized reflection coefficients have peaks at shifted frequencies from the original unmodified coefficient plots. This could be remedied by expanding the objective to prioritize maintaining the same frequency along with minimizing reflection coefficient. 42 Figure 2.19: Multi-objective optimization of reflection coefficients for log periodic dipole array 2.4 Application to SINR Optimization The Applebaum Criterion for SINR Maximization Consider desired and interfering signals incident upon a linear antenna array. The objective is to maximize signal to interference and noise ratio (SINR). Prior work has investigated a method for accomplishing this analytically [40, 41], in which the phases of the element sources are modified to nullify incident interference. This has been expanded upon with use of a genetic algorithm to vary the source phases [42, 43]. These cases have utilized isotropic radiating elements with digital phase shift. Subsequent definitions for this problem originate from these studies. We define that, for this linear antenna array with equal spacing and ๐‘ number of identical elements, the element source amplitudes and phases are defined within the 43 weighting vector ๐’˜ = |๐‘ค ๐‘š |๐‘’ ๐‘—๐œ™๐‘š (2.5) for 0 โ‰ค ๐‘š โ‰ค ๐‘ where ๐‘ค ๐‘š is the magnitude and ๐œ™ ๐‘š is the phase for each element ๐‘š. We likewise define the angle vector ๐’–(๐œƒ) as 1 ๏ฃฎ ๏ฃน ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๐‘—๐œ‹๐‘ ๐‘–๐‘›(๐œƒ) ๏ฃบ ๏ฃฏ ๐‘’ ๏ฃบ ๐’–(๐œƒ) = ๏ฃฏ๏ฃฏ .. ๏ฃบ (2.6) . ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ๐‘’ ๏ฃฏ (๐‘โˆ’1)๐‘—๐œ‹๐‘ ๐‘–๐‘›(๐œƒ) ๏ฃบ ๏ฃบ ๏ฃฐ ๏ฃป The array receives incident desired and interference signals along with element-wise noise, represented as ๐’™(๐‘ก) = ๐’™desired (๐‘ก) + ๐’™interference (๐‘ก) + ๐’™noise (๐‘ก) (2.7) The desired signal is narrow-band and centered at the frequency corresponding to the half-wavelength element spacing. We represent this desired signal as ๐’™ desired (๐‘ก) = ๐‘Ž ๐‘‘ (๐‘ก)๐’–(๐œƒ๐‘‘ ) (2.8) for ๐‘Ž ๐‘‘ (๐‘ก) being some time-varying signal and ๐’–(๐œƒ๐‘‘ ) from the definition of ๐’–(๐œƒ) above with incident angle ๐œƒ๐‘‘ . The interference ๐’™ i (๐‘ก) is likewise introduced in the form of (2.8) by time-varying signals ๐‘Ž ๐‘– (๐‘ก) with varying incident angles ๐œƒ๐‘– . The noise term ๐’™n (๐‘ก) is introduced as an element-wise Gaussian noise added to the system. We define a covariance matrix for these signals ฮฆ๐‘ž = ๐ธ{๐’™ ๐‘ž (๐‘ก)โˆ— ๐’™ ๐‘ž (๐‘ก)๐‘‡ } (2.9) in which ๐‘ž โˆˆ {desired, interference, noise} and ๐ธ{ยท} is the expectation operator. The power received from a signal can be calculated as ๐‘ƒ๐‘ž = ๐’˜ โ„‹ ฮฆ ๐‘ž ๐’˜ (2.10) where the superscript โ„‹ is the Hermitian operator, or conjugate transpose. 44 The SINR is then calculated as the power received from the desired signal compared to that of the undesired interference plus noise. For a desired signal with amplitude ๐‘Ž, we can define the value ๐‘Žยฏ 2 = ๐ธ{|๐‘Ž| 2 }. Then, the following objective, defined as the Applebaum criterion from [40], maximizes SINR with relation to the defined quantities: 2 ๐‘ƒdesired ๐’˜ ๐‘‡ ๐’–(๐œƒ๐‘‘ ) SINR = = ๐‘Žยฏ 2 (2.11) ๐‘ƒundesired ๐’˜ โ„‹ (ฮฆ๐‘– + ฮฆ๐‘› ) ๐’˜ With the assumption that the incident angles ๐œƒ๐‘ž of the desired and interfering signals are known yet the individual amplitudes are indistinguishable upon receipt, it is not possible to directly compute ๐‘Žยฏ 2 and ฮฆ๐‘– + ฮฆ๐‘› . Therefore, to define a cost function ๐ฝ(๐“”) for optimization, the following function has the same maximum as (2.11) when ๐“” is the weighting vector ๐’˜ [42]. 2 ๐’˜ ๐‘‡ ๐’–(๐œƒ๐‘‘ ) ๐ฝ(๐“”) = (2.12) ๐’˜ โ„‹ (ฮฆ๐‘– + ฮฆ๐‘› + ฮฆ๐‘‘ ) ๐’˜ With these definitions in place for the antenna array, it is possible to optimize the array system for SINR maximization using (2.12) as the objective. Creating a Port Extraction Case for SINR Maximization A novel approach to this problem can be formulated by using a simulated EM system and measuring phase shift from port voltages for optimization, compared to previous cases utilizing isotropic elements and digital phase shift and optimization not involving the simulation of an EM system. A linear dipole array will be simulated and optimized for this new approach. The array consists of five dipoles with lengths and separation spacing equal to a half-wavelength, corresponding to a center frequency of 1.375GHz. The dipoles have diameter 4mm. A 1mm gap is introduced between the flanges of each element, over which a circuit spans the connecting edge. This EM system is shown in Figure 2.20. 45 Figure 2.20: Layout of the linear dipole array The source phases are generated using an RC pair and calculated after simulation using equation (2.3) as done for the circuit-only phase optimization scheme in section 2.3. Within the RC pair, the resistance is fixed as 50ฮฉ and the capacitor is connected parallel to the output ports with its value varied. The resulting measured phase will be used as input in (2.12) by means of weighting array ๐’˜. The amplitude of the interference is defined as 20dB greater than that of the desired signal, and the element-wise noise is defined as 20dB less than the desired signal. A single interfering wave is incident upon the linear dipole array, with the incidence angle changing between 40โ—ฆ and โˆ’30โ—ฆ off broadside approximately every 50 time-steps. The desired signal is incident broadside of the array. Figure 2.21 roughly depicts this scenario, with a smaller desired signal incident broadside to the array and time-varying interference from off-broadside angles of the array. Both the desired and interfering signals are defined as modified Gaussian waves calculated using the system center frequency, bandwidth, and time-steps. 46 Figure 2.21: Incident desired and interfering waves on linear dipole system As done in [42], the voltage source amplitudes in the array are fixed with relative weightings to produce a Dolph-Chebyshev pattern with 20dB sidelobe suppression. The coefficients multiplying each element are calculated using Chebyshev polynomial coefficients for a fourth-order Chebyshev polynomial taken from [44]. The method for calculating weighting coefficients comes from [45]. The resulting coefficients are: ๏ฃฎ ๏ฃน ๏ฃฏ0.2956๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ0.6837๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ |๐’˜| = ๏ฃฏ๏ฃฏ 1.0 ๏ฃบ๏ฃบ (2.13) ๏ฃฏ ๏ฃบ ๏ฃฏ0.6837๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ0.2956๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฐ ๏ฃป with 1.0 corresponding to the center element, and 0.2956 to the end elements of the five-element array. SINR is then optimized using a cost function ๐ฝ(๐“”) equal to (2.12) with ๐‘“ (๐“”) being the weighting array ๐’˜ in section 2.4. An approximation of SINR will be made to show relative 47 effectiveness of the optimization scheme. The approximation is made using the power from interfering signals and noise, ๐‘ƒundesired , as the total power across a port in the array minus the calculated power of the defined desired signal. Equation (2.11) is then used to calculate the approximate SINR. Linear System SINR Maximization The linear dipole array is modeled and centered in a radiation box extending 12 ๐œ†max beyond each exterior face of the array. A PML region extends 13 ๐œ†max past the outer faces of the radiation box. The geometry is discretized using tets with an average edge length of ๐œ†min /30 on the dipoles and ๐œ†min /6 in the radiation box and PML region, leading to 202K total mesh elements. The case is simulated using a time-step size of 11ps. The cost function ๐ฝ(๐“”) is set as (2.12) with ๐“” taking value as the phases within the weighting array ๐’˜. To generate phase shifts in the array elements for ๐’˜, an RC pair is used and the phase shift is calculated from the measured port voltages using equation (2.3). Within the RC pair, the resistance values are fixed at 1ฮฉ and the capacitance values are varied between 70pF and 1nF. To counter voltage loss due to the addition of phase shifting elements, an inverting amplifier is added prior to the RC pair to increase the output voltage. The inverting amplifier consists of an ideal operational amplifier with positive terminal grounded, negative terminal preceded by a 50ฮฉ resistor, and an amplifying resistor which varies in value from 10ฮฉ to 1Mฮฉ. Additionally, to generate phase shift up to 180โ—ฆ , the number of RC pairs is varied between one and three. The voltage source has a time-step delay which varies from 0 to 150 time-steps as well. The resulting circuit for a single RC pair is shown in Figure 2.22. 48 Figure 2.22: Linear circuit used for SINR optimization with modified components high- lighted The set of optimization inputs which influence output phase and amplitude for ๐’˜ and subsequent SINR approximation are then the capacitance, the number of RC pairs, the amplifier resistance, and the source delay. The genetic algorithm settings are modified for a population and offspring of 15, 5 generations, and ๐œ‚ ๐‘€ equal to 5. Since the interference changes every 50 time-steps and the total signal is constantly changing in time, the genetic algorithm is reset at intervals of 50 time-steps to prevent convergence to just one interference incidence angle. Using (2.12) as the optimization objective and approximating SINR using the method described in 2.4, the system is optimized and SINR improvement over time is shown in Figure 2.23. The dips to 0dB SINR improvement correspond to the minimum result each time the genetic algorithm increments a generation and the drop in peaks every 50 time-steps correspond to the algorithm being reset. 49 Figure 2.23: Optimized SINR for linear dipole system with time-varying interference 2.5 Conclusion This chapter has demonstrated the use of transient port extraction to conduct linear EM- circuit analysis and optimization using a genetic algorithm. Transient port extraction was shown to produce solutions equivalent to that of self-consistent analysis to solver-tolerance precision. The optimization scheme was developed starting at the circuit-level, converging to analytic solutions for a phase-shifting circuit. Single-objective and multiple-objective optimization were demonstrated for improving reflection coefficient by modifying the feeding circuit. These cases built towards extending a previous investigation into SINR maximization, using an EM-circuit simulation instead of isotropic sources with digital phase shift. While testing the optimization scheme on linear systems is tenable with self-consistent solution methods, the adaption of this scheme to non-linear cases is straightforward with transient port extraction and unrealizable with self-consistent analysis. Chapter 3 details 50 the transition from linear to non-linear systems using the methods developed to this point. 51 CHAPTER 3 PORT PARAMETER EXTRACTION AND OPTIMIZATION OF NON-LINEAR SYSTEMS 3.1 Introduction With the developments made in chapter 2, the use of transient port extraction to analyze tightly-coupled systems of broadband EM devices fed by strongly non-linear circuits is now investigated within this chapter. Optimization of these systems is unrealizable with typical self-consistent solution methods yet tenable with port extraction. Since the EM model is represented in terms of port responses and is circuit-agnostic, the analysis and optimization schemes developed in the previous chapter can reuse the existing port responses with circuits now incorporating non-linear components. Furthermore, the use of non-linear components may necessitate solving the system with a smaller time-step size. An advantage of port extraction is that the extracted response of the EM system can be interpolated for use with a smaller step size. This chapter compares the results of port extraction for a monopole fed by a non-linear amplifier to that of a self-consistent solution. The non-linear amplifier is then applied to the broadband Vivaldi antenna from section 2.2 and the system is optimized to reduce amplifier distortion. Next, the SINR case from the previous chapter is modified to include non-linear circuit components and then optimized. 3.2 Optimization of Non-Linear Systems Incorporating non-linear circuit elements requires the use of an iterative method for convergence of the non-linear voltage-current relation at each time-step. These non-linear components include devices such as diodes and transistors. The linear systems developed in the previous chapter can be augmented with non-linear circuit components to create non-linear cases. Simulating the EM-circuit system in a self-consistent, fully-coupled 52 manner becomes computationally expensive due to repeatedly inverting the coupled EM-circuit matrix within the non-linear iterative process. Repeating this process many times within an optimization scheme is then unrealizable. Utilizing port extraction to create a lower-order circuit-agnostic model in which the EM system is characterized by port responses rather than a large mesh allows the non-linear MNA solve and internal iterative method to be isolated from the EM system, as depicted in the right-hand side of Figure 1.1. Verifying Port-Extraction to Self-Consistent Method for Non-Linear Feed An exponential amplifier is used as the feeding circuit in a non-linear EM-circuit network for comparison of self-consistent and transient port extraction solutions. This circuit has amplification equal to ๐‘‰๐‘–๐‘› ๐‘‰out = โˆ’๐‘… ๐ผsat ๐‘’ ๐‘‰๐‘‡ where ๐‘… is the value of the resistor in parallel with the op amp, ๐ผsat is the saturation current of the diode, ๐‘‰๐‘–๐‘› is the voltage into the diode, and ๐‘‰๐‘‡ is the thermal voltage. The resulting non-linear circuit is shown in Figure 3.1. Figure 3.1: Non-linear amplifying circuit used for comparison between port extraction and self-consistent solutions The EM geometry used will be the strip monopole above a finite ground plane, previously used in section 2.2 and shown in Figure 2.1. 53 The source circuit contains an exponential amplifier which maintains a non-linear relation between amplified voltage and current. To test these methods, the circuit shown in Figure 3.1 is used. The voltage source in the circuit is set to a sinusoid with frequency 1GHz, matching the center frequency of the previous strip monopole case. The circuit values used are the unit amplitude sinusoidal voltage source; a diode with saturation current 1fA, ideality factor 1.6, and thermal voltage 25.85mV; an operational amplifier with positive terminal grounded and negative terminal connected to the diode; a resistor with value 100kฮฉ in parallel with the operational amplifier negative terminal and output; and a 1ฮฉ resistor between operational amplifier output and the port nodes. Simulating this EM-circuit system for 1000 time-steps with a step size of 2.62ps yields an output voltage exponentially amplified and inverted to the input. The results of port extraction and a self-consistent, fully coupled solve are shown in Figure 3.2. The โ„’ 2 error of the data shown is 3.8 ร— 10โˆ’12 , near the set solver tolerance of 10โˆ’12 . Figure 3.2: Comparison of port extraction and self-consistent (denoted coupled) solution for output voltage of strip monopole fed by non-linear amplifier circuit 54 Non-Linear Analysis of a Broadband Vivaldi Antenna The Vivaldi antenna optimized with linear components in section 2.3 is now used to demonstrate optimization of a broadband EM system fed by a non-linear circuit. This is the optimization analysis enabled by port extraction which is otherwise untenable using a fully coupled simulation method. The exponential amplifier in Figure 3.1 will be modified for use with the Vivaldi antenna. The voltage source is defined as a modified Gaussian pulse with unit amplitude, center frequency of 11GHz, and bandwidth of 10GHz. The diode is defined with a saturation current of 1fA, ideality factor of 1.6, and thermal voltage of 26mV. The initial amplifying resistor is set to 1Mฮฉ, and the output resistor at 1ฮฉ. The non-optimized, or "naive", amplified signal of this EM-circuit system is compared to that from the same voltage source with only a 50ฮฉ series resistance, noted as unamplified, in Figure 3.3. Figure 3.3: Time-domain port voltage from non-optimized ("naive") exponential amplifier feeding broadband Vivaldi antenna compared to unamplified signal 55 The port voltage is amplified with strong distortion present. A Fourier transform of these signals is taken to better understand the distortion. The same port voltages are shown in frequency domain in Figure 3.4. Figure 3.4: Frequency-domain comparison of exponential amplifier feeding broadband Vivaldi antenna versus unamplified signal An optimization scheme is developed for this non-linear system. The objective is to minimize distortion of the amplifier as applied to the Vivaldi antenna. One potential solution includes adding a band-pass filter to the amplifying circuit. A passive RC band- pass filter is selected for simplicity. This circuit configuration is shown in Figure 3.5 with modified components highlighted. 56 Figure 3.5: Initial state of exponential amplifier circuit with band-pass filter feeding broadband Vivaldi antenna For this optimization scheme, the vector of inputs ๐“” includes the amplifying resistance, the low-pass filter capacitance, the high-pass filter capacitance, and the numbers of low-pass and high-pass filters cascaded before the port. The resistors following the output of the operational amplifier and in the filter are fixed at 1ฮฉ. The amplifying resistance is varied from 100ฮฉ to 10Mฮฉ, both filter capacitances are individually varied from 1pF to 100pF, and the number of low-pass and high-pass RC pairs individually varied from one to ten. Each low-pass filter uses the same capacitance and resistance, and likewise for the high-pass filters respectively. The cost function is defined as the relative absolute difference between the Fourier transforms of the measured port voltage of the optimized filtered amplifier circuit and the unamplified port voltage. The optimization scheme developed uses a population of 150, offspring of 50, and 30 total generations with mutation coefficient ๐œ‚ ๐‘€ equal to 20. The result of this optimization is shown below, with time-domain comparison in Figure 3.6 and frequency-domain magnitude comparison in 3.7. After optimization, the resulting circuit contains two identical high-pass RC pairs with values 1ฮฉ and 5.4pF, and nine identical low-pass RC pairs with values 1ฮฉ and 6.7pF. The relative absolute difference to the unamplified frequency- domain waveform is decreased from 0.761 for the initial naive amplifier circuit to 0.161 for the optimized filtered amplifier circuit. 57 Figure 3.6: Time-domain comparison of unamplified, naive amplified, and optimized amplified port voltages of broadband Vivaldi antenna 58 Figure 3.7: Frequency-domain comparison of unamplified, naive amplified, and optimized amplified port voltages of broadband Vivaldi antenna Non-Linear SINR Optimization The linear SINR optimization set-up from section 2.4 will be combined with the non- linear circuit in Figure 3.1 to create a non-linear optimization scheme. The changes made from the linear case include using the non-linear exponential amplifier circuit as the source circuit for each element as well as varying the thermal voltage of the diode by means of changing ambient temperature. Practically, this could mean housing the diode in a thermal chamber to control the temperature around the component, or alternatively considered as an approach for an inverse scheme for a temperature sensor to assess ambient conditions. Since port extraction creates a circuit-agnostic EM model, modifying the source circuits to be non-linear is a simple alteration to the optimization scheme. All other components of the case are unchanged. SINR is optimized utilizing the same cost function as section 2.4. The feeding circuit selected is shown in Figure 3.8. Although this amplifier was 59 demonstrated to have strong amplitude distortion, this case will focus on output phase shift from the non-linear amplifier. Varying RC pairs are added to better control phase shift. The vector of inputs ๐“” consists of a voltage source delay of 0 to 150 time-steps, an amplifying resistance between 100ฮฉ and 1Mฮฉ, capacitance between 30pF and 80pF, and diode temperature between 285K and 315K for determining thermal voltage. Figure 3.8: Non-linear exponential amplifier circuit used for SINR optimization The resulting EM and non-linear circuit system is optimized using the same genetic algorithm inputs and scheme as section 2.4. The resulting approximated SINR improvement is shown in Figure 3.9. Every 50 time generations when the interference changes, the best solution per generation within the genetic algorithm dips back towards 0dB improvement, meaning interference is poorly mitigated. Towards the end of each interference window, SINR reaches between 15 and 20dB improvement. Considering computational cost, the post-extraction circuit solve for this system takes an average of 243๐œ‡s, whereas a single self-consistent forward solve takes approximately 4.5s. This improvement in computational cost is about five orders of magnitude. 60 Figure 3.9: Optimized SINR for non-linear circuits feeding dipole system with incident time-varying interference 3.3 Conclusion Optimization with port extraction for non-linear EM-circuit systems has been demon- strated by decreasing distortion from a strongly non-linear amplifier as well as adapting the SINR case from chapter 2 for a non-linear feeding circuit. A comparison was made between port extracted and self-consistent solutions for a monopole fed by a non-linear amplifier, showing identical solutions between the methods to solver-tolerance precision. The non-linear amplifier was then adapted and optimized to reduce distortion in port voltage when feeding a broadband Vivaldi antenna. The relative absolute difference of the non-linear amplifying circuit was reduced from 0.761 to 0.161 when compared to the unamplified port voltage. The SINR optimization scheme developed in chapter 2 was adapted for use in a non-linear system by including the non-linear amplifier within the dipole array circuit system, with a variety of inputs being optimized including diode ther- 61 mal temperature. The results showed success in mitigating the time-varying interference incident upon the dipole system. 62 CHAPTER 4 CONCLUSION AND FUTURE WORK The focus of this work was to investigate using transient port extraction to create a reduced-order EM-circuit model amenable to optimization. Optimization using port extraction was demonstrated by first developing an optimization scheme for linear EM- circuit systems. Linear circuit-only optimization was shown to converge near to the desired analytic solution, subject to relevant solving tolerances such as step size and convergence tolerance. This was then applied to reflection coefficient minimization for a single-objective optimization case using a broadband Vivaldi antenna, and a multi-objective optimization case using a log-periodic dipole array. These linear methods were then combined to demonstrate SINR maximization using Applebaumโ€™s criterion for a linear dipole array with time-varying incident interference. Next, the linear system optimization methods developed were applied to non-linear systems. The result from transient port extraction was compared to that of a self-consistent method for a strip monopole above a finite ground plane fed by a non-linear circuit. To demonstrate the use of transient port extraction for optimization of a tightly-coupled broadband EM system and strongly non-linear feeding circuit, a non-linear amplifying circuit feeding a broadband Vivaldi antenna was optimized to reduce distortion of the port voltage compared to the source. The linear feeding circuits in the SINR case were modified to include non-linear components, and the optimization scheme was repeated with a wider range of variables including diode thermal voltage. The results showed maximization of SINR using non-linear circuits to feed the dipole array and mitigate time-varying incident interference. The work within this thesis has demonstrated the use of transient port parameter extrac- tion to circumvent the computational cost of optimization of these systems with traditional, self-consistent coupled analysis. In one scenario included, the cost for optimization after 63 port extraction of the EM system was five orders of magnitude more efficient than the equivalent self-consistent scheme. This work has potential for further investigation. The current scope of non-linear circuit components was limited to a diode; utilizing other MNA stamps for non-linear components such as transistors opens consideration for more complex non-linear circuits, and subsequently further optimization possibilities. The multi-objective optimization scheme developed focused on an array with similar elements, near-identical circuits, and a common objective in minimizing reflection coefficient. This could be expanded to examine and optimize a system containing different EM devices with varying circuit configurations and objectives. Furthermore, the SINR optimization case utilized Applebaumโ€™s criterion to mitigate time-varying incident interference by means of circuit modification. Using transient port extraction for optimization of tightly-coupled broadband EM devices fed by strongly non-linear circuits could be expanded to other SINR maximization or EM-circuit optimization schemes. 64 BIBLIOGRAPHY [1] K. Li, A.M. Tassoudji, S.Y. Poh, M. Tsuk, R.T. Shin, and J.A. Kong. Fd-td analysis of elec- tromagnetic radiation from modules-on-backplane configurations. IEEE Transactions on Electromagnetic Compatibility, 37(3):326โ€“332, 1995. [2] Er-Ping Li. Electrical Modeling and Design for 3D System Integration. John Wiley & Sons, Ltd, 2012. [3] Yong Wang, D. Gope, V. Jandhyala, and C. . R. Shi. Generalized kirchoffโ€™s current and voltage law formulation for coupled circuit-electromagnetic simulation with surface integral equations. IEEE Transactions on Microwave Theory and Techniques, 52(7):1673โ€“1682, 2004. [4] Paule Charland, Ellen El-Khatib, and John Wolters. The use of deconvolution and total least squares in recovering a radiation detector line spread function. Medical Physics, 25(2):152โ€“160, 1998. [5] Scott Oโ€™Connor, Stephen Hughey, Dan Dault, Andrew J. Pray, Jorge M. Villa-Giron, and Balasubramaniam Shanker. A novel port/network parameter extraction technique for coupling circuits with full-wave time-domain integral equation solvers. IEEE Transactions on Microwave Theory and Techniques, 67(2):553โ€“564, 2019. [6] Omkar H. Ramachandran, Scott Oโ€™Connor, Zane D. Crawford, Leo C. Kempel, and B. Shanker. Port parameter extraction-based self-consistent coupled em-circuit fem solvers. IEEE Transactions on Components, Packaging and Manufacturing Technology, 12(6):1040โ€“1048, 2022. [7] Richard Courant et al. Variational methods for the solution of problems of equilibrium and vibrations. Lecture notes in pure and applied mathematics, pages 1โ€“1, 1994. [8] PL Arlett, AK Bahrani, and OC Zienkiewicz. Application of finite elements to the solution of helmholtzโ€™s equation. In Proceedings of the Institution of Electrical Engineers, volume 115, pages 1762โ€“1766. IET, 1968. [9] P Silvester. A general high-order finite-element analysis program waveguide. IEEE Transactions on Microwave Theory and Techniques, 17(4):204โ€“210, 1969. [10] Bo-nan Jiang, Jie Wu, and Louis A Povinelli. The origin of spurious solutions in computational electromagnetics. Journal of computational physics, 125(1):104โ€“123, 1996. [11] M Hara, T Wada, T Fukasawa, and F Kikuchi. A three dimensional analysis of rf electromagnetic fields by the finite element method. IEEE Transactions on Magnetics, 19(6):2417โ€“2420, 1983. [12] Jie Wu and Bo-nan Jiang. A least-squares finite element method for electromagnetic scattering problems. 1996. 65 [13] BM Azizur Rahman and J Brian Davies. Penalty function improvement of waveguide solution by finite elements. IEEE Transactions on Microwave Theory and Techniques, 32(8):922โ€“928, 1984. [14] Jean-Claude Nรฉdรฉlec. Mixed finite elements in R3 . Numerische Mathematik, 35(3):315โ€“ 341, 1980. [15] Zoltan J Cendes. Vector finite elements for electromagnetic field computation. IEEE Transactions on Magnetics, 27(5):3958โ€“3966, 1991. [16] Hassler Whitney. Geometric integration theory. Courier Corporation, 2012. [17] Rui Wang and Jian-Ming Jin. A symmetric electromagnetic-circuit simulator based on the extended time-domain finite element method. IEEE transactions on microwave theory and techniques, 56(12):2875โ€“2884, 2008. [18] Zheng Lou and Jian-Ming Jin. Modeling and simulation of broad-band antennas using the time-domain finite element method. IEEE transactions on antennas and propagation, 53(12):4099โ€“4110, 2005. [19] Dan Jiao, A Arif Ergin, Balasubramaniam Shanker, Eric Michielssen, and Jian-Ming Jin. A fast higher-order time-domain finite element-boundary integral method for 3-d electromagnetic scattering analysis. IEEE Transactions on Antennas and Propagation, 50(9):1192โ€“1202, 2002. [20] Su Yan and Jian-Ming Jin. A fully coupled nonlinear scheme for time-domain modeling of high-power microwave air breakdown. IEEE Transactions on Microwave Theory and Techniques, 64(9):2718โ€“2729, 2016. [21] Jin-Fa Lee, R. Lee, and A. Cangellaris. Time-domain finite-element methods. 45(3):430โ€“ 442. [22] A. Bossavit. Whitney forms: a class of finite elements for three-dimensional computa- tions in electromagnetism. 135(8):493โ€“500. [23] A. Bossavit and L. Kettunen. Yee-like schemes on a tetrahedral mesh, with diagonal lumping. 12(1-2):129โ€“142. [24] R Hiptmair. Finite elements in computational. Acta Numerica 2002: Volume 11, (11):237, 2002. [25] RN Rieben, GH Rodrigue, and DA White. A high order mixed vector finite element method for solving the time dependent maxwell equations on unstructured grids. Journal of Computational Physics, 204(2):490โ€“519, 2005. [26] Chung-Wen Ho, A. Ruehli, and P. Brennan. The modified nodal approach to network analysis. IEEE Transactions on Circuits and Systems, 22(6):504โ€“509, 1975. [27] W. McCalla and D. Pederson. Elements of computer-aided circuit analysis. IEEE Transactions on Circuit Theory, 18(1):14โ€“26, 1971. 66 [28] B. Gear. Backward differentiation formulas. Scholarpedia, 2(8):3162, 2007. revision #91024. [29] M.M. Hassoun and Pen-Min Lin. A formulation method for including ideal operational amplifiers in modified nodal analysis. In Proceedings of 40th Midwest Symposium on Circuits and Systems. Dedicated to the Memory of Professor Mac Van Valkenburg, volume 2, pages 1161โ€“1164 vol.2, 1997. [30] Jean-Pierre Berenger. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2):185โ€“200, 1994. [31] T. Rylander and Jian-Ming Jin. Perfectly matched layer in three dimensions for the time-domain finite element method applied to radiation problems. IEEE Transactions on Antennas and Propagation, 53(4):1489โ€“1499, 2005. [32] Nathan Mortimore Newmark. Computation of dynamic structural response in the range approaching failure. Department of Civil Engineering, University of Illinois. [33] W. L Wood. Practical time-stepping schemes. Oxford [England] : Clarendon Press ; New York : Oxford University Press, 1990. [34] O. C. Zienkiewicz. A new look at the newmark, houbolt and other time stepping formulas. a weighted residual approach. Earthquake Engineering & Structural Dynamics, 5(4):413โ€“418, 1977. [35] Zane Crawford, Jie Li, Andrew Christlieb, and Balasubramaniam Shanker. Uncon- ditionally stable time stepping method for mixed finite element maxwell solvers. Progress In Electromagnetics Research C, 103:17โ€“30, 01 2020. [36] J. Blank and K. Deb. pymoo: Multi-objective optimization in python. IEEE Access, 8:89497โ€“89509, 2020. [37] Kalyanmoy Deb and Debayan Deb. Analysing mutation schemes for real-parameter genetic algorithms. Int. J. Artif. Intell. Soft Comput., 4(1):1โ€“28, feb 2014. [38] Kalyanmoy Deb, Karthik Sindhya, and Tatsuya Okabe. Self-adaptive simulated binary crossover for real-parameter optimization. In Proceedings of the 9th Annual Conference on Genetic and Evolutionary Computation, GECCO โ€™07, page 1187โ€“1194, New York, NY, USA, 2007. Association for Computing Machinery. [39] Randy Haupt and Douglas Werner. Genetic Algorithms in Electromagnetics. John Wiley Sons, Ltd, 2007. [40] S. Applebaum. Adaptive arrays. IEEE Transactions on Antennas and Propagation, 24(5):585โ€“598, 1976. [41] R.T. Compton. The effect of random steering vector errors in the applebaum adaptive array. IEEE Transactions on Aerospace and Electronic Systems, AES-18(4):392โ€“400, 1982. 67 [42] D.S. Weile and E. Michielssen. The control of adaptive antenna arrays with genetic al- gorithms using dominance and diploidy. IEEE Transactions on Antennas and Propagation, 49(10):1424โ€“1433, 2001. [43] A. Massa, M. Donelli, F.G.B. De Natale, S. Caorsi, and A. Lommi. Planar antenna array control with genetic algorithms and adaptive array theory. IEEE Transactions on Antennas and Propagation, 52(11):2919โ€“2924, 2004. [44] Milton Abramowitz. Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables,. Dover Publications, Inc., USA, 1974. [45] Constantine A. Balanis. Antenna Theory: Analysis and Design. Wiley-Interscience, USA, 2005. [46] Jianming Jin and Douglas J. Riley. Finite Element Analysis of Antennas and Arrays. Wiley-IEEE Press, 2009. 68