Passive reduced-order modeling of electromagnetic systems

Computer methods in applied mechanics and engineering E;omput.-l~letlaoas AplSI.-rv]efh.~ngrg. -lbx) (-1"999)"345-~3B~

Passive reduced-order modeling of electromagnetic systems Andreas C. Cangellaris a'*, Li Zhao b aECE Department, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA bECE Department, University of Arizona, Tucson, AZ 85721, USA

Received 15 January 1998

Abstract Rapid electromagnetic field simulation is now being recognized as a critical component of next-generationcomputer-aided design frameworks that will be used for performanceevaluationand design of the informationprocessing and communicationsystems of the 21st century. Except for very simple systems, computer simulation of electromagnetic interactions is hindered by the very large number of degrees of freedom involved in the discrete model. One potentiallyuseful approach to overcomingthis computationalbottleneckis model order reduction, where parts of the electromagneticmodel are replaced by models of substantiallylower order, yet capable of capturing the electromagnetic behavior of the original subsystemswith sufficientengineeringaccuracy. Possible approaches to the developmentof such model order reduction techniquesare the subject of this paper. The potentialapplicationsand the benefitsfrom such electromagneticmodel order reductionare discussed. The importantissue of passivityof the generated reduced model is considered, and a set of constraintson the state representationof the discrete electromagneticboundary-valueproblem are identifiedfor the reduced order model to be passive. The paper concludeswith a presentationof several numericalexamples from the applicationof model order reductiontechniquesto the analysis of electromagnetic waveguides and antennas. © 1999 Published by Elsevier Science S.A. All rights reserved.

1. Introduction Rapid advances in device switching speeds, combined with the explosive growth of the wireless computing and communication market with its demands for low-voltage, low-power, tightly-integrated, mixed analogdigital Circuits, have brought the need for electromagnetic computer-aided design (EM-CAD) capability to the forefront of modem integrated circuit design. Clock frequencies in the GHz regime, with signal edge speeds less than 50 ps at the chip boundary, are expected from fast processors in the near future. For such switching speeds and clock frequencies, transmission line effects need be taken into account for proper timing analysis, while accurate prediction of displacement currents due to unbalanced interconnects as well as quantification of direct radiation from drivers are necessary for module electromagnetic compatibility assessment. The number of inputs/outputs (I/©s) per chip is climbing to several hundreds, and the prediction of the disturbance of groundand power-plane reference voltage and its impact on false switching is becoming a formidable task as several tens or even hundreds of the I / O s switch simultaneously at the aforementioned speeds. Accurate noise prediction and containment is becoming a major issue for the low-voltage designs that the electronics industry is trying to effect at a consumer market-driven low cost. Such low:voltage designs ( < 1.5 V supplies) require a very tight noise budget, which needs to be effected without over-hardening of the design in order to meet the low-cost demand. Furthermore, market competitiveness implies low cost and very short design cycles. Consequently~, rapid prototyping is essential, and can be effected only by means of a sophisticated electronic CAD environment equipped with rapid electromagnetic modeling capability.

* Corresponding author. 0045-7825/99/$ - see front matter © 1999 Published by Elsevier Science S.A. All rights reserved. PII: S0045-7825(98)001 62-5

346

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

Despite significant advances in electromagnetic modeling and simulation methodologies and computer tool sophistication, EM-CAD tools do not exhibit yet the efficiency needed for module/system design and subsequent design optimization. A brief assessment of state-of-the-art, commercially available EM-CAD tools reveals that their capabilities are currently limited to individual component design and circuits of rather low complexity. However, the compact integration that characterizes the new generation of electronic products requires a fairly rigorous modeling of the interactions between several components for proper design. Except for very simple systems, the distributed linear and nonlinear electromagnetic analysis required for the design and optimization of state-of-the-art and future high-performance information processing and communication systems is hindered by the large number of degrees of freedom involved in the electromagnetic model. Thus, the direct use of time-domain nonlinear differential-equation integration, as implemented in SPICE-like simulators enhanced with transmission-line modeling or even combined with rigorous FDTD-based Maxwell's equations solvers [ 1,2], is prohibitive for circuits of the complexity encountered in realistic integrated electronic systems. Even standard SPICE-like nonlinear circuit simulators for lumped nonlinear circuit analysis cannot handle the complexity of the circuits resulting from the modeling of the interconnects and power/ground distribution networks using lumped inductors, capacitors and resistors. The circuit simulation community has responding to this simulation challenge with the development of model order reduction techniques that allow the replacement of large linear portions of the circuit with a substantially smaller model that approximates sufficiently its external behavior [3-12]. This way, the size of the problem left for the nonlinear simulator to solve becomes substantially smaller, and thus problems of sizes out of the reach of conventional simulators can now be tackled. More recently, model order reduction has been applied to distributed electromagnetic systems [13-17]. More specifically, Pad~ approximations to the response of electromagnetic systems are used to effect the so-called fast frequency sweep, where the system response over a broad range of frequencies is obtained at approximately the cost of factoring the approximating matrix at only a few frequency points. For example, it was demonstrated in [17] that the matrix resulting from the finite difference approximation of Maxwell's equations needed to be factored only at a single frequency when the Lanczos algorithm was used for the generation of the Pad6 approximation. In addition to rapid broadband calculation of the frequency response of an individual electromagnetic component or a specific functional block, the development of macromodels of electromagnetic multiports is essential for design-driven simulation of large electromagnetic systems. When macromodels for multiports are connected together, it is important to keep in mind the fact that interconnections of stable systems may not necessarily result in stable systems; however, interconnections of passive circuits always result in systems that are passive and, hence, (asymptotically) stable [18]. This, then, implies that it is not enough for an electromagnetic macromodel to be stable. What matters, when this macromodel is to be connected with other functional blocks, is for the macromodel to be passive. In this paper, model order reduction methodologies are presented for semi-discrete electromagnetic systems obtained from the spatial discretization of the hyperbolic system of Maxwell's equations. It is shown how the properties of the original semi-discrete system impact the passivity of the developed reduced-order model, and specific criteria for developing passive reduced-order models are established. The presentation concludes with examples from the application of model order reduction to the rapid electromagnetic analysis of waveguide components and antennas.

2. The semi-discrete electromagnetic model While a variety of approaches exist for the spatial discretization of Maxwell's equations, the semi-discrete approximation effected through the use of Yee's lattice [19] is chosen for the purposes of this paper. The reason for this choice will become apparent when the passivity of the semi-discrete model is considered later in the paper. A uniform, rectangular lattice is assumed, defined by equally spaced nodes along the three axes of a Cartesian coordinate system: I along x, J along y and K along z. The total number of nodes in the grid is N = I- J . K. With the definitions U = 1, V= L W= I . J, the nth electric node, corresponding to node (i, j, k) in the grid, is given by

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 3 4 5 - 3 5 8

n= 1 +(i-

1)U+(j-

1)V+(k-

1)W

347

(1)

where i = 1,2 . . . . . /, j = 1,2 . . . . . J, k = 1,2 . . . . . K, and n = 1,2 . . . . . N. It is assumed that the media are linear, passive and time-independent. Thus, Maxwell's curl equations in the Laplace domain assume the form V X E = -stxH

(2)

V × H = s e e + o'E + J~

(3)

where the electric permittivity, e, electric conductivity o-, and magnetic permeability, /z, are, in general, position dependent. The dependence of the electric and magnetic field vectors E and H, as well as the imposed source current density J : on position and the Laplace variable s is suppressed for simplicity. In order to cast the semi-discrete form of Maxwell's equations in a matrix form, we begin with the definition of the following two vectors of discrete unknowns E = [E x, Ey, Ez] +

(4)

H

(5)

=

[Hx,H,, H~] r

where E x is a vector of length N, containing the N E x values on the grid, with similar definitions for the remaining five vectors. Using the definitions in (4), (5) and writing the curl operator in its matrix form

I0 Vx =

-0/0 010y /]

