ELECTROMAGNETIC DEVICE DESIGN SYSTEM FOR FAST FREQUENCY SWEEP BASED ON FINITE ELEMENT METHOD

Abstract

The disclosure provides an electromagnetic device design system including: one or more processors of a machine; and computer-storage medium storing instructions, which when executed by the machine, cause the machine to perform operations for EM sensitivity analysis for an electromagnetic (EM) device. The operations include: initiating physical parameters of the EM device, wherein the EM device comprises multiple ports; performing EM simulation for the EM device at a pre-solution frequency using the finite-element method (FEM); applying single-size matrix Padé via Lanczos (MPVL) method in fast frequency sweep and performing EM simulation for the EM device under excitation at each port to obtain field solution of the EM device at frequencies in a frequency range; calculating S-parameters for the multiple ports of the EM device; and calculating derivatives of the S-parameters with respect to one of the physical parameters in the frequency range.

Claims

1. An electromagnetic device design system comprising: one or more processors of a machine; and computer-storage medium storing instructions, which when executed by the machine, cause the machine to perform operations for EM sensitivity analysis for an electromagnetic (EM) device, the operations comprising: initiating physical parameters of the EM device, wherein the EM device comprises multiple ports; performing EM simulation for the EM device at a pre-solution frequency using the finite-element method (FEM); applying single-size matrix Padé via Lanczos (MPVL) method in fast frequency sweep and performing EM simulation for the EM device under excitation at each port to obtain field solutions of the EM device in a frequency range; calculating S-parameters for the multiple ports of the EM device; obtaining derivatives of a full scattering matrix for the multiple ports with respect to the physical parameters of the EM device in the frequency range; selecting a solution and updating EM device design to replace initial values of the physical parameters of the EM device with a selected solution; performing EM simulation for the EM device at a frequency of the selected solution using FEM; and determining the simulated result satisfies a physical specification of the EM device; wherein: applying single-size MPVL method in fast frequency sweep comprises: transforming a single-size system matrix with a dimension of N×N into a double-size system matrix with a dimension of 2N×2N for omitting second order terms of frequency, wherein N represents a number of elements in a field vector and the single-size system matrix is of a linear combination of global finite-element system matrices; generating a first linear system using the double-size system matrix; representing solving vectors of the first linear system with the global finite-element system matrices using the block matrix inversion method and transforming the first linear system into a second linear system, wherein the second linear system is of the single-size system matrix; and solving the second linear system by performing fast frequency sweep incorporated with MPVL method and obtaining field solution in a frequency range as the following:
x.sub.k≈x.sub.k.sup.q=∥r.sub.0∥v.sub.k.sup.q(I.sub.q−(s−s.sub.0)T.sub.k.sup.q).sup.−1e.sub.1 where s represents a frequency; s.sub.0 represents a pre-solution frequency; q is a reduced order in MPVL method; x.sub.k is a vector of field solution under excitation at port k; x.sub.k.sup.q represents a vector of q.sub.th order reduced field solution under excitation at port k; r.sub.0 represents a solution vector of the first linear system at 0.sup.th MPVL iteration; V.sub.k.sup.q is represented as V.sub.k.sup.q=[v.sub.m].sub.m=1.sup.q=[v.sub.1 v.sub.2 . . . v.sub.q], where v.sub.m is an orthonormal basis vector of the Krylov subspace for model order reduction; I.sub.q is an identity matrix with a dimension of q×q; T.sub.k.sup.q, is a reduced order matrix; and e.sub.1 is represented as e.sub.1=[1 0 . . . 0].sup.T with a dimension of 1×q.

2. An electromagnetic device design system comprising: one or more processors of a machine; and computer-storage medium storing instructions, which when executed by the machine, cause the machine to perform operations for EM sensitivity analysis for an electromagnetic (EM) device, the operations comprising: initiating physical parameters of the EM device wherein the EM device comprises multiple ports; performing EM simulation for the EM device at a pre-solution frequency using the finite-element method (FEM); solving a linear system under excitation at port j to obtain a field solution under excitation at port j by fast frequency sweep; solving an adjoint representation of the linear system under excitation at port k to obtain an adjoint field solution between port k and port j by fast frequency sweep; calculating derivatives of a S-parameter of port j and port k with respect to ϕ.sub.i based on an adjoint sensitivity formula; obtaining derivatives of a full scattering matrix for the multiple ports with respect to the physical parameters of the EM device in the frequency range; selecting a solution and updating EM device design to replace the one of the physical parameters of the EM device with the selected solution; performing EM simulation for the EM device at a frequency of the selected solution using FEM; and determining the simulated result satisfies a physical specification of the EM device; wherein the adjoint sensitivity formula is written as $\frac{\partial S_{k,j}}{\partial {\varphi }_i}=-{\hat{x}}_{k,j}^T\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j,$ where x.sub.j represents a vector of the field solution under excitation at port j; {tilde over (G)}.sub.i represents a derivative of G with respect to ϕ.sub.i; {tilde over (C)}.sub.i represents a derivative of C with respect to ϕ.sub.i; s is a frequency; {circumflex over (x)}.sub.k,j.sup.T represents a transpose vector of the adjoint field solution between port j and port k; S.sub.k,j represents the S-parameter of port j and port k; and ϕ.sub.i represents an i.sub.th one of the physical parameters.

3. An electromagnetic device design system comprising: one or more processors of a machine; and computer-storage medium storing instructions, which when executed by the machine, cause the machine to perform operations for EM sensitivity analysis for an electromagnetic (EM) device, the operations comprising: initiating physical parameters of the EM device wherein the EM device comprises multiple ports; performing EM simulation for the EM device at a pre-solution frequency using the finite-element method (FEM); solving a linear system under excitation at port j to obtain a field solution under excitation at port j by fast frequency sweep; solving the linear system under excitation at port k to obtain a field solution under excitation at port k by fast frequency sweep; calculating derivatives of a S-parameter of port j and port k with respect to ϕ.sub.i based on a self-adjoint sensitivity formula; obtaining derivatives of a full scattering matrix for the multiple ports with respect to the physical parameters of the EM device in the frequency range; selecting a solution and updating EM device design to replace initial values of the physical parameters of the EM device with the selected solution; performing EM simulation for the EM device at a frequency of the selected solution using FEM; and determining the simulated result satisfies a physical specification of the EM device; wherein the self-adjoint sensitivity formula is written as $\frac{\partial S_{k,j}}{\partial {\varphi }_i}=-{\kappa }_{k,j}{x}_k^T\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j,$ where x.sub.j represents a vector of the field solution under excitation at port j; {tilde over (G)}.sub.i represents a derivative of G with respect to ϕ.sub.i; {tilde over (C)}.sub.i represents a derivative of C with respect to ϕ.sub.i; s is a frequency; x.sub.k.sup.T represents a transpose vector of the field solution under excitation at port k; S.sub.k,j represents the S-parameter of port j and port k; ϕ.sub.i represents an i.sub.th one of the physical parameters; and K.sub.k,j represents a correlated coefficient of port j and port k.

