A) Moments of the harmonic balance equations

The moments computation algorithm is applied to a general nonlinear system described by its harmonic balance equations [Reference Gad, Khazaka, Nakhla and Griffith19]. The harmonic balance formulation is derived from the modified nodal analysis (MNA) formulation for a nonlinear system in the time domain which is of the form [Reference Ho, Ruehli and Brennan20]

(4)$${\bi Gx}\lpar t\rpar +{\bi Cx}\lpar t\rpar +{\bi f}\lpar {\bi x}\lpar t\rpar \rpar ={\bi b}\lpar t\rpar \comma \;$$

where $x\lpar t\rpar \in {\rm {\opf R}}^n$ is a vector of n unknown variables consisting of node voltages and branch currents. ${\bi G} \in {{{\rm \opf R}}^{n \times n}}$ and ${\bi C} \in {{\rm \opf R}^{n \times n}}$ are matrices representing the contributions of the linear memory less and memory elements, respectively. ${\bi f}\lpar x\rpar \in {{\rm \opf R}^n}$ and ${\bi b}\lpar t\rpar \in {{\rm \opf R}^n}$ are the vectors of nonlinear equations and the independent sources, respectively.

The harmonic balance equations can be obtained from the formulation in (4) by expressing the vector of variables x(t) as a Fourier Series of H harmonics given by:

(5)$$x\lpar t\rpar ={{\bi A}_0}+\sum\limits_{k=1}^H \, \left({{{\bi A}_k}{\rm cos}\lpar {\omega _k}t\rpar +{{\bi B}_k}{\rm sin}\lpar {\omega _k}t\rpar } \right).$$

In this case, the number of variables in the formulation increases to N _h and the harmonic balance equations can be expressed as:

(6)$$\bar{\bi GX}+\bar{\bi CX}+{\bi F}\lpar {\bi X}\rpar ={{\bi B}_{DC}}+\alpha {{\bi B}_{RF}}.$$

In this formulation, ${\bi X} \in {{\rm \opf R}^{{N_h}}}$ is a vector of unknown cosine and sine coefficients for each of the variables in x(t). The vectors ${{\bi B}_{DC}} \in {{\rm \opf R}^{{N_h}}}$ and ${{\bi B}_{RF}} \in {{\rm \opf R}^{{N_h}}}$ contain the contributions of the dc and the RF signal sources, respectively. In this paper, α refers to the amplitude of the input RF voltage signal. $\bar{\bi G} \in {{\rm \opf R}^{{N_h} \times {N_h}}}$ and $\bar{\bi C} \in {{\rm \opf R}^{{N_h} \times {N_h}}}$ are block matrices representing the contributions of the linear memoryless and memory elements, respectively, and ${\bi F}\lpar {\bi X}\rpar \in {{\rm \opf R}^{{N_h}}}$ is the vector of nonlinear equations which also includes the nonlinear charge expressions.

The moments of the system are defined as the coefficients of the Taylor series expansion of the output solution vector X with respect to the input RF amplitude α. This can be expressed as

(7)$${\bi X}={{\bi M}_0}+{{\bi M}_1}\alpha+{{\bi M}_2}{\alpha ^2}+{{\bi M}_3}{\alpha ^3}+\ldots=\sum\limits_{k=0}^\infty \, {{\bi M}_k}{\alpha ^k}\comma \;$$

where M_k is the kth moment vector and can also be expressed as [Reference Gad, Khazaka, Nakhla and Griffith19]:

(8)$${{\bi M}_k}=\mathop {\left. {\displaystyle{1 \over {k!}}\displaystyle{{{\partial ^k}{\bi X}} \over {\partial {\alpha ^k}}}} \right\vert }\nolimits_{\alpha=0} .$$

It is important to note that the Jacobian matrix of the harmonic balance equations given in (6) is typically large and block dense, therefore making the harmonic balance solution computationally expensive. Alternatively, the moments method does not require performing a harmonic balance simulation, but rather computes the moments defined in (7) using a sparse moments matrix, thereby presenting significant computation cost reduction.

B) Moments computation algorithm

The moments computation algorithm has been presented several times in the literature [Reference Achar and Nakhla21–Reference Chiprout and Nakhla23] and has many useful applications in electronic circuit simulation, such as in model order reduction based simulation methods. In this section, an overview of the moments computation algorithm is presented, with an explanation of what makes the algorithm very computationally attractive for distortion analysis applications. For the full details of the algorithm, the reader is invited to refer to the literature [Reference Tannir and Khazaka17, Reference Achar and Nakhla21–Reference Chiprout and Nakhla23].