O/Oz

0

- O/ax

-O/ay

O/Ox

o

6)

it is straightforward to show that the semi-discrete form of (2), (3) can be written as -A v

-A T

0

T

AT

--AT

"E = - s D h H

[0 A AA] A w

0

-A v

Au

-

.

.H=sD~E+DoE+J~

(7)

(8)

In the above equations, the matrices A u , A ~ , A w are sparse with only two bands having nonzero elements: one band is along the diagonal with all values equal to 1, and the second band at a distance of U, V, W, respectively, to the left of the diagonal with all values equal to - 1 . The matrices De, D n, D~ are diagonal matrices with elements dependent on the electromagnetic properties of the media and the grid size. The system of (7), (8) may be cast in a compact form by defining the vector of state variables X = [H, E] T, and the 3N × 3N matrix, P,

[0 awa] -A

p =

-Aw A~

0 -A u

v

(9)

In addition, since the objective is to develop reduced-order macromodels for multiport electromagnetic systems, the source notation is slightly modified to include the imposed currents used for the excitation of the ports. For this purpose, it is assumed that the electromagnetic system under consideration has p ports, each coinciding with one electric field node, and a constant matrix B of dimension 6N / p is introduced, with nonzero elements only in its bottom 3N rows associated with the electric field nodes in the state vector X. The specific values of these elements will depend, in general, on the source distribution and numerical grid characteristics. Using the vector U(s) to denote the Laplace transforms of the current source waveforms at the p ports, the discrete source term may be cast in the form B U ( s ) . With these definitions, the resulting compact form of (7), (8) is

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

348

E° ro]X=

oO]x[:

o°]x

or, in a yet more compact form, (11)

(G + s C ) X = - B U ( s )

where

Because of the assumed passivity of the media, the matrices D e, D n are (symmetric) positive definite and D~ is (symmetric) non-negative definite. Consequently, C is also (symmetric) positive definite. Next, we introduce a frequency shift, s = s o + s'. Through a straightforward matrix manipulation, (11 ) is cast in the form (I

-

s'a)x

(13)