4. An electromagnetic device design system comprising: one or more processors of a machine; and computer-storage medium storing instructions, which when executed by the machine, cause the machine to perform operations for EM sensitivity analysis for an electromagnetic (EM) device, the operations comprising: initiating physical parameters of the EM device wherein the EM device comprises multiple ports; performing EM simulation for the EM device at a pre-solution frequency using the finite-element method (FEM); solving a linear system under excitation at port j by performing fast frequency sweep with single-size MPVL method to obtain an order reduced field solution under excitation at port j in a frequency range; solving an adjoint representation of the linear system under excitation at port k by performing fast frequency sweep with single-size MPVL method to obtain an order reduced adjoint field solution between port k and port j in the frequency range; calculating derivatives of a S-parameter of port j and port k with respect to ϕ.sub.i based on an adjoint sensitivity formula; obtaining derivatives of a full scattering matrix for the multiple ports with respect to the physical parameters of the EM device in the frequency range; selecting a solution and updating EM device design to replace initial values of the physical parameters of the EM device with the selected solution; performing EM simulation for the EM device at a frequency of the selected solution using FEM; and determining the simulated result satisfies a physical specification of the EM device; wherein: performing fast frequency sweep with single-size MPVL method comprises: transforming a single-size system matrix with a dimension of N×N into a double-size system matrix with a dimension of 2N×2N for omitting second order terms of frequency, wherein N represents the number of elements in a field vector and the single-size system matrix comprises a linear combination of global finite-element system matrices; generating a first linear system using the double-size system matrix; representing solving vectors of the first linear system with the global finite-element system matrices using block matrix inversion method and transforming the first linear system into a second linear system, wherein the second linear system is of the single-size system matrix; and solving the second linear system by performing fast frequency sweep incorporated with MPVL method and obtaining a field solution in a frequency range by the following:
x.sub.k≈x.sub.k.sup.q=∥r.sub.0∥v.sub.k.sup.q(I.sub.q−(s−s.sub.0)T.sub.k.sup.q).sup.−1e.sub.1; and the adjoint sensitivity formula is written as: $\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\hat{x}}_{k,j}^{\mathrm{Tq}}\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j^q,$ where s represents a frequency; s.sub.0 represents a pre-solution frequency; q is a reduced order in MPVL method; x.sub.k is a vector of field solution under excitation at port k; x.sub.k.sup.q represents a vector of q.sub.th order reduced field solution under excitation at port k; r.sub.0 represents a solution vector of the first linear system at 0.sup.th MPVL iteration; V.sub.k.sup.q is represented as V.sub.k.sup.q=[v.sub.m].sub.m=1.sup.q=[v.sub.1 v.sub.2 . . . v.sub.q], where v.sub.m is an orthonormal basis vector of the Krylov subspace for model order reduction; I.sub.q is an identity matrix with a dimension of q×q; T.sub.k.sup.q is a reduced order matrix; e.sub.1 is represented as e.sub.1=[1 0 . . . 0].sup.T with a dimension of 1×q; {tilde over (G)}.sub.i represents a derivative of G with respect to ϕ.sub.i; {tilde over (C)}.sub.i represents a derivative of C with respect to ϕ.sub.i; {circumflex over (x)}.sub.k,j.sup.qT represents a transpose vector of qth order reduced adjoint field solution between port j and port k; S.sub.k,j represents the S-parameter of port j and port k; and ϕ.sub.i represents an i.sub.th one of the physical parameters.

5. An electromagnetic device design system comprising: one or more processors of a machine; and computer-storage medium storing instructions, which when executed by the machine, cause the machine to perform operations for EM sensitivity analysis for an electromagnetic (EM) device, the operations comprising: initiating physical parameters of the EM device wherein the EM device comprises multiple ports; performing EM simulation for the EM device at a pre-solution frequency using the finite-element method (FEM); performing fast frequency sweep with single-size MPVL method to obtain an order reduced field solution under excitation at port j of the EM device in a frequency range; performing fast frequency sweep with single-size MPVL method to obtain an order reduced field solution under excitation at port k of the EM device in a frequency range; calculating derivatives of a S-parameter of port j and port k with respect to ϕ.sub.i based on a self-adjoint sensitivity formula; obtaining derivatives of a full scattering matrix for the multiple ports with respect to the physical parameters of the EM device in the frequency range; selecting a solution and updating EM device design to replace initial values of the physical parameters of the EM device with the selected solution; performing EM simulation for the EM device at a frequency of the selected solution using FEM; and determining the simulated result satisfies a physical specification of the EM device; wherein: performing fast frequency sweep with single-size MPVL method comprises: transforming a single-size system matrix with a dimension of N×N into a double-size system matrix with a dimension of 2N×2N for omitting second order terms of frequency, wherein N represents the number of elements in a field vector and the single-size system matrix comprises a linear combination of global finite-element system matrices; generating a first linear system using the double-size system matrix; representing solving vectors of the first linear system with the global finite-element system matrices using block matrix inversion method and transforming the first linear system into a second linear system, wherein the second linear system is of the single-size system matrix; and solving the second linear system by performing fast frequency sweep incorporated with MPVL method and obtaining field solution of the EM device in a frequency range as the following:
x.sub.k≈x.sub.k.sup.q=∥r.sub.0∥V.sub.k.sup.q(I.sub.q−(s−s.sub.0)T.sub.k.sup.q).sup.−1e.sub.1, and the self-adjoint sensitivity formula is written as $\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\kappa }_{k,j}{x}_k^{\mathrm{qT}}\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j^q,$ where s represents a frequency; s.sub.0 represents a pre-solution frequency; q is a reduced order in MPVL method; x.sub.k is a vector of the field solution under excitation at port k; x.sub.k.sup.q represents q.sub.th order reduced solution vector under excitation at port k; r.sub.0 represents a solution vector of the first linear system at 0.sup.th MPVL iteration; V.sub.k.sup.q is represented as V.sub.k.sup.q=[v.sub.m].sub.m=1.sup.q=[v.sub.1 v.sub.2 . . . v.sub.q], where v.sub.m is an orthonormal basis vector of the Krylov subspace for model order reduction; I.sub.q is an identity matrix with a dimension of q×q; T.sub.k.sup.q, is a reduced order matrix; e.sub.1 is represented as e.sub.1=[1 0 . . . 0].sup.T with a dimension of 1×q; x.sub.j.sup.q represents q.sub.th order reduced solution vector under excitation at port j; {tilde over (G)}.sub.i represents a derivative of G with respect to ϕ.sub.i; {tilde over (C)}.sub.i represents a derivative of C with respect to ϕ.sub.i; x.sub.k.sup.qT represents a transpose vector of qth order reduced field solution under excitation at port k; S.sub.k,j represents the S-parameter of port j and port k; ϕ.sub.i represents an i.sub.th one of the physical parameters; and κ.sub.k,j represents a correlated coefficient of port j and port k.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0155] FIG. 1 is a flow chart for the operations performed by the invented device design system in the embodiment of this microwave filter.