The moments computation algorithm is about determining the unknown moment vectors, M_k, when the vector of unknown variables X is expressed using (7). Just as is the case with Taylor series, the moments expansion has a finite radius of convergence, thereby limiting their application to circuits that are mildly nonlinear in the signal path. The moments expansion for general mildly nonlinear amplifier circuits is therefore carried out at the DC operating point, which must be computed first as the ‘zero’ moment vector, as will be shown next. To account for stronger nonlinearities that occur outside the RF signal path, such as when simulating mixer circuits with a high power local oscillator, the moments expansion must be carried out at the local oscillator operating point, as shown in [Reference Tannir and Khazaka17].

The main steps of the moments computation algorithm for amplifier circuits are as follows:

• − The zero moment vector, M₀, is obtained by finding the solution of the system described by (6) with the RF amplitude (α) set to zero. This is equivalent to obtaining the dc solution of the system.

• − The first moment vector, M₁, is obtained by using one LU decomposition to solve the relation given by

  (9)$${\bf \Phi }{{\bi M}_1}={{\bi B}_{RF}}\comma \;$$
  where the matrix Φ is the sparse moment computation matrix which has the same structure as a harmonic balance Jacobian, but contains only dc spectral components. It is defined as
  (10)$${\bf \Phi }=\mathop {\left. {\bar{\bi G}+\bar{\bi C}+\displaystyle{{\partial {\bi F}\lpar {\bi X}\rpar } \over {\partial {\bi X}}}} \right\vert }\nolimits_{\alpha=0} .$$
  Note that the matrices $\bar{\bi G}$ and $\bar{\bi C}$ are very sparse, while the Jacobian of the nonlinear vector F(X) is given by
  (11)$$\displaystyle{{\partial {\bi F}\lpar {\bi X}\rpar } \over {\partial {\bi X}}}=\left[{\matrix{ {\displaystyle{{\partial {F_1}} \over {\partial X1}}} & \cdots & {\displaystyle{{\partial {F_1}} \over {\partial Xn}}} \cr \vdots & \ddots & \vdots \cr {\displaystyle{{\partial {F_n}} \over {\partial X1}}} & \cdots & {\displaystyle{{\partial {F_n}} \over {\partial Xn}}} \cr } } \right].$$

  In a standard harmonic balance Jacobian, each $\displaystyle{{\partial {F_j}} \over {\partial Xi}}$ term is, when present, a full block matrix which is not the case with the moments matrix since it is evaluated with the RF amplitude set to zero [Reference Gad, Khazaka, Nakhla and Griffith19].

• − The remaining moment vectors, M_n, are found by solving the following recursive relation

  (12)$${\bf \Phi }{{\bi M}_n}=- n\sum\limits_{j=1}^{n - 1} \, \lpar n - j\rpar {{\bi T}_j}{{\bi M}_{n - j}}.$$

• − In this relation, T_j is a matrix that contains the jth moment of the derivatives of the nonlinear equations.

As can be seen from (9) to (12), the computation of the moment vectors is a solution of a set of linear algebraic equations where the left-hand-side matrix (Φ) is the same throughout. A summary of the moments computation algorithm is given in the flowchart of Fig. 2.

Fig. 2. Moments computation algorithm flowchart.

C) Link between the moment vectors and IP ₃

The input–output relation of a nonlinear circuit can be expressed using a Volterra series as [Reference Maas24]

(13)$$x\lpar t\rpar ={{\bi H}_0}+{{\bi H}_1}\lsqb {v_{in}}\lpar t\rpar \rsqb +{{\bi H}_2}\lsqb {v_{in}}\lpar t\rpar \rsqb +{{\bi H}_3}\lsqb {v_{in}}\lpar t\rpar \rsqb +\ldots\comma \;$$

where H_n is the nth Volterra operator. When the values of the Volterra kernels (${{\bi H}_n}\lpar j{\omega _1}\comma \; \ldots\comma \; j{\omega _n}\rpar $) can be analytically obtained at specific frequency values, we can say that the IP ₃ of amplifier circuits can be determined using

(14)$$I{P_3}=\sqrt {\displaystyle{4 \over 3}\displaystyle{{\vert {{\bi H}_1}\lpar j{\omega _1}\rpar \vert } \over {\vert {{\bi H}_3}\lpar j{\omega _1}\comma \; j{\omega _1}\comma \; - j{\omega _2}\rpar \vert }}} .$$

In [Reference Tannir and Khazaka16–Reference Tannir and Khazaka18], it was shown that the kernels required to determine IP ₃ according to (14) were obtained numerically from the first and third moment vectors. This allowed us to express the relation for computing IP ₃ as

