Analysis of microstrip antennas by multilevel matrix decomposition algorithm

European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2000 Barcelona, 11-14 September 2000  ECCOMAS 1 ANALYSIS OF MICROSTRIP ANTENNAS BY MULTILEVEL MATRIX DECOMPOSITION ALGORITHM Josep Parrón*, Juan M. Rius*, Alex Heldring*, Eduard Úbeda*, and Juan R. Mosig † *Department of Signal Theory and Communications (TSC) Universitat Politècnica de Catalunya (UPC) Jordi Girona 1-3, Mòdul D3, Campus Nord UPC, 08034 Barcelona, Spain e-mail: [email protected] , web page: http://www-tsc.upc.es † Laboratory of Electromagnetics and Acoustics (LEMA) École Polytechnique Fédérale de Lausanne (EPFL) LEMA, EL-B, EPFL, Ecublens, CH1015, Lausanne, Switzerland e-mail: [email protected] , web page: http://lemawww.epfl.ch Key words: Integral equations, Numerical methods, Method of Moments, Microstrip antennas, Electromagnetic scattering. Abstract. Integral equation methods (IE) are widely used in conjunction with Method of Moments (MoM) discretization for the numerical analysis of microstrip antennas. However, their application to large antenna arrays is difficult due to the fact that the computational requirements increase rapidly with the number of unknowns N. Several techniques have been proposed to reduce the computational cost of IE-MoM. The Multilevel Matrix Decomposition Algorithm (MLMDA) has been implemented in 3D for arbitrary perfectly conducting surfaces discretized in Rao, Wilton and Glisson linear triangle basis functions . This algorithm requires an operation count that is proportional to N·log2N. The performance of the algorithm is much better for planar or piece-wise planar objects than for general 3D problems, which makes the algorithm particularly well-suited for the analysis of microstrip antennas. The memory requirements are proportional to N·logN and very low. The main advantage of the MLMDA compared with other efficient techniques to solve integral equations is that it does not rely on specific mathematical properties of the Green's functions being used. Thus, we can apply the method to interesting configurations governed by special Green's functions like multilayered media. In fact, the MDA-MLMDA method can be used at the top of any existing MoM code. In this paper we present the application to the analysis of large printed antenna arrays. 2 1INTRODUCTION Integral equation methods (IE)1 are widely used in conjunction with method of moments (MoM) discretization2 for the numerical analysis of electromagnetic radiation and scattering. However, their application to large antenna arrays is difficult due to the fact that the computational requirements increase rapidly with the number of unknowns N. The full system of equations resulted when discretizing integral equations requires storage memory proportional to N 2 and an operation count in the direct solution proportional to N 3. One of the most recent techniques proposed to reduce the computational cost of IE-MoM uses a sparse representation of the LU factorization3. This approach, called LU sparse integral factored representation (LUSIFER), requires CPU time proportional to N γ , with γ ranging from 2.20 to 2.89, depending on the geometry of the problem. On the other hand, most of the efficient techniques are commonly based on an iterative solution of the system of equations with conjugate gradient (CG) or biconjugate gradient (BiCG) algorithms, where the matrix-vector products are optimized in one of these two ways. 1. Exploiting physical or mathematical properties of the matrix: the fast Fourier transform4, the fast multipole method (FMM)5, multilevel algorithms6,7, are some examples of this technique. 2. Using sets of basis and testing functions that radiate narrow beams, and thus produce impedance matrices that contain many small elements. Here, we can include impedance matrix localization (IML)8,9, and wavelet expansions10,11. A different way to obtain an approximately sparse impedance matrix is the IE-MEI12, a combined field integral equation (CFIE) formulation with different testing functions in the electric and magnetic field parts, such that most of the elements in the impedance matrix are negligible. In our opinion, the multilevel fast multipole algorithm (MLFMA)6 and the multilevel matrix decomposition algorithm (MLMDA)7 are particularly interesting due to the low computational requirements for large problems and their flexibility in the application to 3D arbitrary surfaces modeled by Rao, Wilton and Glisson (RWG) linear triangles basis functions13. Both require a number of operations per iteration for very large objects, on the order of N·log2N, which for large N, is much less than the N 2 operation count required by direct matrix-vector multiply, and memory storage on the order of N·logN. The MLFMA has been implemented very efficiently in 3D6,14. The main advantage of the MLMDA compared with the other techniques is that it does not rely on specific mathematical properties of the Green's functions being used (free space, complex exponential, Hankel...). Thus the whole procedure can be easily applied to interesting configurations governed by special Green's functions like multilayered media. In fact, the MDA-MLMDA method can be used at the top of any existing MoM code, in 2D or 3D. This algorithm has also shown to have much better performance for planar or piece-wise planar objects than for general 3D problems, so it is particularly well-suited for the analysis of large printed antennas. Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 3 2THE ALGORITHM To solve large printed antenna problems we use the electric field integral equation (EFIE) in the frequency domain. The EFIE discretized by MoM may be expressed in matrix form as1,2 – [Ei]=[Z][I] (1) where [I] are the coefficients of the induced current discretized in RWG basis functions (unknowns), [Ei] is the discretization of the incident field and [Z] is the impedance matrix. This matrix includes the Green's functions with all the information about the multilayered media as shown in15. Note that all the procedure that follows is independent of which Green's functions are being used. In the CG and the BiCG, Eq. (1) is solved iteratively for the induced current coefficients [I]. In each iteration, the main computational efforts to obtain the kth estimation of the induced current [I (k)] are the matrix-vector products [Z][I (k–1)] and [Z]H[R (k–1)], where [R (k–1)]=[Z][I (k–1)]+ [Ei](2) is the previous estimation of the residual. Using direct matrix-vector multiplication, the operation count and the memory requirements for each iteration are proportional to N 2 if N is large. 2.1 Multilevel subdivision of the object In order to reduce the huge operation count in the direct matrix-vector multiply, we divide the object into many non-overlapping subdomains. Considering a pair of subdomains, we can compute the field due to the RWG basis functions at the source subdomain tested by the RWG weighting functions at the observation subdomain as [Em]=[Zmn][In](3) where n and m, respectively, are the indexes of the RWG basis functions in the source and observation subdomains and [Zmn] is a submatrix of the impedance matrix. If all possible pairs of subdomains are considered, the matrix-vector multiplication [Z][I (k–1)] in (1) may be obtained as the addition of submatrix operations of the form (3). As can be seen in 3,7,16, when source and observation subdomains are not touching each other, the field [Em] due to the source tested at the observation domain has a much smaller number of degrees of freedom than the number of elements in [Em]. Therefore, this field can be computed in substantially fewer than NbMb operations (Nb and Mb are the number of original RWG functions in the source and observation boxes, respectively). In fact, the larger the subdomains are or the further apart they are, the larger is the saving in the computation. It is thus convenient to subdivide the object into a set of boxes, as large as possible, such that most of them are not touching one another. Figure 1 shows the subdivision of the box enclosing the object into smaller boxes at multiple levels, in the form of an octal tree. The largest boxes not touching each other are at level 2, while the smallest boxes are at level L. Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 4 Level 2 Level 3 Level 1 S Figure 1: Multilevel subdivision of the box enclosing the object Hence, we will introduce a recursive procedure which will begin at level 2 and will stop at the finest level L. There are two possible cases for each few non-empty source and observation boxes which belong to the same subdivision level l (2£ l£L): 1. Boxes are touching one another or are the same: then they are subdivided into level l+1 boxes, except if we have already reached the finest level, l=L. If this is the case, direct submatrix-vector multiplication (3) is performed, requiring the previous computation of the corresponding matrix terms zmn. 2. Boxes are not touching each other: then product (3) is computed very efficiently by the algorithms that will be introduced in the next sections. 2.2 Matrix Decomposition Algorithm (MDA) In order to compute interaction between two non-touching boxes at the same level in a efficient way, we proceed as follows. Equivalent RWG sources are defined for each box as illustrated in Figure 2. In general, they must be located in the boundary of the box (Fig. 2a), but if all of the original RWG basis functions in the box are contained in a plane (as it happens in microstrip antennas), they may be located in the boundary of the rectangle resulting from the intersection between the box and the plane (Fig. 2b). a)b) Figure 2: a) Equivalent RWG sources defined for a box, evenly distributed at vertex, edges and faces, b) if all original RWG basis functions are contained in a plane, equivalent RWG sources are only necessary at the vertices and edges of the rectangle. Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 5 Now let us use p and q to refer to the indexes of the equivalent RWG functions at the boundary of these source and observation boxes, Q will be the number of equivalent RWG functions in each box. As we said before, in the computation of (3), the floating-point operation count in the direct submatrix-vector multiplication is NbMb. Alternatively, one can find the current coefficients [Ip] at the equivalent RWG in the source box that produce the same field as the original sources [In] in the equivalent RWG at the observation box (indexes q)17. [Eq]=[Zqn ][In]=[Zqp][Ip] [Ip]=[Zqp]–1[Zqn][In](4) At the original RWG in the observation box (indexes m), the equivalent sources [Ip] produce the same field as the original sources: [Em]=[Zmp][Ip]= [Zmp] ][Zqp]–1[Zqn][In](5) The operation count is now proportional to QNb+cQ3+MbQ(6) where, in our implementation, c is on the order of 1/500. If the number of equivalent RWGs small, Q<<Nb and Q<<Mb, this alternative approach is much faster than the direct submatrixvector product (3). The minimum number of equivalent sources and size of the smallest box al level L to obtain accurate computation will be discussed later in section 2.4. 2.3 Multilevel Matrix Decomposition Algorithm (MLMDA) In the analysis of electrically large objects, for boxes in coarser levels Q can become very large, in that case another algorithm, as proposed in 7, can be used to improve the efficiency of MDA. This algorithm, that will be only outlined here, computes the field due to the original RWG in the source box, tested by the original RWG in the observation box, in three steps. 1. At the finest source box level ls=L, and the coarser observation box level lo=Lb, for all the source boxes, use Eq. (4) to compute the equivalent sources [Ip] that produce in the observation box the same field [Eq] as the original sources [In]. 2. At the next source level ls=ls21, compute [Ip] for all the source boxes that produce the same field as the original sources in the observation boxes at the next observation level lo=lo+1 for each pair of source-observation boxes. This step is repeated until the coarsest source level ls=Lb and the finest observation level lo=L are reached. 3. Finally, use (5) to compute the field [Em] due to [Ip] at the equivalent sources in the source level Lb, tested by the original RWG weighting functions contained in the boxes at the observation level L. It must be noted that a different set of equivalent sources [Ip] in the source box is used for each different observation box. Here, the number of equivalent sources Q needed is smaller than MDA and is maintained constant during the process. The operation count results to be proportional to N·log2N. Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 6 2.4 Optimum implementation and numerical considerations For each pair of non-touching boxes at the same level (far field boxes), we may estimate the operation count for the three algorithms described here: the direct submatrix-vector product is the fastest for small boxes and nearly empty boxes, the MDA is the best for medium-size boxes, while the MLMDA is the most efficient for large boxes. In our implementation, the fastest of the three algorithms is used for each pair of source and observation boxes. Source and observation boxes that are at finest level touching one another or are the same box (near field boxes), produce matrices [Zmn] that are saved in hard disk, loaded each iteration and used to compute the near field contribution by direct submatrix-vector product (3). When saved near field exceeds the available hard disk capacity, the remaining near field boxes will be recomputed in each iteration by direct multiplication as non-touching boxes at the same level. It must be noted that the most significant contribution in the computation of (1) comes from this near field. Then, as contribution of MDA and MLMDA algorithms to the overall computation is small, we can increase the permitted error of these algorithms without degrading significantly the final result. For the MDA, the minimum number of equivalent sources Q to obtain an accurate computation with Eq. (4) and (5) has been theoretically studied in 16. The results that will be presented in this paper, corresponding to the case when all of the original RWG basis functions in source and field boxes are contained in a plane (Fig. 2b), have been computed using 40 os d SS2 12Q .       =λ (7) rounded to the closest multiple of 4. Ss and So are the sizes of the source and observation boxes, respectively, d is the distance between the centers of the boxes and λ is the wavelength. Comparing this choice of Q, determined through several numerical experiments, with other presented in an early paper18, it can be seen as time performance has been improved very much with no losing of accuracy. In the case of MLMDA, besides the choice of Q (smaller than in MDA), some more improvement can be made if we realize that there is no reason to stop the MLMDA recursion at the finest level in which the object has been subdivided. In fact, as number of boxes of level L inside a box of level Lb (Lb<L) is 2(L–Lb), the number of operations (4) and (5) will grow very fast if we increase L. But the higher L the smaller error we get, we have found out that a good ratio time-error is achieved for L=Lb+2. In any case, MDA algorithm (with the new choice of Q) has shown to have better time performance than MLMDA in a so wide range of box sizes, that we haven't found yet a case big enough that MLMDA could be used improving MDA performance. Optimum size of the box at the finest level (Sb) depends a lot on the hard disk speed of our computer, if Sb is too small then little near field will be saved in hard disk and lots of far field small boxes will computed with MDA or direct multiplication, this may be slower than loading them from hard disk. When Sb is too big, much near field will be saved, including some big boxes that would be faster if we computed them by MDA instead of loading from hard disk. Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 7 To summarize, we would say that the faster is your hard disk the bigger you can pick up Sb. In our computer this value is around λ. Another important thing must be taken into account and it is numerical precission when computing impedance matrices, there are two main reasons for this, the first one is that we are using an iterative method so errors can propagate easily and give us some trouble with convergence, the other one is that matrix [Zqp] in Eq. (4) is often ill-conditioned so small errors in its computation can produce big errors when inverting it. 3RESULTS In order to validate MDA-MLMDA accuracy and performance, we have run as benchmark some of the results published in19. All results have been obtained using a Pentium II 400 MHz. computer with 384 Mb of RAM memory running a high level MATLAB 5 language. Let us first consider a 4×1 series-fed microstrip array, as depicted in Fig. 3. The number of RWG unknowns was 1552, the size of the smallest box in the multilevel subdivision was 0.5λ. The relative RMS error versus the standard MoM solution was 1.8 10-3 in the induced current, 1.4 10-3 and 1.0 10-3 in the E and H-plane radiation patterns (Fig. 4). Figure 3: Series-fed 4×1 microstrip array: L = 10.08 mm, W = 11.79 mm, d1 = 3.93 mm, d2 = 1.93 mm, L1 = 23.6 mm, L2 = 13.4 mm, L3 = 12.32 mm, thickness of substrate 1.5748 mm, εr = 2.1, f = 9.42 GHz. -80 -60 -40 -20 0 20 40 60 80 -30 -25 -20 -15 -10 -5 0 dB θ MLMDA MOM Figure 4: Series-fed 4×1 array: E plane radiation pattern Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 8 Obviously, in this first example the classic MoM is faster than our MDA-MLMDA, so to test the capability of the method let us analyze the radiation from microstrip corporate fed planar arrays. In Fig. 5 can be seen the current distributions on the 8×4 microstrip corporatefed planar array, Fig. 6 shows radiation patterns for different array configurations, the agreement with results published in19 is excellent. Finally, Table 1 summarizes computational requirements for these array configurations. There are two result columns for each configuration, in the first column, we have used a coarse discretization, this is enough for obtaining good results in radiation patterns, if we are interested in an accurate computation of input impedance or current distributions, a finer mesh should be used, second column shows the requirements for this case. 4CONCLUSIONS The MDA-MLMDA algorithm has shown excellent performance for the analysis of large microstrip arrays. The error of this method compared with MoM for small problems is negligible.The memory requirements are very low: almost all the computer memory is used for storing near field terms and the sparse incomplete LU descomposition of the preconditioning matrix for BiCG. Antenna discretizations of more than 50,000 unknowns can be solved in a few hours in a PC computer, using MATLAB 5 language, showing a much better time performance than those methods based in CG-FFT. Figure 5: The current distribution on the 8×4 microstrip corporate-fed planar array. (a) The x-directed current. (b) The y-directed current. Josep Parrón, Juan M. Rius, Alex Heldring, Eduard Úbeda, and Juan R. Mosig. 9 Figure 6: Radiation patterns for different configurations of microstrip corporate fed planar arrays.