[0156] FIG. 2 is a perspective view of a four-pole waveguide filter according to one embodiment of the disclosure;

[0157] FIG. 3 is another perspective view of a four-pole waveguide filter according to one embodiment of the disclosure;

[0158] FIG. 4 is graph depicting Derivatives of ReS.sub.11 and ImS.sub.11 with respect to h.sub.1 at ϕ.sub.1 for the four-pole waveguide filter with the frequency swept from 10.5 GHz-11.5 GHz;

[0159] FIG. 5 is graph depicting Derivatives of ReS.sub.21 and ImS.sub.21 with respect to h.sub.1 at ϕ.sub.1 for the four-pole waveguide filter with the frequency swept from 10.5 GHz-11.5 GHz;

[0160] FIG. 6 is graph depicting Derivatives of ReS.sub.11 and ImS.sub.11 with respect to h.sub.c1 at ϕ.sub.1 for the four-pole waveguide filter with the frequency swept from 10.5 GHz-11.5 GHz; and

[0161] FIG. 7 is graph depicting Derivatives of ReS.sub.21 and ImS.sub.21 with respect to h.sub.c2 at ϕ.sub.1 for the four-pole waveguide filter with the frequency swept from 10.5 GHz-11.5 GHz.

DETAILED DESCRIPTION

[0162] For an EM device design system, an adjoint formula and a self-adjoint formula that are respectively integrated with fast frequency sweeps are generated to improve EM sensitivity analysis for EM device designs. Further, a single-size MPVL method is proposed for model order reduction and is used in fast frequency sweep for the present EM sensitivity analysis. The present EM sensitivity analysis method obtains the same accuracy as the techniques of the prior art, while taking much less time by avoiding repetitively solving large systems of EM equations for different frequencies.

[0163] The single-size MPVL method incorporated in fast frequency sweep is derived as follows:

[0164] the Helmholtz equation for the full wave EM simulation is formulated as:

[00005]$\begin{array}{cc}\nabla \times \left(\frac{1}{\mu }\nabla \times E\right)-{\omega }^2ϵE=-j\omega J,& \left(1\right)\end{array}$

[0165] where ∈ and μ represent the permittivity and the permeability of the medium, respectively; ω represents the angular frequency; J represents the electric current source of the EM problem; E represents the electric field intensity to be solved; finite element method (FEM) is one of the commonly used methods to solve this Helmholtz equation; assuming that the EM simulation is performed on an EM device comprising multiple ports; the FEM equation in a general form is written as:

(K.sub.0+sK.sub.1+s.sup.2K.sub.2)x=sb.sub.j,  (2)

[0166] where s represents the complex frequency corresponding to ω; x represents the unknown field solution vector containing the unknown values of E for FEM; K.sub.0, K.sub.1, and K.sub.2 represent global finite-element system matrices, which are dependent on the physical parameters of the EM device but independent on frequency; b.sub.j represents a vector describing the excitation at the port j; to obtain the full scattering matrix of the multi-port structure, the S-parameters can be written as