(15)$$I{P_3}=\sqrt {\displaystyle{{\vert {{\bi M}_1}\lsqb {\omega _1}\rsqb \vert } \over {\vert {{\bi M}_3}\lsqb 2{\omega _1} - {\omega _2}\rsqb \vert }}}\comma \;$$

where M₁[ω ₁] is the numerical entry in the first moment vector at the fundamental frequency, and M₃[2ω ₁ − ω ₂] is the numerical entry in the third moment vector at the third-order intermodulation frequency.

A) Extended formulation of system equations

Consider a nonlinear circuit excited by one or more input tones and containing a MEMS variable capacitor. The electrical steady-state solution can be obtained using the harmonic balance approach by expressing the periodic solution as a truncated series of sine and cosine functions at the harmonics of the inputs as well as at the intermodulation products. This results in a set of nonlinear algebraic equations in the form given by (6). In order to account for the mechanical characteristics of the MEMS device, this formulation has to be expanded to include the MEMS mechanical dynamic model nonlinearities. Therefore, additional variables for the charge of the nonlinear capacitor (q(t)) and the displacement of the membrane (z(t)) must be accounted for alongside the regular unknown node voltages and currents of the MNA formulation [Reference Ho, Ruehli and Brennan20]. For a general system containing a MEMS variable capacitor, the set of variables will therefore include, among other variables,

(16)$${x_1}\lpar t\rpar =V\lpar t\rpar \comma \;$$

(17)$${x_2}\lpar t\rpar =z\lpar t\rpar \comma \;$$

where x ₁(t) and x ₂(t) will be two variables in the unknown solution vector x(t). The charge Q defined in (2) is also expressed as an unknown variable, q(t), in the solution vector x(t), as a function of the applied voltage and the displacement as given by

(18)$${x_3}\lpar t\rpar =q\lpar t\rpar =q\lsqb {x_1}\lpar t\rpar \comma \; {x_2}\lpar t\rpar \rsqb .$$

The electromechanical equation in (1) can now be rewritten as a function the variables x ₁(t) and x ₂(t) as:

(19)$$f\lpar t\rpar \equiv f\left[{{x_1}\lpar t\rpar \comma \; {x_2}\lpar t\rpar \comma \; \displaystyle{{d{x_2}\lpar t\rpar } \over {dt}}\comma \; \displaystyle{{{d^2}{x_2}\lpar t\rpar } \over {d{t^2}}}} \right]=0.$$

This new set of nonlinear dynamic equations will complement the regular electric circuit equations for all the MEMS capacitors in the system. It is important to observe from these relations that difficulties in performing a full simulation will arise due to the presence of singularities at z = d. In [Reference Rizzoli, Gaddi, Iannacci, Masotti and Mastri13], it is shown that an introduction of a carefully selected additional state variable x _A(t) is sufficient to remove all singularities.

The combination of the new MEMS variables and expressions of (16–19) with the regular MNA formulation of nonlinear electronic components given by (4) results in a new set of second-order differential algebraic equations as follows

(20)$${\bi Gx}\lpar t\rpar +{\bi C\dot x}\lpar t\rpar +{\bi S\ddot x}\lpar t\rpar +{\bi f}\lpar {\bi x}\lpar t\rpar \rpar ={\bi b}\lpar t\rpar \comma \;$$

where x(t) is the vector of unknown variables and can be expressed using a Fourier Series as:

(21)$${\bi x}\lpar t\rpar ={{\bi A}_0}+\sum\limits_{k=1}^H \, \left({{{\bi A}_k}{\rm cos}\lpar {\omega _k}t\rpar +{{\bi B}_k}{\rm sin}\lpar {\omega _k}t\rpar } \right).$$

Similarly, the first and second derivatives of x(t) with respect to time can also be expressed as:

(22)$${\bi \dot x}\lpar t\rpar =\sum\limits_{k=1}^H \, \left({{\omega _k}{{\bi B}_k}{\rm cos}\lpar {\omega _k}t\rpar - {\omega _k}{{\bi A}_k}{\rm sin}\lpar {\omega _k}t\rpar } \right)\comma \;$$

(23)$${\bi \ddot x}\lpar t\rpar =\sum\limits_{k=1}^H \, \left({ - \omega _k^2 {{\bi A}_k}{\rm cos}\lpar {\omega _k}t\rpar - \omega _k^2 {{\bi B}_k}{\rm sin}\lpar {\omega _k}t\rpar } \right)\comma \;$$

where ω ₁ → ω _k are the harmonics of the operating frequency. The matrices G, C, f(x(t)) and b(t) are identical to those defined in (14). The new matrix ${\bi S} \in {{\rm \opf R}^{n \times n}}$ contains the coefficients of the second-order differentials associated with the MEMS equation of motion and is of the same dimensions as G and C.