= RV(s')

where A = -(G + soC)-'C,

R = -(G + soC)-lB

(14)

Defining a desired output vector as Y(s) = FX(s), where F is a selector matrix, we have Y(s) = F ( I - s ' A ) - IR U (s)

(15)

For the case of multiports, the number of outputs is the same with the number of inputs. More specifically, with the selection of the current density as the excitation at each port, the electric field vector (or, equivalently, the voltage at the port, if the definition of the voltage makes sense) becomes the output quantity. For such cases, the selector matrix F is a p × 6N matrix, with constant, non-zero entries at the right half of each of its rows, corresponding to the electric field unknowns in the state vector X. In particular, except for a scaling constant factor associated with the proper definition of the observed output quantity at the port, the selector matrix F is simply the transpose of the source matrix B defined earlier. Thus, the electromagnetic characterization of mulfiports is effected in terms of a transfer-function matrix, H(s'), which, in view of (15), is given by H ( s ' ) = B T ( I -- s ' A ) - I R

(16)

From (15) it is clear that the eigenvalues of A correspond to the poles of the transfer function matrix H(s'). The objective of model order reduction techniques is to generate rational function approximations (Pad6 approximations) to the elements of H ( s ' ) of order much smaller than 6N. At this point, it is appropriate to consider this approximation objective carefully, in view of the fact that what is of primary interest in the context of this discussion is the model order reduction of distributed electromagnetic systems. Continuous electromagnetic systems are, of course, of infinite order; hence, the number of poles required for the exact representation of their delay properties is infinite. However, when dealing with a numerical approximation of Maxwell's equations, what is of interest is the reduced-order modeling of the resulting discrete electromagnetic system. These discrete systems have finite order, and their transfer functions have finite bandwidth which is controlled by the grid size and the order of approximation of the spatial derivatives. Consequently, the pole-residue representation of the transfer function of the approximated electromagnetic system will be of finite order. As already pointed out in the Introduction, the Asymptotic Waveform Evaluation (AWE) technique was the first one to be applied for such a purpose, for both circuits containing lumped elements as well as distributed electromagnetic circuits. However, because of the poor conditioning of the matrix resulting from the momentmatching phase of the algorithm, the originally anticipated increase in approximation accuracy as the order of the Pad6 approximation was increased is not achieved. As a result, only a limited number of poles (typically less than ten) can be extracted accurately using the original AWE algorithm. Several enhancements such as complex frequency hopping [5] and multipoint Pad6 [8] have been proposed and used to increase the order of the Pad6 approximation generated by the AWE algorithm.

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

349

The aforementioned limitation of the AWE process motivated research in alternative methods for the generation of Pad6 approximation of high orders [20,21,10,22]. In the following, we review two of the most successful algorithms that have resulted from these research efforts.

3. The Block Arnoidi algorithm For simplicity, let M denote the total number of state variables. Once again, the number of ports in the multiport which is to be reduced is denoted by p. The Block Arnoldi algorithm reduces the system matrix A in (14) to a q × q block upper Hessenberg matrix,/-/q such that

AV=VHq,

vTv=I

(17)

For the purposes of model order reduction it is desirable to have q << M. The M x q matrix V forms an orthonormal basis for the Krylov subspace

Kr(A, R, Ft) -~ colsp[R, AR, A 2R . . . . . A~R]

(18)

where ~ = LqlpJ with the operator L'J understood as the truncation below to the nearest integer. In order to explain how model order reduction is effected through the Block Arnoldi algorithm, consider the change of variables X = VX

(19)

where X is the reduced state vector of length q. This change of variables constrains the original state vector X to the q-dimensional subspace spanned by the columns of the matrix V. Substituting (19) into (13) and multiplying on the left by V T we obtain

( v T v - - s ' V T A V ) ~ = VTRU(s ') which, in view of (17), yields t

^

(I - s Hq)X = vTRU(s ')

(20)

Finally, the response of the reduced system is calculated as

~'(s') = B TV(I -- s 'H q)-- 1V TRU(s)

(21 )

Clearly, the poles of the reduced order system are the reciprocal of the eigenvalues of Hq. A complete pole/residue decomposition can be obtained by the eigendecomposition of Hq, Hq = TAqT -~ and subsequent substitution in (21)

~'(S') = BTVT(I - s A'q ) -1 T -1 V TRU(s)

(22)

These last two expressions are the ones used for fast frequency sweep. The inversion of ( I - s'Aq) is trivial since Aq is a diagonal matrix. From the above development it is clear that the Pad6 approximation is generated directly, without the need for calculating and matching the moments. Hence, the Arnoldi algorithm does not suffer from the numerical ill-conditioning that hinders the AWE process. As shown in [20], a qth-order Pad6 approximation of the transfer function of the reduced system generated by the Block Arnoldi process has its first r~ = Lq [pJ moments equal to those of the original transfer function. From a computational point of view, the major cost of this algorithm is associated with the LU decomposition of the matrix A. The construction of the Krylov subspace is then effected by q forward-backward substitutions. Despite its sparsity, the matrix A can be very large, considering that boundary value problems of practical interest could involve several hundreds of thousands, or even several millions of unknowns. Clearly, this is an issue that needs to be considered carefully when such model order reduction (or fast frequency sweep) algorithms are evaluated from a computational efficiency point of view, and are compared to other algorithms

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

350

that provide for broadband frequency response calculation (e.g. direct time-domain integration algorithms such as FDTD).

4. The Block Lanczos algorithm An alternative to the Block Arnoldi algorithm is the Block Lanczos algorithm. As originally proposed by Feldmann and Freund in the context of the Matrix Pad6 Via Lanczos (MPVL) process for model order reduction [22], the Block Lanczos algorithm reduces the system matrix A in (14) to a q × q block tridiagonal matrix, Tq, such that

W TA V = D T q ,

wrV=D

(23)

where the M × q matrix V forms a basis for the Krylov subspace Kr(A, R, Lq/pJ), and the M x q matrix W forms a matrix for the Krylov subspace Kr(A T, B, Iq/pJ). The two matrices V and W are generated to be bi-orthogonal. Hence, the matrix D is a diagonal matrix. The details of the algorithm have been given in [22] and will not be repeated here. The algorithm matches the first L2q/pJ moments of the resulting reduced-order transfer matrix to those of the original transfer matrix. In addition to the dominating LU decomposition of the matrix A, the block Lanczos algorithm requires 2q forward-backward substitutions for the generation of the matrices V and W to effect a Pad6 approximation of order q. In a manner similar to that for the block Arnoldi algorithm, the model order reduction effected by the Lanczos algorithm may be illustrated by introducing the change of variables

X=VX

(24)

where X is the reduced state vector of length q. Substituting (24) into (13), multiplying on the left by W v, and making use of (23) we obtain (I -- s ' T q ) f ~ : D - J W T R U ( s ')

(25)

Finally, the response of the reduced system is calculated as

Y(s') = B TV(I --

s ' T q ) - 1D -

IWTRU(s' )

(26)

This last equation is the one used for fast frequency sweep. Since q is rather small and Tq is block tri-diagonal, the inversion of the matrix (I - s'Tq) is very fast. Furthermore, the arrays V and W that seem to be needed for the calculation of the fast frequency sweep do not need to be generated and stored. Clearly, storing these arrays would be rather costly from a computer memory resources point of view, since their dimension M can be of the order of several hundreds of thousands for realistic problems. However, as explained in [22], during the Lanczos process, the arrays V and W are constructed in such a manner that it is

.-- [ao' ] where the matrices a I and a 2 are upper triangular matrices. Substituting these last expressions in (26), and using the bi-orthogonality of V and W, one obtains

y(s') = [aXzOlD(l - s'T q)-lD-'[aO1]U(s ')

(28)

Clearly, all arrays present in this final form of the response matrix are of dimension q × q. These are the only arrays needed to be stored for fast frequency sweep purposes.

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

351

5. Passivity of the discrete model The expressions given in (21) or (22) and (28) can be used directly for the purposes of fast frequency sweep. However, appropriate post-processing is required before such expressions can be used reliably as reduced-order model representations of the electromagnetic systems under study. As already mentioned in the Introduction, reduced-order model representations of passive electromagnetic multiports are extremely useful when these multiports constitute parts of more complex functional blocks. For the purposes of design-driven simulation at the functional block level, a network-oriented simulation approach is used. The incorporation of the reduced-order models for the electromagnetic multiports in the overall network-oriented circuit simulator may be effected either through recursive convolution, utilizing the poleresidue representations of the elements of the transfer function matrix [23,24], or through the direct incorporation of the state-space representation of the reduced system in the circuit simulator [22]. In either case, the passivity of the reduced system needs to be verified in order to avoid non-physical instabilities in the subsequent simulation of the overall circuit. It is important to observe that it is meaningless to talk about passivity of the reduced-order model without establishing first the passivity of the original discrete model. Therefore, our discussion of model passivity will begin with the examination of the passivity of the system of (11) which results from the specific discretization of the system of Maxwell's equations. Our analysis will make use of the following useful results [25]: THEOREM 1. The transfer function matrix l?l(s) o f a passive (linear, solveable, time-invariant) network is positive-real; that is (a) Each element of t?l(s) is analytic f o r Re(s) > 0. (b) ffI(s*) = ffl*(s) f o r Re(s) > 0. (c) flh(S) = H *T(s) + H(s) i> 0 f o r Re(s) > 0. THEOREM 2. If a matrix, A, is positive-real, then so is its inverse, A - 1, if it exists. THEOREM 3. l f A is positive-real and A h = A *v + A > O f o r Re(s) > O, then A 1 exists. THEOREM 4. If 10 is a real constant m × n matrix and A(s) is an m × m positive-real matrix, then 10 VA10 is an n × n positive-real matrix. With the output defined as Y ( s ) = B T X ( s ) , the transfer function of (11) is given by H(s) = - B T ( G + s C ) - t B

(29)

According to Theorem 1, the discrete approximation to the system of Maxwell's equations will be passive if H(s) is positive-real. However, the matrix B in (29) is a real constant matrix. Thus, in view of Theorems 2, 3 and 4, to prove the passivity of the discrete system of (11) it suffices to show that the matrix S = G + sC is positive-real and S h > O. This takes us back to Theorem 1 in which the three requirements for a matrix to be positive-real are given. First, we note that matrices G and C are real. Hence, requirements (a) and (b) are automatically satisfied. To prove that requirement (c) is also satisfied, strengthened by the additional requirement that S h > 0, we need to show that z*T(S *T + S)Z > 0 for Re(s) > 0 and for any complex vector z. Setting s = a +rio, one obtains after some straightforward matrix algebra Z*T(S *T + S)z = z*T[(G + G T) + a(C + cT)]z

(30)

where use has been made of the fact that C is symmetric. Finally, using the fact that C + C T = 2C, the above equation may be cast in the form z*T(S *T + S)z = z*T[(G + G T) + 2aC]z

(31)

Since C is positive definite, it follows immediately that 2az*TCz > 0 for a > 0. Furthermore, using the fact that G is skew-symmetric (see (12)), it is straightforward to show that

A.C, Cangellaris, L. Zhao / Comput. Methods Appl, Mech. Engrg. t69 (1999) 345-3.58

352

But D,~ is a non-negative definite matrix; hence, z*T(G + G'r)z >IO. Thus, we conclude that the product in (44) is positive definite for any complex vector z and for Re(s) > 0; hence, the discrete system of (11) is passive. It is important to point out that critical to this proof of passivity of the discrete system was the fact that the matrix G was skew-symmetric. Clearly, this skew-symmetry of G is a direct consequence of the uniformity of the (orthogonal) cartesian grid used for the discretization, as well as the staggering of the electric and magnetic field nodes. From a numerical integration point of view, it is extremely useful to be able to validate that the (semi-discrete) system of state equations resulting from the numerical approximations of the spatial derivatives is passive, since passivity guarantees stability and thus a stable numerical solution will always be achieved with a stable integration algorithm. Proof of passivity of the semi-discrete system when unstructured grids are used is rather cumbersome since, in most cases, it can be effected only through an eigenvalue analysis. Taking a closer look at the structure of the matrix G, it becomes apparent that the specific form of the numerical approximation of the two curl operators in Maxwell's system plays an important role on the passivity of the approximation. In order to account for the general case, let Pm and Pe denote the matrix approximations resulting from the discretization of the curl of the magnetic and electric field, respectively, over the entire computational domain. Then, the matrix G assumes the general form

0

G--- _p,,,

Pc]

(33)

D~

and thus, the requirement for passivity of the approximation centers around the properties of the matrix

Gh=G+GT= [

0

(p _ p V ) ]

(eo - e ~ ) ~

2D~

j

(34)

More specifically, considering that D~ is non-negative definite, the passivity of the numerical model depends solely on the properties of the matrix Q = Pe - PVm"TO elaborate, let us expand the term Z*T(Gh)Z, where Z is written in the form

It is then

Q lrLz~J ,l = z*TQzz + [(z*TQz2)T]* + 2z*TD°zz 2DJ = 2 Re{z*~Qz~} + 2z*~Ooz~

(36)

where use has been made of the fact that Q is real. The second term in the last equation is always non-negative. Considering that it is common to assume loss-free media for numerical wave simulation purposes, the passivity of the numerical model depends on whether the first term in (36) is non-negative for any z. Clearly, this is T ensured when Q = 0 or, equivalently, Pe = Pro. As already seen in the first part of this section, this occurs naturally when the (structured) Yee's lattice is used for the discretization. For other discretization choices, the passivity of the discrete model can be evaluated by examining the matrix P e - PVm. If this matrix is non-negative, the discrete model is certainly passive and its numerical integration (with a numerically stable integration scheme) will lead to a stable numerical solution. Before concluding this session we must point out the fact that the aforementioned discussion did not take into account the discrete approximation of any absorbing boundary conditions used for grid truncation purposes. In their discrete form such conditions lead to modifications in the matrices Pe and P,,, as well as the matrix C in (11), and thus impact passivity. However, the discussion of this topic is beyond the scope of this paper.

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. L)2grg. 169 (1999) 345-358

353

6. Development of the passive reduced-order model Both the Arnoldi and Lanczos algorithms presented above lead to numerically robust processes for the generation of stable, high-order Pad6 approximations of the responses of discrete electromagnetic systems. However, the passivity of the resulting reduced-order models is not guaranteed; thus, it has to be checked before one can proceed with the incorporation of the reduced-order model in a network-oriented circuit simulator. However, there exists an alternative approach to the generation of the reduction of the original discrete system that leads directly to a passive, reduced-order model. It is based on the use of congruence transformations to effect the direct reduction of the matrices G and C in (11), with the matrix V, generated by the Block Arnoldi process, being the transformation matrix. This approach was first proposed and demonstrated in [26] for the case of model order reduction of RLC circuits. Following this approach, (11) is multiplied on the left by V v and using (19) yields

[(VTGV) + s(VVCV)]f~ = -(VTB)U(s)

(37)

while the output vector is also written in terms of the state vector/~"

f = (BTV)X

(38)

The above equations may be cast in the compact form

(~ + s ~ ) k = -BU(s)

(39)

],9 = B T2

(40)

where the reduced order matrices C,, C and ]~ are given by the transformations

= V TGV = vTcv =

(4l)

VTB

The passivity of the reduced system can be proven in a manner similar to that used in the previous section to prove the passivity of the original discrete model. More specifically, to prove passivity it suffices to show that the transfer function, l~(s), of the reduced system (39), (40) given by W(s) =/~T(G + s C ) - ' / ~

(42)

is positive real. Or, in view of the theorems for positive-real matrices given in the previous section, and the fact that the matrix/~ in (42) is a real constant matrix, it suffices to show that the matrix P~ = G + sC is positive real and Ah > 0. This is done next. If the frequency shift, s 0, is taken real, the transformation matrix V is real. Thus, all matrices in (41) are real and requirements (a) and (b) in Theorem 1 of the previous section are automatically satisfied. To prove that requirement (c) is also satisfied, strengthened by the additional requirement that -4h > 0, we need to show that z*T(A *T + A)z > 0 for Re(s) > 0 and for any complex vector z. Setting s = a +jto, and using (41), one obtains after some straightforward matrix algebra z*T(A *x +A)z = z*TvT[(G + G T) + a(C + cT)]Vz

(43)

where use has been made of the fact that C is symmetric. Finally, with the change of variables w = Vz and using the fact that C + C T = 2C, the above equation may be cast in the form z*T(A *T +A)z = w*T[(G + G T) + 2aC]w

(44)

Since C is positive definite, it follows immediately that 2aw*TCw > 0 for a > 0. Finally, as it was shown in the previous section (see (32)), G + G T is a non-negative definite matrix. Thus, we conclude that the product in (44) is positive definite for any complex vector z and for Re(s)> 0; hence, the reduced system of (39), (40) is passive. From a computational complexity point of view, this algorithm is essentially the same with the Block-Arnoldi

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. t69 (1999) 345-358

354

algorithm discussed earlier. Clearly, the congruence transformation matrix, V, needs to be generated and stored before the model-order reduction process is effected through (41). Considering that the dimension of V is M × q, with M in the order of several hundreds of thousands for realistic problems, this is an expensive proposition from a memory storage point of view. Thus, despite its success with the passive reduction of RLC circuits, the feasibility of the application of this algorithm to passive reduction of three-dimensional electromagnetic systems is questionable. The obvious exception is the case of passive model order reduction of networks containing transmission line systems [27]. In this case the distributed electromagnetic systems are one-dimensional; hence, the number of degrees of freedom in the discrete model is in the order of several hundreds rather than hundreds of thousands.

7. Numerical examples The numerical examples presented in this section concern the development of Pad6 approximations for three-dimensional electromagnetic systems through the Lanczos algorithm. Through numerous simulations, it has been observed that the Lanczos algorithm allows the development of high-order Pad6 approximations with high accuracy using just a single expansion frequency point. For the electromagnetic applications of interest, the expansion point, % is taken to be s o = 2wfm,~~, where fm,x is the maximum frequency in the bandwidth of interest. As far as stopping criteria for the Lanczos algorithm are concerned, we refer the reader to the theory set forth in [28]. The following examples involve both electromagnetic propagation and radiation cases. Consequently, the FD-TD grid used for the simulation has to be truncated properly in order to facilitate accurate solutions with negligible (spurious) numerical reflection noise. This is effected using perfectly matched layers (PML) implemented through the generalized theory for PML (GT-PML) [29]. The implementation of PML in a way that is compatible with the Lanczos algorithm has been discussed in detail in [30] and will not be repeated here. The first example is for the electromagnetic analysis of the three-stub waveguide filter shown in Fig. 1. The cross-sectional symmetry of the filter has been used to reduce the size of the problem by a factor of two. More specifically, the filter is operating at the TE ~0 mode; hence, a perfect magnetic conductor plane is placed on the plane of symmetry and, thus, only one half of the original structure needs to be modeled. The excitation current source is defined on a cross-sectional plane of the feeding waveguide, and is designed to be consistent with the

0

-10

- 6 0

-4 ~17.23q 1-17.23q I2.86

2.01

2.86

Fig. 1. A three-stub waveguide filter. All dimensions are in mm.

.

10

.

.

.

i

11

,

,

,

,

f

12

,

,

,

,

i

,

13

,

,

,

i

14

,

,

,

,

15

Frequency (GHz) Fig. 2. Magnitude of the reflection coefficient for the three-stub waveguide filter of Fig. 1. Solid line: PVL approximation of order 200 at expansion frequency f = 15 GHz; dashed-line: FEM solution from [31l; dotted line: measurement,

355

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

TElo fields. The filter is terminated by means of a six-layer PML on the source side and a four-layer PML on the output side. The PMLs are designed for reflection coefficients of 10 -5. The frequency range of interest remains below 15 GHz; consequently, only the TE~0 mode propagates. Clearly, the lower bound in the frequency range of interest is taken to be greater than the cut-off frequency of the TE~o mode. In particular, the bandwidth over which the response of the filter is calculated is 10.4-15 GHz. Fig. 2 depicts the reflection coefficient at the source, $1~, for the waveguide filter of Fig. 1. The solid line depicts the response obtained from a Pad6 approximation of order 200. The expansion point was taken to be s o = 2~r X 15 X 109 rad/s. The dashed line is obtained by solving the problem at several frequencies over this frequency band by means of the finite element method [31]. Finally, the dotted line depicts the measured response. Good agreement is observed, especially in the lower portion of the frequency band. The second example deals with the modeling of the microstrip square patch antenna of Fig. 3. The patch is placed on a conductor-backed dilectric substrate of thickness 3.175 mm and relative dielectric constant of 2.5. It is probe-fed at a point along one of the two planes of symmetry of the patch, at a distance of 14.224 cm from the edge. The source is implemented in the finite difference grid through a perfectly conducting wire that is placed at the position of the feeding-probe and connects the patch to the ground plane. One segment of this wire, of length equal to one finite-difference cell, is removed, and the electric field at this segment is used to effect a voltage source. A four-layer PML is used to truncate the computational domain. It is designed for a reflection coefficient of 10 .4 at normal incidence, and it is placed two cells above the patch and three cells away from its edges. Fig. 4 compares the calculated reflection coefficient, S~ for the antenna with measured data (dotted line) [32]. The solid line is for the numerical solution obtained from a Pad6 approximation of order 200. The expansion point for the Pad6 approximation was taken to be s o -= 2~r × 5 X l09 rad/s. The dashed line is for FD-TD results obtained through the Fourier transform of the transient response of the antenna excited by a short Gaussian pulse [32]. For the numerical solutions, the reference impedance used for the calculation of S ~ was taken to be equal to the real part of the input impedance at the first resonant point. Considering that there is some uncertainty in the actual value of the substrate dielectric constant of the patch used for the measurements (Er = 2.5-+0.05), the agreement between measured and calculated data is very good. The third example deals with the electromagnetic analysis of a narrow-band, band-pass stripline filter. The geometry of the filter is shown in Fig. 5. Because of the transverse electromagnetic field (TEM) nature of the fundamental propagating mode in the stripline, transmission line analysis is frequently used for the design of such filters. A critical parameter in such a transmission line-based design is the effective gap capacitance that 10

A

ff

-10

-20

-30

-40

i

i

I

I

2

Frequency (GHz)

Fig. 3. A probe-fed patch antenna (a =41.275 mm, h = 3.175 mm, s = 14.224 mm, ~ = 2.5).

Fig. 4. Magnitude of the reflection coefficient of the patch antenna of Fig. 3. Solid line: PVL approximation of order 200 at expansion frequency f = 5 GHz; dotted line: measurement.

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg, 169 (1999) 345-358

356

i-q Fig. 5. Geometry of a balanced strip-line filter (b = 3.2 mm, w = 2.66 ram, d = 1,515 mm, l = 9,648 mm, ~ = 2.2).

effects the coupling of the center resonant segment of the filter to the input and output lines [33]. Even small inaccuracies in the gap capacitance value may lead to unacceptable shifts in the center frequency of the pass band. Therefore, even though transmission line theory may be used for the purposes of achieving an initial design, design iteration using a full-wave solver is desirable for design optimization. Before attempting the development of a Pad6 approximation for the response of this bandpass filter, the propagation characteristics of the input and output striplines were calculated over the frequency bandwidth of 0 - 2 0 GHz. Considering that the insulating dielectric is homogeneous, it is definitely an overkill to use a full-wave approach to obtain the propagation constant (which is simply given by /3 = ~ ) and the characteristic impedance (which, for this case of a balanced stripline, can be obtained simply from closed form expressions [34]). However, the finite-difference, full-wave solver was used to develop a Pad6 approximation for the transmission coefficient for a segment of this stripline, in order to examine the accuracy of the Lanczos algorithm in generating very broadband Pad6 approximations of electromagnetic systems. Fig. 6 compares the propagation constant obtained from a Pad6 approximation of order 180 (solid line), generated from a single expansion frequency point at s o = 2at × 20 × 109 rad/s, with the analytic solution (dotted line). The agreement is excellent up to 17 GHz. Fig. 7 compares the reflection coefficient for the band-pass filter of Fig. 5, calculated using a Pad6 approximation of order 180 (solid line), generated from a single expansion frequency point at s o = 2~r × 20 × 109 rad/s, with that obtained from transmission line theory (dotted line). For the transmission line theory case the gap capacitance was calculated by approximation formulation in [35]. The observed shift in the center frequency of the pass-band is caused by the inaccuracies in the calculation of the edge coupling capacitance (or, more precisely, the equivalent circuit that needs to be used to model the edge coupling effects) in the transmission line model. As expected, this inaccuracy is worse at higher frequencies. 700

. . . .

i

. . . .

I

. . . .

i

. . . .

10

""

600

500

-10 A

II1 ~

-20

-30

C

-40

0

. ~

J i l l

5

10 Frequency

15

20

(GHz)

Fig. 6. Dispersion diagram for the stripline used in the filter of Fig. 5. Solid line: PVL approximation of order 180 at expansion frequency f = 20 GHz; dotted line; analytic solution.

i

0

=

=

=

I

5

. . . .

I

. . . .

10

I

15

,

•

•

=

I

20

Frequency(GHz) Fig. 7. Reflection coefficient for the stripline filter of Fig. 5. Solid line: PVL approximation of order 180 at expansion frequency f = 20 GHz; dotted line: transmission line theory.

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

357

8. Conclusions

Model order reduction methodologies for rapid electromagnetic simulation have been presented in this paper. The emphasis was on the use of Arnoldi and Lanczos algorithms for the generation of the Pad6 approximation of the electromagnetic responses. The numerical stability of these algorithms allows for the calculation of Pad6 approximations of very high orders, thus providing accurate, continuous frequency responses over fairly broad bandwidths. For example, for one of the numerical examples presented in this paper, a Pad6 approximation of order 180, generated using a single frequency expansion point, was used to calculate the response of a narrow-band, band-pass stripline filter over a 20 GHz bandwidth. For the purposes of this paper, the discrete electromagnetic model was generated using Yee's lattice for the coupled system of Maxwell's curl equations. Thus, the number of unknowns used in the discrete approximation is equal to that used when the same problem is solved using the FD-TD method. It is well-known that the FD-TD method can be used to generate broadband responses of electromagnetic systems through impulsive excitations. In addition, several digital signal processing methods can be used to increase the efficiency and accuracy of such an approach [1]. The advantage of calculating the electromagnetic response through the generation of a Pad6 approximation rather than impulsive FD-TD is that the Pad6 approximation leads to a continuous response over the desired frequency band without any need for additional post-processing of the numerically extracted response. On the contrary, the frequency resolution, Af, of the response calculated through a fast Fourier transformation of the FD-TD generated transient response, is restricted by the time step, At, used for the numerical integration. The broader the frequency range over which the response is sought, the shorter the exciting pulse and, thus, the smaller the time step needs to be. In addition, the presence of small geometric features in the electromagnetic structure may reduce even further the time step due to the Courant stability conditions. Since Af = 1/(N At), finer frequency resolution requires an increase in the number of time steps N. Fine frequency resolution is very important when resonant structures such as filters and antennas are analyzed for the purposes of design. Thus, it is practically impossible for FD-TD to provide for the continuous resolution obtained by the Pad6 approximation, unless, of course, one uses a method such as Prony's method to generate a pole/residue representation of the calculated FD-TD transient response. Again, this is a post-processing action, and requires a fairly long transient FD-TD simulation for good accuracy. So it appears that the development of Pad6 approximations of electromagnetic system responses via Lanczos or Arnoldi algorithms has some advantages over FD-TD-based approaches. However, it must be pointed out that the application of these algorithms requires the LU decomposition of a matrix of very large dimension. For all simulations presented here, the LU decomposition was effected through a direct method. However, for problems involving several hundreds of thousands, or even millions of unknowns, direct LU decomposition becomes prohibitive and iterative methods need be used. A brief review of the Lanczos or Arnoldi processes reveals that the iterative solution of the discrete system needs to be performed at each and every iteration of the Lanczos or Arnoldi process. Fast convergence of the iterative solver used becomes then essential for these approaches to remain competitive with FD-TD like approaches. Furthermore, the impact of such an iterative process on the accuracy of the developed Pad~ approximation is not well-understood yet. Another important issue discussed in this paper was the development of passive reduced-order models of large electromagnetic systems, for subsequent use in network-oriented simulators, was discussed. The conditions for passivity were reviewed and the important observation was made that passivity of the discrete model needs to be ensured before any attempt is made to develop a passive reduced-order model. In this context, specific constraints on the discrete forms of the curl operators in Maxwell's equations were derived which, if satisfied, will ensure the passivity of the discrete model. It was shown that these constraints are always met when Yee's discretization approach on uniform grids is used. Finally, a methodology was presented for the development of guaranteed-passive reduced-order models of discrete electromagnetic systems.

References [1] A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Norwood, MA, 1995). [2] Y.-S. Tsuei, A.C. Cangellaris and J.L. Prince, Rigorous electromagnetic modeling of chip-to-package (first-level) interconnections, IEEE Trans. Components, Hybrids, and Manufact. Tech. 16 (1993) 876-883.

358

A.C. Cangellaris, L. Zhao / Comput. Methods Appl. Mech. Engrg. 169 (1999) 345-358

[3] L.T. Pillage and R.A. Rohrer, Asymptotic waveform evaluation for timing analysis, IEEE Trans. Computer-Aided Des. CAD-9 (1990) 352-366. [4] T.K. Tang and M.S. Nakhla, Analysis of high-speed VLSI interconnects using the asymprotic waveform evaluation technique, IEEE Trans. Computer-Aided Des. CAD-1 l (1992) 341-352. [5] E. Chiprout and M.S. Nakhla, Asymptotic Waveform Evaluation And Moment Matching for Interconnect Analysis (Kluwer Academic Publishers, Dordrecht, 1994). [6] S.-Y. Kim, N. Gopal and L.T. Pillage, Time-domain macromodels for VLSI interconnect analysis, IEEE Trans. Computer-Aided Des. CAD-13 (1994) 1257-1270. [7] E. Chiprout and M.S. Nakhla, Analysis of interconnect networks using complex frequency hopping (CFH), IEEE Trans. ComputerAided Des. CAD-14 (1995) 186-200. [8] M. Celik, O. Ocali, M. Tan and A. Atalar, Pole-zero computation in microwave circuits using multipoint Pad6 approximation, IEEE Trans. Circuits Syst.-I: Fundamental Theory and Applications CAS-42 (1995) 6-13. [9] M. Celik and A.C. Cangellaris, Simulation of dispersive multiconductor transmission lines by Pad6 approximation via the Lanczos process, IEEE Trans. Microwave Theory Tech. 44 (1996) 2525-2535. [10] P. Feldmann and R.W. Freund, Efficient linear circuit analysis by Pad6 approximation via the Lanczos process, IEEE Trans. Computer-Aided Des. CAD-14 (1995) 639-649. [11] K.J. Kems, I.L. Wemple and A.T. Yang, Stable and efficient reduction of substrate model networks using congruence transforms, IEEE/ACM Proc. ICCAD (1995) 207-214. [12] L.M. Silveira, M. Kamon and J. White, Efficient reduced-order modeling of frequency-dependent coupling inductances associated with 3-D interconnect structures IEEE/ACM Proc. Design Automation Conference (1995) 376-380. [13] R.F. Remis and P. Van den Berg, A modified Lanczos reduction method for the computation of transient wavefields, Proc. IEEE Antennas and Propagation Society International Symposium, Vol. 3, Baltimore, Maryland (1996) 2088-2091. [14] J.E. Bracken and Z.J. Cendes, Transient analysis via electromagnetic fast-sweep methods and circuit models, Proc. 13th Annual Review of Progress in Applied Computational Electromagnetics, Vol. I (1997) 172-179. [15] M.A. Kolbedhari, M. Nakhla, R. Achar and M. Srinivasan, Solution of EM problems using reduced-order models by complex frequency hopping, Proc. 13th Annual Review of Progress in Applied Computational Electromagnetics, Vol. I (1997) 165-171. [16] M. Celik and A.C. Cangellaris, Simulation of multiconductor transmission lines using Krylov subspace order-reduction techniques, IEEE Trans. Computer-Aided Design 16 (1997) 485-496. [17] A.C. Cangellaris and L. Zhao, Reduced-order modeling of electromagnetic systems with Pad6 via Lanczos approximations, Proc. 13th Annual Review of Progress in Applied Computational Electromagnetics, Vol. I (1997) 148-155. [18] R.A. Rohrer and H. Nosrati, Passivity considerations in stability studies of numerical integration algorithms, IEEE Trans. on Circuits and Systems CAS-28 (1981) 857-866. [19] K.S. Yee, Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media, IEEE Trans. Antennas Propagat. 14 (1966) 302-307. [20] D.L. Boley, Krylov space methods on state-space control models, Circuits Syst. Signal Process 13 (1994) 733-758. [21] K. Gallivan, E. Grimme and P. Van Dooren, Pad6 approximation of large-scale dynamic systems with Lanczos method, Proc. 33rd Conf. on Decision and Control, Buena Vista, FL (1994) 443-448. [22] P. Feldmann and R.W. Freund, Reduced-order modeling of large linear subcircuits via a Block Lanczos algorithm, IEEE/ACM Proc. Design Automation Conference (1995) 474-479. [23] A. Semlyen and A. Dabuleanu, Fast and accurate switching transient calculations on transmission lines with ground return using recursive convolutions, IEEE Trans. Power. Appl. Syst. PAS-94 (1975) 561-571. [24] W.T. Beyene and J.E. Schutt-Aint, Accurate frequency-domain modeling and efficient circuit simulation of high-speed packaging interconnects, IEEE Trans. Microwave Theory Tech. 45 (1997) 1941-1947. [25] R.W. Newcomb, Linear Multiport Synthesis (McGraw-Hill, New York, 1966). [26] A. Odabasioglu, M. Celik and L. T. Pileggi, PRIMA: Passive reduced-order interconnect macromodeling algorithm, to appear in the 1997 Proceedings of ICCAD. [27] A.C. Cangellaris and M. Celik, Order reduction of high-speed interconnect electrical models: the issue of passivity, Proc. IEEE Symposium on IC/Package Design Integration, Santa Cruz, California, 1998. [28] R.D. Slone, A computationally efficient method for solving electromagnetic interconnect problems: The Pad6 approximation via the Lanczos process with an error bound, Master's Thesis, Department of Electrical Engineering, University of Kentucky, Lexington, KY, 1997. [29] L. Zhao and A.C. Cangellaris, GT-PML: Generalized theory of perfectly matched layers and its application to the reflectionless truncation of Finite-Difference Time-Domain grids, IEEE Trans. Microwave Theory Tech. 44 (1996) 2555-2563. [30] L. Zhao and A.C. Cangellaris, Reduced-order modeling of electromagnetic field interaction in unbounded domains truncated by perfectly matched layers, Microwave and Optical Technology Letters 17 (1998) 62-66. [31] G.C. Lizalek, J.R. Ruehl and J.R. Brauer, Multi-mode s-parameter computation using finite element and perfectly matched absorbers, Proc. 12th Annual Review of Progress in Applied Computational Electromagnetics, Vol. lI (1996) 946-953. [32] G. Aguirre, Methodologies for modeling radiated emissions from printed circuit boards and packaged electronic systems. Ph.D. Dissertation, Department of Electrical and Computer Engineering, University of Arizona, 1996. [33] P.A. Rizzi, Microwave Engineering Passive Circuits (Prentice Hall, New Jersey, 1988). [34] D.M. Pozar, Microwave Engineering (John Wiley & Sons, Inc., New York, 1998). [35] G. Matthaei, L. Young and E.M.T. Jones, Microwave Filters, Impedance-Matching Networks, and Coupling Structures (Artech House Books, Dedham, MA, 1980).