[00006]$\begin{array}{cc}\begin{array}{c}{S}_{k,j}\left(s\right)={\kappa }_{k,j}{b}_k^{T}x-{\delta }_{k,j}\\ =s{\kappa }_{k,j}{b_k^T\left(K_0+{\mathrm{sK}}_1+s^2{K}_2\right)}^{-1}{b}_j-{\delta }_{k,j},\end{array}& \left(3\right)\end{array}$$\mathrm{where}$$\begin{array}{cc}{\delta }_{k,j}=\{\begin{array}{cc}1,& k=j\\ 0,& k\ne j\end{array},& \left(4\right)\end{array}$

[0167] κ.sub.k,j is a correlated coefficient of port j and port k, which is constant for a constant EM device and is dependent on the powers incident upon the port k and port j of the EM device.

[0168] To use MPVL method on FEM system, the order of frequency terms contained in the EM equations would be no more than a first order. Therefore, Equation (3) is reformatted as:

[00007]$\begin{array}{cc}{S}_{k,j}\left(s\right)={{\hat{\beta }}_{k,j}^T\left(G+\mathrm{sC}\right)}^{-1}{\beta }_j,& \left(5\right)\end{array}$$\mathrm{where}$$\begin{array}{cc}{\hat{\beta }}_{k,j}=\left[\begin{array}{c}s{\kappa }_{k,j}{b}_k\\ 0_{N\times 1}\end{array}\right],G=\left[\begin{array}{cc}{0}_N& I_N\\ K_0& K_1\end{array}\right],C=\left[\begin{array}{cc}-I_N& 0_N\\ 0_N& K_2\end{array}\right],{\beta }_j=\left[\begin{array}{c}{0}_{N\times 1}\\ b_j\end{array}\right].& \left(6\right)\end{array}$

[0169] where I.sub.N and 0.sub.N represent the identity matrix and zero matrix in .sup.N×N, respectively. Here, N is the number of elements in x. Let ϕ represent a vector of the physical parameters of the EM device, where ϕ.sub.i (i=1, . . . , p) represent the i.sub.th element in ϕ, and p is the total number of physical parameters.

[0170] To use MPVL method to perform the fast frequency sweep, one specific frequency is selected and defined as the pre-solution frequency; with the suitable selection of the pre-solution frequency, the accurate EM solutions can be obtained over the entire frequency band by solving the large linear system once at the pre-solution frequency; let s.sub.0 represent the pre-solution frequency for MPVL method; let K.sub.s represent a single-size system matrix at the pre-solution frequency, which has a dimension of N×N and is written as

K.sub.s=K.sub.0+s.sub.0K.sub.1+s.sub.0.sup.2K.sub.2,  (7)

[0171] let A.sub.s represent a double-size system matrix having a dimension of 2N×2N and is defined as

[00008]$\begin{array}{cc}{A}_s=G+s_{0}C=\left[\begin{array}{cc}-s_0{I}_N& I_N\\ K_0& K_1+s_0{K}_2\end{array}\right].& \left(8\right)\end{array}$

[0172] where A.sub.s is calculated at the pre-solution frequency s.sub.0; let q be the order of the reduced model using MPVL; let T.sub.k.sup.q be defined as the reduced order matrix with the elements t.sub.l,m, i.e.,

[00009]$T_k^q={\left[t_{l,m}\right]}_{q\times q},$

where l=1, 2, . . . , q and m=1, 2, . . . , q; let {w.sub.m}.sub.m=1.sup.q be defined as a set of unit vectors which are orthonormal to each other; let V: be defined as the matrix containing the orthonormal basis of the Krylov subspaces for the model order reduction, where V.sub.k.sup.q=[v.sub.m].sub.m=1.sup.q=[v.sub.1 v.sub.2 . . . v.sub.q]; let r.sub.m be defined as the vector for the calculation of V.sub.k.sup.q and T.sub.k.sup.q in the m.sub.th MPVL iteration; r.sub.m is calculated by solving a first linear system as:

[00010]$\begin{array}{cc}{A}_s{r}_m=\{\begin{array}{cc}{\beta }_k,& m=0,\\ -{\mathrm{Cv}}_m,& m=1,2,.Math.,q.\end{array}& \left(9\right)\end{array}$

[0173] using MPVL as the model order reduction technique to perform the fast frequency sweep; calculating the LU factors of the matrix A.sub.s;

[0174] because the size of A.sub.s is twice the size of K.sub.s, the computational cost for the LU factorization of A.sub.s takes much longer time than that of Ks;

[0175] therefore, a single-size MPVL method is derived to perform the LU factorization of the single-size system matrix K.sub.s, instead of the double-size system matrix A.sub.s; using the block matrix inversion method, Equation (9) is reformulated as:

[00011]$\begin{array}{cc}\begin{array}{c}{r}_0=\frac{1}{s_0}\left[\begin{array}{cc}{K}_s^{-1}{K}_0-I& s_0{K}_s^{-1}\\ s_0{K}_s^{-1}{K}_0& s_0^2{K}_s^{-1}\end{array}\right]\left[\begin{array}{c}{0}_{N\times 1}\\ b_k\end{array}\right]\\ =\left[\begin{array}{c}{I}_N\\ s_0{I}_N\end{array}\right]K_s^{-1}{b}_k,\end{array}& \left(10\right)\end{array}$$\mathrm{and}$$\begin{array}{cc}\begin{array}{c}{r}_m=\frac{1}{s_0}\left[\begin{array}{cc}{K}_s^{-1}{K}_0-I& s_0{K}_s^{-1}\\ s_0{K}_s^{-1}{K}_0& s_0^2{K}_s^{-1}\end{array}\right]\left[\begin{array}{cc}{I}_N& 0_N\\ 0_N& -K_2\end{array}\right]v_m\\ =\frac{1}{s_0}\left[\begin{array}{c}{I}_N\\ s_0{I}_N\end{array}\right]{K_s^{-1}\left[\begin{array}{c}{K}_0\\ -s_0{K}_2\end{array}\right]}^T{v}_m\\ -\frac{1}{s_0}\left[\begin{array}{cc}{I}_N& 0_N\\ 0_N& 0_N\end{array}\right]v_{m,}\end{array}& \left(11\right)\end{array}$

[0176] let u.sub.m be defined as the solution vector for a second linear system as

[00012]$\begin{array}{cc}{K}_s{u}_m=\{\begin{array}{cc}{b}_k,& m=0\\ {\frac{1}{s_0}\left[\begin{array}{c}{K}_0\\ -s_0{K}_2\end{array}\right]}^T{v}_m,& m=1,2,.Math.,q.\end{array}& \left(12\right)\end{array}$

[0177] substituting (12) into (10) and (11),

[00013]$\begin{array}{cc}{r}_0=\left[\begin{array}{c}{u}_0\\ s_0{u}_0\end{array}\right],& \left(13\right)\end{array}$$\mathrm{and}$$\begin{array}{cc}{r}_m=\left[\begin{array}{c}{u}_m\\ s_0{u}_m\end{array}\right]-\frac{1}{s_0}\left[\begin{array}{cc}{I}_N& 0_N\\ 0_N& 0_N\end{array}\right]v_m,& \left(14\right)\end{array}$

[0178] where m represents the MPVL iteration index, i.e., m=1, 2, . . . , q; let x.sub.k.sup.q represent the q.sub.th order reduced vector (q×1 vector), the solution vector x.sub.k under excitation at port k of the device can be formulated as,

x.sub.k≈x.sub.k.sup.q=∥r.sub.0∥V.sub.k.sup.q(I.sub.q−(s−s.sub.0)T.sub.k.sup.q).sup.−1e.sub.1,  (15)

[0179] where e.sub.1[1 0 . . . 0].sup.T∈.sup.q, is the unit vector in .sup.q; I.sub.q is an identity matrix in .sup.q×q.

[0180] according to Equation (3), the S-parameters can be obtained by substituting the solution vector x.sub.k with the Equation (15) and be formulated as

S.sub.k,j(S)=∥r.sub.0∥{circumflex over (β)}.sub.k,j.sup.TV.sub.k.sup.q(I.sub.q−(s−s.sub.0)T.sub.k.sup.q).sup.−1e.sub.1−δ.sub.k,j,  (16)

[0181] where, S-parameters corresponding to a frequency range can be evaluated by performing fast frequency sweep integrated with the single-size MPVL method.

[0182] The following is derivation of the adjoint EM sensitivity analysis using fast frequency sweep model order reduction technique.

[0183] Letting {tilde over (G)}.sub.i and {tilde over (C)}.sub.i represent the derivatives of G and C with respect to a physical parameter of the device ϕ.sub.i, respectively. {tilde over (G)}.sub.i and {tilde over (C)}.sub.i can then be derived as,

[00014]$\begin{array}{cc}{\tilde{G}}_i=\left[\begin{array}{cc}{0}_N& 0_N\\ {\tilde{K}}_0^i& {\tilde{K}}_1^i\end{array}\right],{\tilde{C}}_i=\left[\begin{array}{cc}{0}_N& 0_N\\ 0_N& {\tilde{K}}_2^i\end{array}\right].& \left(17\right)\end{array}$

[0184] where {tilde over (K)}.sub.0.sup.i, {tilde over (K)}.sub.1.sup.i, and {tilde over (K)}.sub.2.sup.i represent the sensitivity of K.sub.0, K.sub.1, and K.sub.2 with respect to the physical parameter ϕ.sub.i, respectively; to calculate the sensitivity of S.sub.i, j, the sensitivities of the global finite-element system matrices w.r.t. physical parameters {tilde over (K)}.sub.0.sup.i, {tilde over (K)}.sub.1.sup.i, and {tilde over (K)}.sub.2.sup.i need to be firstly obtained;

[0185] obtaining the vector b.sub.j from the EM excitation and/or the boundary conditions on the j.sub.th port; the structures of the ports normally do not change during the EM design; thus b.sub.j is a constant vector w.r.t. to a physical parameter ϕ.sub.i in the EM problem, i.e.,

[00015]$\begin{array}{cc}\frac{\partial b_j}{\partial {\varphi }_i}=0_{N\times 1}.& \left(18\right)\end{array}$

[0186] based on (5), formulating the derivatives of S.sub.k,j as,

[00016]$\begin{array}{cc}\frac{\partial S_{k,j}}{\partial {\varphi }_i}=-{{\hat{\beta }}_{k,j}^T\left(G+\mathrm{sC}\right)}^{-1}\left({\tilde{G}}_i+s{\tilde{C}}_i\right){\left(G+\mathrm{sC}\right)}^{-1}{\beta }_j.& \left(19\right)\end{array}$

[0187] letting x.sub.j represent the solution vector for a linear system that is written as

(G+sC)x.sub.j=β.sub.j.  (20)

[0188] calculating the solution vector x.sub.j in (20) by MORe techniques such as MPVL method;

[0189] defining {circumflex over (x)}.sub.k,j as a vector of adjoint field solution of the adjoint representation of the linear system (20), the adjoint representation is written as follows,

(G+sC).sup.T{circumflex over (x)}.sub.k,j={circumflex over (β)}.sub.k,j.  (21)

[0190] the adjoint solution vector {circumflex over (x)}.sub.k,j in (21) is calculated similarly to x.sub.j by MORe techniques such as MPVL method;

[0191] substituting (20) and (21) into (19) to obtain the adjoint sensitivity formula of S.sub.k,j w.r.t. ϕ.sub.i as,

[00017]$\begin{array}{cc}\frac{\partial S_{k,j}}{\partial {\varphi }_i}={\hat{x}}_{k,j}^T\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j,& \left(22\right)\end{array}$

[0192] Equation (22) is the adjoint sensitivity formula which is well fit for fast frequency sweep. MORe technique can be used to calculate x.sub.j and {circumflex over (x)}.sub.k,j; since the format of (22) is derived using MPVL, x.sub.j and {circumflex over (x)}.sub.k,j can be formulated using MPVL method,

[00018]$\begin{array}{cc}\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\hat{x}}_{k,j}^{qT}\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)x_j^q.& \left(23\right)\end{array}$

[0193] The following is derivation of the self-adjoint EM sensitivity analysis using fast frequency sweep model order reduction technique.

[0194] The disclosure further provides a self-adjoint EM sensitivity analysis. Because G and C are formulated as 2×2 block matrices, by substituting (6) into (20) and using the 2×2 block matrix inversion method, we can obtain,

[00019]$\begin{array}{cc}\begin{array}{c}{x}_k=\frac{1}{s}\left[\begin{array}{cc}{K}^{-1}{K}_0-I_N& s{K}^{-1}\\ s{K}^{-1}{K}_0& s^2{K}^{-1}\end{array}\right]\left[\begin{array}{c}{0}_{N\times 1}\\ b_k\end{array}\right]\\ =\left[\begin{array}{c}{K}^{-1}{b}_k\\ s{K}^{-1}{b}_k\end{array}\right],\end{array}\mathrm{where}& \left(24\right)\end{array}$$\begin{array}{cc}K=K_0+s{K}_1+s^2{K}_2.& \left(25\right)\end{array}$

[0195] because K.sub.0, K.sub.1, and K.sub.2 are symmetrical matrices, substituting (6) into (21) and using the 2×2 block matrix inversion method to obtain,

[00020]$\begin{array}{cc}\begin{array}{c}{\hat{x}}_{k,j}=\frac{1}{s}\left[\begin{array}{cc}{K}^{-1}{K}_0-I_N& s{K}^{-1}{K}_0\\ s{K}^{-1}& s^2{K}^{-1}\end{array}\right]\left[\begin{array}{c}s{\kappa }_{k,j}{b}_k\\ 0_{N\times 1}\end{array}\right]\\ ={\kappa }_{k,j}\left[\begin{array}{c}{K}^{-1}{K}_0{b}_k-b_k\\ s{K}^{-1}{b}_k\end{array}\right].\end{array}& \left(26\right)\end{array}$

[0196] defining ξ.sub.k as the vector of weighted differences between {circumflex over (x)}.sub.k,j and x.sub.k, formulated as,

[00021]$\begin{array}{cc}{\xi }_k=\frac{{\hat{x}}_{k,j}}{{\kappa }_{k,j}}-x_k=\left[\begin{array}{c}{K}^{-1}\left(K_0{b}_k-b_k\right)-b_k\\ 0_{N\times 1}\end{array}\right].& \left(27\right)\end{array}$

[0197] formulating {circumflex over (x)}.sub.k,j by the weighted summation of x.sub.k and ξ.sub.k as,

{circumflex over (x)}.sub.k,j=κ.sub.k,jx.sub.k+κ.sub.k,jξ.sub.k,  (28)

[0198] substituting (17) into (22) to obtain,

[00022]$\begin{array}{cc}\begin{array}{c}\frac{\partial S_{k,j}}{\partial {\varphi }_i}=-{\kappa }_{k,j}\left(x_k^T+{\xi }_k^T\right)\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)x_j\\ =-{\kappa }_{k,j}\left[x_k^T\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)+{\xi }_k^T\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)\right]x_j.\end{array}& \left(29\right)\end{array}$