The expression in (20) should be expressed in the frequency domain in order to be able to develop the moments computation algorithm. The new set of equations is therefore defined as

(24)$${\bf \Omega }{\bi X}+{\bi F}\lpar {\bi X}\rpar ={{\bi B}_{DC}}+\alpha {{\bi B}_{RF}}\comma \;$$

where the matrix ${\bf \Omega }=\bar{\bi G}+\bar{\bi C}+\bar{\bi S}$. The structures of the $\bar{\bi G}$, $\bar{\bi C}$, and $\bar{\bi S}$ matrices are as follows. $\bar{\bi G} \in {{\rm \opf R}^{{N_h} \times {N_h}}}$ is a block matrix $\bar{\bi G}=\left[{{{\bi G}_{ij}}} \right]$ representing the contribution of the linear memoryless elements of the network to the frequency components. The blocks ${{\bi G}_{ij}} \in {{\rm \opf R}^{{N_b} \times {N_b}}}$ are diagonal matrices given by

(25)$${{\bi G}_{ij}}={\rm diag}\lpar {g_{ij}}\comma \; \cdots\comma \; {g_{ij}}\rpar \comma \;$$

with g _ij being the corresponding (i, j)th entry in the G matrix in (6). Similarly, $\bar{\bi C} \in {{\rm \opf R}^{{N_h} \times {N_h}}}$ is a block matrix $\bar{\bi C}=\left[{{{\bi C}_{ij}}} \right]$ representing the contribution of the linear memory elements of the network to the frequency components. The blocks ${{\bi C}_{ij}} \in {{\rm \opf R}^{{N_b} \times {N_b}}}$ are diagonal matrices given by

(26)$${{\bi C}_{ij}}={c_{i\comma j}}\left[{\matrix{ 0 & 0 & 0 & \cdots & 0 & 0 \cr 0 & 0 & {{\omega _1}} & \cdots & 0 & 0 \cr 0 & { - {\omega _1}} & 0 & \cdots & 0 & 0 \cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \cr 0 & 0 & 0 & \cdots & 0 & {{\omega _k}} \cr 0 & 0 & 0 & \cdots & { - {\omega _k}} & 0 \cr } } \right]\comma \;$$

with c _ij being the corresponding entry in the C matrix in (20). Finally the matrix $\bar{\bi S} \in {{\opf R}^{{N_h} \times {N_h}}}$ is a block matrix $\bar{\bi S}=\left[{{{\bi S}_{ij}}} \right]$ with each ${{\bi S}_{ij}} \in {{\rm \opf R}^{{N_b} \times {N_b}}}$ block being a diagonal matrix with the following structure

(27)$${{\bi S}_{ij}}={s_{i\comma j}}\left[{\matrix{ 0 & 0 & 0 & \cdots & 0 & 0 \cr 0 & { - \omega _1^2 } & 0 & \cdots & 0 & 0 \cr 0 & 0 & { - \omega _1^2 } & \cdots & 0 & 0 \cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \cr 0 & 0 & 0 & \cdots & { - \omega _k^2 } & 0 \cr 0 & 0 & 0 & \cdots & 0 & { - \omega _k^2 } \cr } } \right]\comma \;$$

with s _ij being the corresponding entry in the S matrix in (5). Note that N _b varies with the number of harmonics (H) and is equal to 2H + 1.

B) Moments computation for MEMS circuits

In this section, we derive the moments computation algorithm using the new formulation given in the previous subsection. In this way, we show how to compute the moments for circuits that contain both the electrical nonlinearities associated with semiconductor circuits, in addition to the nonlinearities associated with MEMS variable capacitors. The moments expansion will be performed around the DC operating point for the same reasons that were previously discussed in Section III.

The derivation of the moments computation algorithm begins by substituting (7) into (24), which results in the following expression:

(28)$$\Omega \sum\limits_{k=0}^q \, {{\bi M}_k}{\alpha ^k}+\sum\limits_{k=0}^q \, {{\bi F}_k}{\alpha ^k} - {{\bi B}_{DC}} - \alpha {{\bi B}_{RF}}=0.$$

The terms F_k are the Taylor expansion coefficients of F(X) with respect to α given by:

(29)$${\bi F}\lpar {\bi X}\rpar =\sum\limits_{k=0}^q \, {{\bi F}_k}{\alpha ^k}.$$

It is useful to also express the derivative of this nonlinear vector with respect to the solution vector X (i.e. the Jacobian matrix $\displaystyle{{\partial F\lpar X\rpar } \over {\partial X}}$) as a Taylor series expansion with respect to α as follows