[0199] because

[00023]$\begin{array}{cc}\begin{array}{c}{\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)}^T{\xi }_k=\left[\begin{array}{cc}{0}_N& {\stackrel{\sim }{K}}_0^i\\ 0_N& {\stackrel{\sim }{K}}_1^i+s{\stackrel{\sim }{K}}_2^i\end{array}\right]\left[\begin{array}{c}{K}^{-1}\left(K_0{b}_k-b_k\right)-b_k\\ 0_{N\times 1}\end{array}\right]\\ =0_{2N\times 1},\end{array}& \left(30\right)\end{array}$

[0200] simplifying the formulation for calculating the derivative (29) as,

[00024]$\begin{array}{cc}\frac{\partial S_{k,j}}{\partial {\varphi }_i}=-{\kappa }_{k,j}{x}_k^T\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)x_j.& \left(31\right)\end{array}$$\begin{array}{cc}\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\kappa }_{k,j}{x}_k^{qT}\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)x_j^q.& \left(32\right)\end{array}$

[0201] Equation (20) is the self-adjoint sensitivity formula which is well fit for fast frequency sweep; MORe techniques can be used to calculate x.sub.j; since the format of (31) is derived using MPVL, x.sub.j can be formulated using MPVL method.

[0202] By using Equations (12)-(14), the disclosure performs the LU factorization and forward/backward substitutions of the original system K.sub.s (single-size system matrix) instead of A.sub.s (double-size system matrix), to reduce calculation time r.sub.m, m=0, 1, . . . , q. Because the MPVL method takes a long time to calculate r.sub.m, it is simplified for efficient x.sub.k calculation. Note that, according to Equations (12)-(14), the disclosure uses only one LU factorization of the original system K.sub.s, whereas the MPVL method uses q+1 forward/backward substitutions.

[0203] Because Equation (31) is derived using MPVL, x.sub.j and {circumflex over (x)}.sub.k,j can also be formulated using the MPVL method,

[00025]$\begin{array}{cc}\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\kappa }_{k,j}{x}_k^{qT}\left({\stackrel{\sim }{G}}_i+s{\tilde{C}}_i\right)x_j^q,& \left(32\right)\end{array}$

[0204] where x.sub.j.sup.q is calculated using fast frequency sweep with MORe technique as in equation (15). Note that, in the single-size MPVL method, {w.sub.m}.sub.m=1.sup.q makes each λ.sub.m equal to non-zero. When performing MPVL method, Equation (32) calculates the derivatives of S-parameters for the self-adjoint EM sensitivity analysis. Fast frequency sweep uses a reduced-order model of order q to estimate the S-parameters of a p port structure and the derivatives thereof, and involves only one LU factorization and p (q+1) forward/backward substitutions. The computational time of the self-adjoint EM sensitivity analysis is not related to the number of physical parameters and increases slightly with the number of frequencies to be solved.

[0205] The comparison of the different methods of EM sensitivity analysis of a four-pole waveguide filter is shown in Table 1. A flow chart for the operations performed by the invented device design system in the embodiment of this microwave filter is shown in FIG. 1. The detailed number of LU factorizations and FB substitutions of the original FEM matrix for the proposed adjoint EM sensitivity analysis using fast frequency sweep is illustrated in Table 1. The disclosure further provides details of the self-adjoint EM sensitivity analysis using discrete frequency sweep, the sensitivity analysis using the finite difference method (i.e., perturbing variables, re-simulating EM and computing the solution difference) for both fast frequency sweep and discrete frequency sweep. For illustration purpose, the disclosure further provides the details of EM simulation using both fast frequency sweep and discrete frequency sweep. In Table 1, n represents the number of physical parameters; n.sub.f represents the number of frequencies; q represents the order for MPVL method; p represents the number of ports of the microwave filter. From Table 1, if the EM application only simulates at 1 or 2 frequencies, the sensitivity analysis using self-adjoint EM sensitivity analysis with discrete frequency sweep is performed. If the EM application only has 1 or 2 physical parameters, the sensitivity analysis using finite difference method with fast frequency sweep is performed. If the EM application has more physical parameters and needs to simulate at more frequencies, the adjoint EM sensitivity analysis using fast frequency sweep obtains the results faster than the existing sensitivity analysis methods. The larger number of physical parameters and frequencies the EM application has, the more significant advantage the proposed method will have.

TABLE-US-00001 TABLE 1 Comparison of the Number of LU Factorizations and F/B Substitutions of FEM Matrix for Different Methods No. of LU No. of F/B Method Factorizations Substitutions Proposed Self-Adjoint 1 p(q + 1) EM Sensitivity Analysis Using Fast Frequency Sweep Proposed Adjoint EM 1 p(2q + 2) Sensitivity Analysis Using Fast Frequency Sweep Self-Adjoint EM n.sub.f p .Math. n.sub.f Sensitivity Analysis Using Discrete Frequency Sweep Finite-Difference n + 1 p(q + l)(n + 1) Sensitivity Using Fast Frequency Sweep Finite-Difference n.sub.f(n + 1) p .Math. n.sub.f(n + 1) Sensitivity Analysis Using Discrete Frequency Sweep Fast Frequency Sweep 1 p(q + 1) Simulation Discrete Frequency n.sub.f p • n.sub.f Sweep Simulation

[0206] n—Number of Physical Parameters;

[0207] n.sub.f- Number of Frequencies at a Frequency range;

[0208] q—Order for MPVL method;

[0209] p—Number of ports of the microwave filter.

[0210] Remarks: Note that, the self-adjoint EM sensitivity method takes the same number of LU decomposition and F/B substitutions as EM simulation using fast frequency sweep using finite element method. In another word, after performing EM simulation using fast frequency sweep, no extra LU or F/B time is needed to calculate the EM sensitivities wr.t. all variables. Only a small amount of extra time of matrix multiplications is needed for EM sensitivity analysis after EM simulation.

[0211] 4. Sensitivity Analysis for Four-Pole Waveguide Filter