(30)$${\bi T}\lpar \alpha \rpar =\displaystyle{{\partial {\bi F}\lpar {\bi X}\rpar } \over {\partial {\bi X}}}=\sum\limits_{k=0} \, {{\bi T}_k}{\alpha ^k}\comma \;$$

where T_n are referred to as the moments of the nonlinear Jacobian matrix. The zero moment vector, M₀, is obtained by setting the value of α in (28) to zero which yields the relation:

(31)$$\Omega {{\bi M}_0}+{\bi F}\lpar {{\bi M}_0}\rpar ={{\bi B}_{DC}}.$$

Note that the solution to equation (31) is essentially the computation of the dc solution for the circuit, which has a relatively small computation cost. To solve for the remaining moments (M_n; n ≥ 1), coefficients of equal powers of α are equated on both sides of (28). Equating the coefficients of the first power of α results in:

(32)$$\Omega {{\bi M}_1}+{{\bi F}_1}={{\bi B}_{RF}}.$$

It is useful here to rewrite ${{\bi F}_1}=\displaystyle{{\partial F} \over {\partial \alpha }}\left\vert _{\alpha=0} \right. {\rm as}\, {F_1}=\displaystyle{{\partial F} \over {\partial X}} \cdot \displaystyle{{\partial X} \over {\partial \alpha }} _{\vert \alpha=0}={{\bi T}_0}{{\bi M}_1}$. Substituting this expression into (32) then yields:

(33)$$\underbrace {\lpar {\bf \Omega }+{{\bi T}_0}\rpar }_{\bf \Psi }{{\bi M}_1}={{\bi B}_{RF}}.$$

The first moment can now be obtained using one LU decomposition to solve (33). It is important to note that the matrix Ψ = (Ω + T₀) is a sparse matrix which is already evaluated when obtaining the dc solution of the nonlinear system. Therefore, it has a similar structure to a Jacobian matrix, but is significantly more sparse. In addition, a comparison of Ψ in (33) with Φ in (10) shows that the sparsity pattern of the original moments computation matrix is preserved for systems that contain MEMS nonlinearities. To obtain the remaining moments, the coefficients of α ⁿ on both sides of (28) are equated to obtain:

(34)$$\Omega {{\bi M}_n}+{{\bi F}_n}=0\quad \quad n\gt 1.$$

The evaluation of F_n in (34) remains unchanged from the original moments computation algorithm described in Section III.B and is given by:

(35)$${{\bi F}_n}={{\bi T}_0}{{\bi M}_n}+n\sum\limits_{j=1}^{n - 1} \, \lpar n - j\rpar {{\bi T}_j}{{\bi M}_{n - j}}.$$

Combining the expression in (35) with that of the relation in (34) results in the following recursive relation for computing the nth moment vector:

(36)$$\underbrace {\lpar {\bf \Omega }+{{\bi T}_0}\rpar }_{\bf \Psi }{{\bi M}_n}=- n\sum\limits_{j=1}^{n - 1} \, \lpar n - j\rpar {{\bi T}_j}{{\bi M}_{n - j}}.$$

The T_j terms are as defined in (30). Note that since F(X) and X are vectors, T_j will be block matrices of the form

(37)$${{\bi T}_j}=\mathop {\displaystyle{{\partial {\bi F}\lpar {\bi X}\rpar } \over {\partial {\bi X}}}}\nolimits_j=\left[{\matrix{ {\displaystyle{{\partial F1} \over {\partial X_{1j}}}} & \cdots & {\displaystyle{{\partial F1} \over {\partial X_{nj}}}} \cr \vdots & \ddots & \vdots \cr {\displaystyle{{\partial Fn} \over {\partial X_{1j}}}} & \cdots & {\displaystyle{{\partial F_n} \over {\partial X_{nj}}}} \cr } } \right]\comma \;$$

where each $\displaystyle{{\partial F_k} \over {\partial X_i}}$ term is a matrix in itself. For example, if we consider the first block matrix $\displaystyle{{\partial F_1} \over {\partial X_i}}$, then each coefficient of its Taylor series expansion with respect to α is given by

(38)$$\displaystyle{{\partial F_1} \over {\partial X_1\, j}}={\Gamma ^{ - 1}}\left[{\matrix{ {\displaystyle{{\partial {f_1}\lpar {x_1}\lpar {t_1}\rpar \rpar } \over {\partial x_1\, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,}}}{_j} & {} & 0 \cr {} & \ddots & {} \cr 0 & {} & {\displaystyle{{\partial f_1\lpar {x_1}\lpar {t_s}\rpar \rpar } \over {\partial x_1}}}{_j} \cr } } \right]\Gamma\comma \;$$