[0212] The example under consideration is the EM sensitivity analysis of a four-pole waveguide filter. The tuning elements are penetrating posts of square cross section placed at the center of each cavity and each coupling window, shown in FIG. 2. Let h.sub.1, h.sub.2, h.sub.3 represent the heights of the tuning posts in the coupling windows. Let h.sub.c1 and h.sub.c2 represent the heights of the square cross section placed in the center of the resonator cavities. Let w.sub.c and w.sub.p represent the widths of the square cross section placed at the center of the resonant cavity and the coupling window respectively. Let d.sub.12 and d.sub.23 represent the diameters of resonant cavities. The iris widths w.sub.1, w.sub.2, w.sub.3 and w.sub.4 are pre-tuned to achieve the desired coupling bandwidth for each section of the filter. The input and output waveguides, as well as the resonant cavities, are standard WR-75 waveguides (a=19.05 mm and b=9.525 mm). For this example, the design space vector is ϕ=[d.sub.12 d.sub.23 h.sub.1 h.sub.2 h.sub.3 h.sub.c1 h.sub.c2 w.sub.1 w.sub.2 w.sub.3 W.sub.p W.sub.c].sup.T The sensitivity analysis is performed at the nominal design point, i.e., (P.sub.i=[10.5 13.3 3.39 4.15 3.62 3.30 2.98 8.7 5.1 5.1 2.0 4.0].sup.T mm. The frequency range for this filter is selected as 10.5 GHz-11.5 GHz with the step of 20 MHz.

[0213] The derivatives of S-parameters for the microwave filter device is calculated by the proposed self-adjoint EM sensitivity analysis using fast frequency sweep for the EM simulation of this filter example. The order for the reduced order model is 16, i.e., q=16. Through the instant EM device design system disclosed in this embodiment for the microwave filter, physical specifications, such as passband bandwidth, return loss, and passband ripple can be satisfied.

[0214] For comparison purposes, the disclosure further provides the self-adjoint EM sensitivity analysis using discrete frequency analysis, the sensitivity analysis using the finite difference method for fast frequency sweep and discrete frequency analysis. Three cases for the sensitivity analysis with different number of frequency points are performed for the comparison: Case 1 with 11 frequency points and Case 2 with 51 frequency points. The comparison of the different methods of sensitivity analysis for the four-pole waveguide filter example is shown in Table 2.

TABLE-US-00002 TABLE 2 Comparison of CPU of different sensitivity analysis methods for the four- pole waveguide filter No. of LU No. of F/B Time of LU Total Sensitivity method factorizations substitutions Factorizations Time 11 Frequencies Proposed Self- 1 34 2.1 min 4.1 adjoint EM min sensitivity analysis Using Fast frequency sweep Proposed Adjoint 1 68 2.1 min 6.2 EM Sensitivity min Analysis Using Fast frequency sweep Self-adjoint EM 11 22 22.3 min 23.9 sensitivity analysis min Using Discrete Frequency Analysis Finite-Difference 13 442 26.4 min 41.2 Sensitivity Analysis min Using Fast frequency sweep Finite-Difference 143 286 4.7 h 4.9 h Sensitivity Analysis Using Discrete Frequency Sweep 101 Frequencies Proposed Self- 1 34 2.1 min 4.1 adjoint EM min sensitivity analysis Using Fast Frequency Sweep Proposed Adjoint 1 68 2.1 min 6.2 Sensitivity Analysis min Using Fast frequency sweep Self-adjoint EM 102 204 3.3 h 3.4 h sensitivity analysis Using Discrete Frequency Sweep Finite-Difference 13 442 26.4 min 41.2 Sensitivity Analysis min Using Fast Frequency Sweep Finite-Difference 1326 2652 44.3 h 45.1 Sensitivity Analysis h Using Discrete Frequency Sweep

[0215] Because 12 physical parameters are used in this example, 13 complete evaluations of the S-parameters (once at nominal values and 12 perturbations for 12 different physical parameters) are needed for the sensitivity analysis using the finite different method. For the self-adjoint EM sensitivity analysis using discrete frequency analysis, the cost increases as the number of frequency points for the EM design increases. From Table 2, the self-adjoint EM sensitivity analysis uses only one LU factorization and 34 forward/backward substitutions to obtain the derivatives of the S-parameters w.r.t 12 different physical parameters. The LU factorizations and forward/backward substitutions are the most time-consuming part during the overall sensitivity analysis process. The number of LU factorizations and forward/backward substitutions used by the proposed method is much fewer than that used by the existing methods listed in Table 2, therefore, the disclosed method takes much less time than the existing methods.

[0216] FIG. 4 and FIG. 5 show detailed comparisons of the derivatives of S.sub.11 and S.sub.21 with respect to h.sub.1 for the four-pole waveguide filter with the frequency swept from 10.5 GHz-11.5 GHz, respectively. As shown in FIG. 3 and FIG. 4, the derivatives calculated using the self-adjoint EM sensitivity analysis match very well with the derivatives calculated using the existing self-adjoint method and finite difference methods. To further verify the accuracy of the proposed self-adjoint EM sensitivity analysis method, the disclosure further provides the calculation of the derivatives of S-parameters w.r.t. other physical parameters. Detailed comparisons of the derivatives of S.sub.11 and S.sub.21 with respect to h.sub.c2 for this four-pole waveguide filter are shown in FIG. 6 and FIG. 7, respectively. From the comparisons, the self-adjoint EM sensitivity analysis using fast frequency sweep achieves the correct solution with a very low computational cost.

[0217] Appendix 1: Single-Size MPVL Method for EM sensitivity Analysis

[0218] 1) Define δ.sub.k,j, κ.sub.k,j, e.sub.1, I.sub.q, K.sub.0, K.sub.1, K.sub.2, b.sub.k, b.sub.j. Define frequency s. Define pre-solution frequency s.sub.0.

[0219] Define {w.sub.m}.sub.m=1.sup.q. Set K.sub.s=K.sub.0+s.sub.0K.sub.1+s.sub.0.sup.2K.sub.2

[0220] 2) Obtain u.sub.0 by solving the linear system K.sub.su.sub.0=b.sub.1. Set

[00026]$r_0=\left[\begin{array}{c}{u}_0\\ s_0{u}_0\end{array}\right]\mathrm{and}r=r_0.$

[0221] For m=1, 2, q do:

[0222] 3) Set

[00027]$v_m=\frac{r}{.Math.r.Math.}\mathrm{and}{\lambda }_m=w_m^T{v}_m.$

[0223] 4) If m>1, set

t.sub.m,m−1=∥r∥.

[0224] 5) Solving the linear system to obtain um,

[00028]$K_s{u}_m={\frac{1}{S_0}\left[\begin{array}{c}{K}_0\\ -s_0{K}_2\end{array}\right]}^T{v}_m.$

[0225] Calculate

[00029]$r=r_m=\left[\begin{array}{c}{u}_m\\ s_0{u}_m\end{array}\right]-\frac{1}{S_0}\left[\begin{array}{cc}{I}_N& 0_N\\ 0_N& 0_N\end{array}\right]v_m.$