where t ₁ to t _s are time sample points that are equally spaced over the fundamental period. Note that frequency mapping and truncation methods [Reference Kundert, White and Sangiovanni-Vincentelli11] are used in order to handle quasi-periodic inputs efficiently using the Fast Fourier Transform, and Γ is the Inverse Direct Fourier Transform matrix. Also note that the matrix vector multiplication with Γ can be done efficiently by taking advantage of the Fast Fourier Transform algorithm. Once evaluated, the Taylor coefficient $\displaystyle{{\partial F1} \over {\partial X1j}}$ is entered in T_j at the location corresponding to $\displaystyle{{\partial F1} \over {\partial X1\, \, \, \, }}$.

From (37) and (38), it can be seen that the evaluation of the T_j matrices reduces to the evaluation of the moments for the derivatives of the nonlinear functions, i.e. the evaluation of $\displaystyle{{\partial f\lpar x\rpar } \over {\partial x\, j}}$ for each function f and each variable x. Next, we show how these derivatives are determined for the MEMS nonlinearities.

C) Computation of the derivatives for the MEMS nonlinearities

The computation of the moments requires the computation of the coefficients of the series expansions of the partial derivatives for each of the nonlinear functions found in the nonlinear MEMS dynamical model described by (1–3) so they can be evaluated using (38). Note that the derivation of these expressions is done only once, and are then stored in the simulator to then be evaluated numerically in an efficient manner going forward.

The two primary nonlinear functions of interest are:

(39)$${F_1}\lpar Q\rpar =\displaystyle{{Q^2} \over {2{{\rm \epsilon }_0}A}}\comma \;$$

(40)$${F_2}\lpar V\comma \; z\rpar =\displaystyle{{{{\rm \epsilon }_0}A} \over {d - z}}V.$$

Using the variable assignment of Section IV.A., these nonlinear functions can be expressed as:

(41)$${f_1}=\displaystyle{{x_3^2 } \over {2{{\rm \epsilon }_0}A}}\comma \;$$

(42)$${f_2}=\displaystyle{{{{\rm \epsilon }_0}A} \over {d - {x_2}}}{x_1}.$$

We need to determine the derivatives of these functions with respect to the variables, x _1,x _2, and x ₃, in addition to the required expressions for the evaluation of their series expansion coefficients. To simplify the presentation, we define a new variable g ^(m,n) which is the partial derivative of function f _m with respect to variable x _n. With this new variable, the task at hand becomes determining the expressions required for computing each $g_j^{\lpar m\comma n\rpar } $ term, which are the coefficients of the series expansions of g ^{(m, n)} with respect to the RF amplitude α defined as:

(43)$${g^{\lpar 1\comma 3\rpar }}=\displaystyle{{\partial {f_1}} \over {\partial x3}}=\displaystyle{{{x_3}} \over {{{\rm \epsilon }_0}A}}=\sum\limits_{i=0} \, g_i^{\lpar 1\comma 3\rpar } {\alpha ^i}\comma \;$$

(44)$${g^{\lpar 2\comma 1\rpar }}=\displaystyle{{\partial {f_2}} \over {\partial x1}}=\displaystyle{{{{\rm \epsilon }_0}A} \over {\lpar d - {x_2}\rpar }}=\sum\limits_{i=0} \, g_i^{\lpar 2\comma 1\rpar } {\alpha ^i}\comma \;$$

(45)$${g^{\lpar 2\comma 2\rpar }}=\displaystyle{{\partial {f_2}} \over {\partial x2}}=\displaystyle{{{{\rm \epsilon }_0}A} \over {{{\lpar d - {x_2}\rpar }^2}}}{x_1}=\sum\limits_{i=0} \, g_i^{\lpar 2\comma 2\rpar } {\alpha ^i}.$$

Note that the variables x ₁(t), x ₂(t), and x ₃(t) can also be represented using series expansions in the time domain as:

(46)$${x_1}=\sum\limits_{i=0} \, m_i^{\lpar 1\rpar } {\alpha ^i}\comma \;$$

(47)$${x_2}=\sum\limits_{i=0} \, m_i^{\lpar 2\rpar } {\alpha ^i}\comma \;$$

(48)$${x_3}=\sum\limits_{i=0} \, m_i^{\lpar 3\rpar } {\alpha ^i}.$$

We begin the derivations with the first function shown in (41). A comparison of this function with its partial derivative expression given in (43) allows us to express f ₁ in terms of g ^(1,3) and x ₃ as follows:

(49)$${f_1}=\displaystyle{1 \over 2}{x_3}{g^{\lpar 1\comma 3\rpar }}.$$