[0226] 6) For l=1, 2, . . . , m do: Set

[00030]$t_{l,m}=\frac{w_l^{T}r}{{\lambda }_l}\mathrm{and}r=r-t_{l,m}{v}_l.$

[0227] End For

[0228] End For

[0229] 7) Set

T.sub.j.sup.q=[t.sub.l,m].sub.q×q,V.sub.j.sup.q=[v.sub.1v.sub.2 . . . v.sub.q].

x.sub.j≈∥r.sub.0∥V.sub.j.sup.q(I.sub.q−(s−s.sub.0)T.sub.j.sup.q).sup.−1e.sub.1.

[0230] 8) Calculate

[00031]$S_{k,j}\left(s\right)\approx s.Math.r_0.Math.{{\kappa }_{k,j}\left[\begin{array}{c}{b}_k\\ 0_{N\times 1}\end{array}\right]}^T{V_j^q\left(I_q-\left(s-s_0\right)T_j^q\right)}^{-1}{e}_1-{\delta }_{k,j}.$

[0231] Appendix 2: Adjoint Electromagnetic Sensitivity Analysis using Single-Size MPVL Method

[0232] 1) Define κ.sub.k,j, e.sub.1, I.sub.q, {tilde over (G)}.sub.i, {tilde over (C)}.sub.i, K.sub.0, K.sub.1, K.sub.2, {circumflex over (β)}.sub.k,j, and b.sub.j.

[0233] Define frequency s. Define pre-solution frequency s.sub.0.

[0234] Define {w.sub.m}.sub.m=1.sup.q. Set K.sub.s=K.sub.0+s.sub.0K.sub.1+s.sub.0.sup.2K.sub.2.

[0235] 2) Obtain u.sub.0 by solving the linear system K.sub.su.sub.0=b.sub.1.

[0236] Set

[00032]$r_0=\left[\begin{array}{c}{u}_0\\ s_0{u}_0\end{array}\right]\mathrm{and}r=r_0.$

[0237] For m=1, 2, . . . , q do:

[0238] 3) Set

[00033]$v_m=\frac{r}{.Math.r.Math.}\mathrm{and}{\lambda }_m=w_m^T{v}_m.$

[0239] 4) If m>1, set

t.sub.m,m-1=∥r∥.

[0240] 5) Solving the linear system to obtain um,

[00034]$K_s{u}_m={\frac{1}{s_0}\left[\begin{array}{c}{K}_0\\ -s_0{K}_2\end{array}\right]}^T{v}_m.$

[0241] Calculate

[00035]$r=r_m=\left[\begin{array}{c}{u}_m\\ s_0{u}_m\end{array}\right]-\frac{1}{s_0}\left[\begin{array}{cc}{I}_N& 0_N\\ 0_N& 0_N\end{array}\right]v_{m.}$

[0242] 6) For l=1, 2, . . . , m do: Set

[00036]$t_{l,m}=\frac{w_l^{T}r}{{\lambda }_l}\mathrm{and}r=r-t_{l,m}{v}_l.$

[0243] End For

[0244] End For

[0245] 7) Set

T.sub.k.sup.q=[t.sub.l,m].sub.q−q,V.sub.k.sup.q=[v.sub.1v.sub.2 . . . v.sub.q],

x.sub.j.sup.q=∥r.sub.0∥V.sub.j.sup.q(I.sub.q−(s−s.sub.0)T.sub.j.sup.q).sup.−1e.sub.1.

[0246] 8) Repeat 2)-7) by changing b.sub.j to {circumflex over (β)}.sub.k,j to calculate {circumflex over (x)}.sub.k,j.sup.qT

[0247] 9) Calculate

[00037]$\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\hat{x}}_{k,j}^{q^T}\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j^q.$

[0248] Appendix 3: Self-Adjoint Electromagnetic Sensitivity Analysis using Single-Size MPVL Method

[0249] 1) Define κ.sub.k,j, e.sub.1, I.sub.q, {tilde over (G)}.sub.i, {tilde over (C)}.sub.i, K.sub.0, K.sub.1, K.sub.2, b.sub.k, and b.sub.j.

[0250] Define frequency s. Define pre-solution frequency s.sub.0.

[0251] Define {w.sub.m}.sub.m=1.sup.q. Set K.sub.s=K.sub.0+s.sub.0K.sub.1+s.sub.0.sup.2K.sub.2.

[0252] 2) Obtain u.sub.0 by solving the linear system K.sub.su.sub.0=b.sub.j.

[0253] Set

[00038]$r_0=\left[\begin{array}{c}{u}_0\\ s_0{u}_0\end{array}\right]\mathrm{and}r=r_0.$

[0254] For m=1, 2, . . . , q do:

[0255] 3) Set

[00039]$v_m=\frac{r}{.Math.r.Math.}\mathrm{and}{\lambda }_m=w_m^T{v}_m.$

[0256] 4) If m>1, set

t.sub.m,m-1=∥r∥.

[0257] 5) Solving the linear system to obtain um,

[00040]$K_s{u}_m={\frac{1}{s_0}\left[\begin{array}{c}{K}_0\\ -s_0{K}_2\end{array}\right]}^T{v}_m.$

[0258] Calculate

[00041]$r=r_m=\left[\begin{array}{c}{u}_m\\ s_0{u}_m\end{array}\right]-\frac{1}{s_0}\left[\begin{array}{cc}{I}_N& 0_N\\ 0_N& 0_N\end{array}\right]v_m.$

[0259] 6) For l=1, 2, . . . , m do: Set

[00042]$t_{l,m}=\frac{w_l^{T}r}{{\lambda }_l}\mathrm{and}r=r-t_{l,m}{v}_l.$

[0260] End For

[0261] End For

[0262] 7) Set

T.sub.j.sup.q=[t.sub.l,m].sub.q×q,V.sub.j.sup.q=[v.sub.1v.sub.2 . . . v.sub.q],

x.sub.j.sup.q=∥r.sub.0∥V.sub.j.sup.q(I.sub.q−(s−s.sub.0)T.sub.j.sup.q).sup.−1e.sub.1.

[0263] 8) Repeat 2)-7) by changing b.sub.j to b.sub.k to calculate x.sub.k.sup.q.

[0264] 9) Calculate

[00043]$\frac{\partial S_{k,j}}{\partial {\varphi }_i}\approx -{\kappa }_{k,j}{x}_k^{q^T}\left({\tilde{G}}_i+s{\tilde{C}}_i\right)x_j^q.$