Taking the derivative of (49) with respect to α then yields the expression:

(50)$$\displaystyle{{\partial {f_1}} \over {\partial \alpha }}=\displaystyle{1 \over 2}{g^{\lpar 1\comma 3\rpar }}\displaystyle{{\partial {x_3}} \over {\partial \alpha }}+\displaystyle{1 \over 2}\displaystyle{{\partial {g^{\lpar 1\comma 3\rpar }}} \over {\partial \alpha }}.$$

At this point, we make use of the chain rule so we can express $\displaystyle{{\partial {f_1}} \over {\partial \alpha }}$ as:

(51)$$\displaystyle{{\partial {f_1}} \over {\partial \alpha }}=\displaystyle{{\partial {f_1}} \over {\partial x3}}\displaystyle{{\partial {x_3}} \over {\partial \alpha }}={g^{\lpar 1\comma 3\rpar }}\displaystyle{{\partial {x_3}} \over {\partial \alpha }}.$$

We can then substitute (51) into (50) and rearrange to obtain:

(52)$${g^{\lpar 1\comma 3\rpar }}\displaystyle{{\partial {x_3}} \over {\partial \alpha }}=\displaystyle{1 \over 2}{g^{\lpar 1\comma 3\rpar }}\displaystyle{{\partial {x_3}} \over {\partial \alpha }}+\displaystyle{1 \over 2}{x_3}\displaystyle{{\partial {g^{\lpar 1\comma 3\rpar }}} \over {\partial \alpha }}\comma \;$$

(53)$${x_3}\displaystyle{{\partial {g^{\lpar 1\comma 3\rpar }}} \over {\partial \alpha }}={g^{\lpar 1\comma 3\rpar }}\displaystyle{{\partial {x_3}} \over {\partial \alpha }}.$$

Finally, we replace g ^(1,3) and x ₃ in (53) with their series expansions found in (43) and in (48) to obtain:

(54)$$\sum\limits_{i=0} \, {m_i}{\alpha ^i}\sum\limits_{i=1} \, ig_i^{\lpar 1\comma 3\rpar } {\alpha ^{i - 1}}=\sum\limits_{i=0} \, g_i^{\lpar 1\comma 3\rpar } {\alpha ^i}\sum\limits_{i=1} \, i{m_i}{\alpha ^{i - 1}}.$$

From this expression, each $g_j^{\lpar 1\comma 3\rpar } $ term required for the moments algorithm as given by (38) is found accordingly. The general expressions for evaluating $g_n^{\lpar 1\comma 3\rpar } $ are therefore:

(55)$$g_0^{\lpar 1\comma 3\rpar }=\displaystyle{{m_0^{\lpar 3\rpar } } \over {{{\rm \epsilon }_0}A}}\comma \;$$

(56)$$g_n^{\lpar 1\comma 3\rpar }=\displaystyle{1 \over {nm_0^{\lpar 3\rpar } }}\left({ng_0^{\lpar 1\comma 3\rpar } m_n^{\lpar 3\rpar }+\sum\limits_{i=1}^{n - 1} \, g_{n - i}^{\lpar 1\comma 3\rpar } m_i^{\lpar 3\rpar } \lpar 2i - n\rpar } \right).$$

A fundamentally similar procedure is followed for the second function f ₂ shown in (42). The difference here is that f ₂ is a function of two variables and therefore we need to determine the partial derivatives with respect to both x ₁ and x ₂. The relation between f ₂ and the derivatives g ^(2,1) and g ^(2,2) are:

(57)$${f_2}={x_1}{g^{\lpar 2\comma 1\rpar }}\comma \;$$

(58)$${f_2}=\lpar d - {x_2}\rpar {g^{\lpar 2\comma 2\rpar }}.$$

Notice that the partial derivative with respect to x ₁ results in a constant term and therefore only $g_0^{\lpar 2\comma 1\rpar } $ will be nonzero in value. As for determining the $g_j^{\lpar 2\comma 2\rpar } $ terms, we start by taking the derivative of (58) with respect to α which yields

(59)$$\displaystyle{{\partial {f_2}} \over {\partial \alpha }}=d\displaystyle{{\partial {g^{\lpar 2\comma 2\rpar }}} \over {\partial \alpha }} - \left({{g^{\lpar 2\comma 2\rpar }}\displaystyle{{\partial {x_2}} \over {\partial \alpha }}+{x_2}\displaystyle{{\partial {g^{\lpar 2\comma 2\rpar }}} \over {\partial \alpha }}} \right).$$

By once again employing the chain rule, $\displaystyle{{\partial {f_2}} \over {\partial \alpha }}$ can be expressed as

(60)$$\displaystyle{{\partial {f_2}} \over {\partial \alpha }}=\displaystyle{{\partial {f_2}} \over {\partial x2}}\displaystyle{{\partial {x_2}} \over {\partial \alpha }}={g^{\lpar 2\comma 2\rpar }}\displaystyle{{\partial {x_2}} \over {\partial \alpha }}\comma \;$$

which can then be substituted into the expression of (59) to obtain:

(61)$$2{g^{\lpar 2\comma 2\rpar }}\displaystyle{{\partial {x_2}} \over {\partial \alpha }}=\lpar d - {x_2}\rpar \displaystyle{{\partial {g^{\lpar 2\comma 2\rpar }}} \over {\partial \alpha }}.$$

Replacing the variables x ₂ and g ^(2,2) in (61) by their series expansions given in (45) and (47) we then obtain:

(62)$$\left({d - {m_i}\sum\limits_{i=0} \, {\alpha ^i}} \right)\sum\limits_{i=1} \, i{g_i}{\alpha ^{i - 1}}=2\sum\limits_{i=0} \, {g_i}{\alpha ^i}\sum\limits_{i=1} \, i{m_i}{\alpha ^{i - 1}}.$$

This relation is used to determine the individual $g_j^{\lpar 2\comma 2\rpar } $ terms required for determining the moments, which are evaluated using:

(63)$$g_0^{\lpar 2\comma 2\rpar }=\displaystyle{{{{\rm \epsilon }_0}A} \over {{{\lpar d - m_0^{\lpar 2\rpar } \rpar }^2}}}m_0^{\lpar 1\rpar }\comma \;$$

(64)$$g_n^{\lpar 2\comma 2\rpar }=\displaystyle{1 \over {n\lpar d - m_0^{\lpar 2\rpar } \rpar }}\left({2ng_0^{\lpar 2\comma 2\rpar } m_n^{\lpar 2\rpar }+\sum\limits_{i=1}^{n - 1} \, g_i^{\lpar 2\comma 2\rpar } m_{n - i}^{\lpar 2\rpar } \lpar 2n - i\rpar } \right).$$

Note that the CPU cost of computing these derivatives is negligible when compared to the overall algorithm. The expressions for the derivatives are stored in the simulator along with the derivatives of other nonlinear expressions present in the circuit model, such as diode nonlinearities, MOSFET saturation current voltage relations and others. Note that a summary of the relations to determine the moment expansions of common standard functions present in most nonlinear expressions is provided in Table 2 [Reference Griffith and Nakhla25].

Table 2. Formulas for derivatives of standard functions.

Finally, we remind the reader that once the moment vectors have been evaluated, the computation of a figure of merit like IP ₃ is accomplished according to the relation given by (15) using the values of numerical entries at specific locations in the moments which we have now shown how to compute. A summary of the main steps of the algorithm is given in Fig. 3.

Fig. 3. Summary of the IP ₃ computation algorithm using moments.

D) Effect of using a more refined MEMS model

In this paper, we have developed the proposed moments method for efficient distortion analysis of circuits containing MEMS variable capacitors on the basis of the simplified 1-D MEMS model that was presented in Section II. For studying intermodulation distortion in addition to the effects of Q and V on the dynamic response of MEMS switches, the 1-D model works well and provides a lot of useful information. The simplified model, however, has its limitations. For a more accurate analysis of intermodulation distortion in addition to other phenomena, such as pull-down, power handling and noise, more refined models can be used. These refinements can include, for example, an expression that accounts for the nonlinear behavior of the spring in the mass–spring–damper model [Reference Innocent, Wambacq, Donnay, Tilmans, Sansen and Man8]. Another modification can be made to include the effects of the contact force, F _c, between the moving plate and the dielectric layer, in addition to the electrostatic force [Reference Rebeiz3]. Similarly, we can also include a general nonlinear expression to more accurately describe the value of the damping coefficient, b, in the refined model [Reference Bao, Yang, Sun and French26]. These additional refinements would come at the expense of additional model complexities, which will mainly appear when determining the additional expressions for the derivatives of the MEMS nonlinearities required to compute the moments. This can be accomplished with the aid of Table 2 for most nonlinear expressions, with the resulting derivatives stored in the simulator as part of the device model.

If the additional expressions in the refined models are only functions of the existing variables in the formulation, the computation time of the proposed moments method will not be significantly affected since the size of the system formulation will remain the same. On the other hand, if new variables are introduced, though the computation time will increase, this will also be equally the case with the reference harmonic balance method since both approaches are based on similar formulations, thereby preserving the CPU cost advantage of the proposed method.