Journal of Sensors
Volume 2015 (2015), Article ID 125621, 10 pages
http://dx.doi.org/10.1155/2015/125621

Research Article

A New MEMS Stochastic Model Order Reduction Method: Research and Application

Bian Xiangjuan,¹ Youping Gong,² Chen Guojin,² and Lv Yunpeng²

¹School of Faculty & Technology, Zhejiang International Studies University, Hangzhou 310012, China
²School of Mechanical Engineering Hangzhou Dianzi University, Hangzhou 310018, China

Received 18 November 2014; Revised 16 March 2015; Accepted 16 March 2015

Academic Editor: Hiroki Kuwano

Copyright © 2015 Bian Xiangjuan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Modeling and simulation of MEMS devices is a very complex tasks which involve the electrical, mechanical, fluidic, and thermal domains, and there are still some uncertainties that need to be accounted for during the robust design of MEMS actuators caused by uncertain material and/or geometric parameters. According to these problems, we put forward stochastic model order reduction method under random input conditions to facilitate fast time and frequency domain analyses; the method makes use of polynomial chaos expansions in terms of the random input variables for the matrices of a finite element model of the system and then uses its transformation matrix to reduce the model; the method is independent of the MOR algorithm, so it is seamlessly compatible with MOR method used in popular finite element solvers. The simulation results verify the method is effective in large scale MEMS design process.

1. Introduction

MEMS are attractive for many applications because of their small size and weight, which allow systems to be miniaturized [1]. Modeling and simulation is at the basis of the prediction of the device behavior and optimization of its performance. The complexity of the simulation analysis adds to the complexity of the model. MEMS characterization often requires computationally expensive analysis such as transient analysis. Moreover, the derived models are often nonlinear. The primary source of nonlinearity is the electrostatic force, which couples the electrical and mechanical energy domains. Simulation of nonlinear models makes use of iterative time consuming algorithms. To improve the simulation computation efficiency, MOR (model order reduction) method was used to extract the lower order ODE system that reproduces the input/output behavior with good accuracy. The method is thus solely based on the mathematical properties of the original system and is, therefore, formal, robust, and in great part automatic [2, 3]. These properties render the use of mathematical model order reduction more and more popular in the study of MEMS devices. Variations during fabrication lead to uncertain material and/or geometric parameters causing a significant impact on MEMS device performance [4]. While there have been significant advances in numerical simulation methods that allow better understanding of the underlying multiphysics [5–7], they, however, assume that the geometrical and physical properties of the device are known in a deterministic sense. Recently there have been efforts towards developing computational methods that can handle input uncertainties. A standard way to deal with these uncertainties is brute-force Monte Carlo simulations; these basically involve considering a large number (typically greater than 10000) of realizations (or samples) of the geometry and solving the deterministic problem for each one of these realizations [8]. Several techniques have been developed, for improving convergence, for example, Latin hypercube sampling, [9, 10], the quasi-Monte Carlo (QMC) method, and the Markov chain Monte Carlo method [11]. The statistics such as the mean and standard deviation of the required output quantity such as the electrostatic force are then generated. This is a very time consuming process. Also, there is some uncertainty research work about MEMS: Kong et al. [12] studied the performance variability of a ceramic MEMS actuator under random variations in the shape of the actuator and the air gap in the condenser. Han and Kwak [13] presented the use of robust optimization during the design of a microgyroscope using MC simulations to compare predicted yields. Liu et al. [14] presented a robust design method to minimize the sensitivity of a laterally vibrating resonator against width variations due to fabrication errors. Stochastic model order reduction can offer an efficient way for optimization of dynamic problems [15–17]. This has led to a lot of efforts towards developing MOR algorithms for variational and parametric uncertainty analysis. In [18], a variational balanced truncation method was introduced for model reduction of variable geometry interconnects. In [19], a method for model order reduction of RCL interconnects was described. Various algorithms for parametric model order reduction have been proposed in [20–25].

In this paper, we input parameter uncertainty. We put forward a new method for MEMS stochastic MOR computing process in the presence of input parameter uncertainty. Firstly, we expressed the random input variables as standard random variables by KL series expansion method; then, we expressed the random output variables as standard random variable polynomials by PC method; while obtaining the polynomials of output variables, performing deterministic MOR for each system to generate the corresponding transformation matrix at each point on the Smolyak sparse grid, at last, calculate the mean, standard deviation and other statistics of the system response. The remainder of this paper is organized as follows. In the next section, we represent random input variables matrix of MEMS models; Section 3 is dedicated to the PC representation of the output random variables; polynomial chaos expansion coefficient solution by sparse grid method was described in Section 4; stochastic reduced order model is described in Section 5. The computational gain provided by the sparse PC expansions simulation is illustrated in Section 6 by numerical examples.

2. Representing Input Variable Stochastic of MEMS Model

2.1. MEMS Stochastic Variables

The probability of the random events can be described as a positive number less than 1; the number which describes the degree of random event is called the sample space Ω, denoted as . Probability measure events may occur only giving a measure of the standard and, for the sample point of the sample space, a similar standard is still not given. It has two basic futures: (1) MEMS stochastic variable is a real single-valued function of sample points; (2) to any real number , is a random event. When we get the definition of random variable, then a basic event can be represented by a determined value, denoted as ; and any event can use domain value of random variables; we can note the following: . The likelihood of the event occurring with probability can be described as .

2.2. MEMS Stochastic Process

MEMS stochastic process refers to a cluster of random variables defined in the parameter set; every point at the parameter set is corresponding to a random variable. One-dimensional stochastic process can be viewed as a natural extension of random vectors.

If a stochastic process was denotes as , to describe its probabilistic properties, the most important aspect is the distribution function of random variables at time ; we denote the function as

If we consider the relationship of the stochastic process between any two random variables, then we can get the following expression:

And if there are variables in MEMS system, then the stochastic process can be described as

2.3. MEMS Stochastic Field

MEMS stochastic field is the natural expansion of stochastic process in special field; to a stochastic field, its basic parameters is special variable , so we can define the stochastic field as series of random variables on a field set. And there are corresponding variables at every point of the field set. In fact, we just take into consideration special variables as basic parameter’s stochastic field denoted as , where is defined area. To a uniform stochastic field, because of its mean function is a constant, so it can be turned into the following form:

where is mean field function and is zero-mean stochastic field. Obviously, the covariance function of is the same as the correlation function of . And we can use limited distribution function to represent the probability of random structure, like the following expression:

where is value limit of variable.

2.4. Stochastic Input Variables of MEMS Model

The stochastic problem of MEMS can be described as the following equation:

where represent the solution of response-problem; is random force; is the MEMS material stochastic field; and is the stochastic event. From expression (5), we discrete the stochastic field and realize the conversion from stochastic field to discrete random variable set. Then we will use KL expansion to make these variables independent from each other. The KL expansion can be written as

where is the mean of the random process, and are correlation function, are the eigenvalues and eigenfunctions, respectively, where and are the spatial coordinates. forms a set of uncorrelated random variables. and form the eigenvector-eigenvalue pair of the covariance kernel that satisfy

At the same time, there is a relationship between the cofeature value function of the variance function:

In practice the expansion in (7) is truncated after a finite number of terms , which leads to a truncation error . As compared to other expansion methods, which use some orthonormal functions , the KL expansion is optimal in the sense that the mean-square error is minimized.

3. MEMS Output Variables by Polynomial Chaos Expansion

The polynomial chaos expansion is a spectral expansion of the random process in terms of the orthogonal polynomials in multidimensional random variables. Let be a set of orthonormal Gaussian random variables. Using this, the polynomial chaos expansion of a second-order output variables random process can be described as

where denotes the polynomial chaos of order in terms of the multidimensional Gaussian random variables. For convenience, (11) can be written as

where there is a one-to-one correspondence between the functions and and between the coefficients. Finite number of random variables is used in the expansion to represent finite number of random parameters in the system. Also, the order of the polynomial used in the expansion is restricted to . Thus, the expansion in (12) can now be written as

The total number of terms included in the polynomial chaos expansion depends on both the dimensionality and the highest order of the multidimensional polynomials used, which can be given as

The eigenvalues and coefficients of eigenvectors are random processes; they belong to ; we write the random eigenvalues and eigenvectors using the Karhunen-Loeve truncated expansion with Hermite polynomial chaos:

where are coefficients of the truncated expended, is family of Hermite polynomial chaos. Introducing a classical normalization condition, , we can get

We also can express the other input random element of the problem with the same approximation on .

4. Polynomial Chaos Expansion Coefficient Solution by Sparse Grid Method

The most important aspect of polynomial chaos method is to calculate the each coefficient of the polynomial chaos expansion. On a multidimensional model containing random variables polynomial chaos unfolds; the output at some point can be calculated with its polynomial coefficients chaos which unfolds. In a vector space expansion, if we us as base, and every sample has corresponding point, we call these points as configuring points. There are many methods to compute the configuring points, such as tensor product method and Stroud-2 method. We use sparse grid method which derived from Smolyak method. This approach can significantly reduce the required number of points and this algorithm provides a linear combination of the tensor product of the interpolation error which can be made very close to full tensor product method. Firstly, let: ; then

where is a linear function and is determined by the limited point function value. When and especially, there is

And all polynomial chaos expansion coefficient solution is just processed by these grids points, and the point’s sets are

where the number of points used in each dimension vary depending on the method of polynomials and quadrature.

5. Stochastic Reduced Order Model

The stochastic reduced order model method’s basic theory is when we get the each point on the Smolyak sparse grid; then we can calculate the full-finite element system matrices at the point, and we can perform deterministic MOR for each system to generate the corresponding transformation matrix. Using this matrix information, we can compute the augmented system using the coefficients.

5.1. Deterministic Reduced Order Model

MEMS devices can be modeled by partial differential equations (PDEs). To simulate such models, spatial discretization via, for example, finite element discretization is necessary, which results in a system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs). After spatial discretization, the number of degrees of freedom usually is very high. Therefore, it is time consuming to simulate the large-scale systems of ODEs or DAEs. Developed from well-established mathematical theories and robust numerical algorithms, model order reduction (MOR) has been recognized as being very efficient in reducing the simulation time of large-scale systems. Through MOR, a small system of ODEs with reduced number of equations (reduced model) is derived. The reduced model is simulated instead, and the solution of the original PDEs or ODEs can then be recovered from the solution of the reduced model. As a result, the simulation time of the original large-scale system can be shortened by several orders of magnitude. The reduced model as a whole can also replace the original system and be reused for many times during the design process, which can save much time further. To illustrate the MOR formulation procedure and moment matching we will take the first-order system’s model order reduction shown in (1) as an example and rewrite it as follows:

where , , . Here the matrices , , and are positive semidefinite matrices. The model order reduction based method can be written as follows.

Step 1. We build block Krylov subspace: .

Step 2. To implement Arnoldi process, specific algorithm is as follows:

(a)given matrix, , , and QR decomposition of matrix , let ;(b)compute , where ;(c)compute , where , ;(d)then compute the following parameter intrun: , where , .

Step 3. According to the transformation matrix, then get the original system of reduced order system:

where , , , .

PRIMA algorithm fully considered the system structure characteristics, so it keeps the passivity of the system and keeps the reduced order model system on the basis of structural characteristics, which has higher calculation accuracy.

5.2. Stochastic Reduced Order Model

To get stochastic reduced order model of MEMS, the first thing is to represent “randomness” in the original system matrices and the state vector; this is done by making use of polynomial chaos expansion. Suppose two input variables are assumed to be random and uniformly distributed. Considering a linear expansion in both of the input random variables, and using one-dimensional polynomial chaotic maps, the in (13) can be written as one-dimensional Hermite polynomial:

Assuming the input of a linear variable, there are

where is reduced order model coefficient matrix of deterministic model. are two-dimensional random space orthogonal polynomials. Once the coefficient matrix was determined by sparse grid method, the stochastic mapping reduced order model can be written as (11):

Let

Equation (18) may be cast in the following “deterministic” form:

Equation (19) is a deterministic system which can be easily used for time and frequency domain analyses. Then we can use deterministic model order reduction method to solve the stochastic reduced order model.

5.3. Computation of the Final Results: Expected Values and Standard Deviations

Once we have solved the systems numerically, that is to say, when we have got the matrices and , then we plug back the results into (11) and (12) to get an approximation of the probability density functions of the random eigenvalues and eigenvectors. As we have seen before this is the solution to the full problem with just one chosen mistuning parameter projected onto the subspace. Nevertheless we are not directly interested in this solution. Indeed, practically what is useful is the average value and standard deviation of each eigenvalue and eigenvector of the mistuned systems. But this is easy to get from the approximations of the probability density functions we have got. We directly use the linearity of the operator and the orthogonality of polynomial chaos on (16) to get

From (27) we get

Hence using the linearity and orthogonality properties as previously mentioned, we get

Finally using (28) and (29), we directly get the standard deviation:

We do exactly the same thing for eigenvectors and we directly get the expected value from the linearity and orthogonality properties applied on (28):

6. Numerical Studies: MEMS Examples

In order to illustrate the model reduction technique, we give an example of the well-known beam structure in a fluid environment. Figure 1 shows the beam structure; when a voltage is applied, the top plate of the structure bends downwards due to the resultant electrostatic force. Also when the beam bends, the pressure distribution of the ambient air under the beam increases. This pressure increase produces a backward pressure force that damps the beam motion. The beam structure has been used in many sensor applications.

Figure 1: Single-ended fixed electrostatic actuating microbeam structure.

The beam can be modeled by coupling the Euler beam equation with the electrostatic force and the Reynolds squeeze-film damping equation as follows:

Among them is the immunity of the direction, is young’s modulus, —moment of inertia’s is residual stress (Figure 5), is the density, is dielectric constant, is the air damping of the pressure, is the ambient pressure, is the air viscosity coefficient, is the distance from the plate under initial state, and is air mean free path (see Table 1).

Table 1: The parameters of single-ended fixed electrostatic actuating microbeam structure.

(beam length)
610

(beam width)
40

(beam thickness)
2.2

(initial gap distance)
2.3

(material Young’s modulus)
149 GPa

(material density)
2330 kg/m³

(air viscosity)
1.82 × 10⁻⁵ kg/(m·s)

(mean-free path of air)
0.064

Let , , substitute them into (3) and in and using Taylor series expansion, and take position, speed, and pressure as state variables like

we can get the micro beam for original system state equation:

To realize the stochastic reduced order model simulation, we firstly established the corresponding finite element model of the beam, to get the initial mathematical model such as rigidity matrix and mass matrix, which can be used as random input parameters, do standard KL series expansion with these matrix and do polynomial chaos expansion with output matrix; after expansion finish, we use Smolyak method to decide which point should be used as reduced order model simulation; then we use SPRMI method at the Smolyak point; at last, we can get the probability distribution from order reduction model simulation results. In MEMS, because micromechanical usually driven by electrostatic force, the machine is usually a conductor in static electric field, subject to the effect of electrostatic force deformation occurring, and the structure size and the location of the components will influence the distribution of electric field in turn, so it is a kind of strong coupling. Here we adopt finite element method combined with the boundary element method to solve the force electric coupling problem of MEMS, solving process as follows.

Step 1. Under the condition of the plate without displacement, giving an initial voltage on the plate and solving problems of electrostatic field with the boundary element, we can get the distribution of electrostatic force of plate.

Step 2. Solve the displacement of plate under electrostatic force through the finite element method; if it reaches precision displacement, skip to Step 4; otherwise, go to Step 3.

Step 3. Based on the displacement of the plate to form the boundary of the electrostatic field, use the boundary element method to solve the problem of electrostatic field, obtain the distribution plate of the electrostatic force, and let ; go to Step 2.

Step 4. The process is over.

We use mesh as shown in Figure 2, where represents the number of inner grid points in the direction and is the number of inner grid points in the direction, according to the two-dimensional model establishment of beam parameters as shown in Figure 2. The beam was divided with 2D quadrilateral element PLANE82, and the air medium was divided with PLANE121; displacement constraints were applied to fix the microbeam anchors in ANSYS (Figure 4); each degree of freedom beam base and left side of beam were restricted; the finite element model is shown in Figure 3. Then we applied effect of voltage on the beam structure. After the mesh is generated, we then project the unknowns and onto the mesh points and apply the trapezoidal rule to discretize the special integral operator and the central difference method to discretize the special derivative operators as follows:

Accelerometer beam clamped at one end, so the boundary conditions are ; substitute (36) into (35) to discrete; then map the state vector to the various nodes; we can convert discrete equations into dimensional state space equation (9), where ,

where , , and are the microdisplacement of the beam end and is the input voltage of the system.

Figure 2: Structure model of microbeam.

Figure 3: Finite element model of microbeam.

Figure 4: The displacement of microbeam.

Figure 5: The stress of microbeam.

Applying the boundary conditions, the finite element model is solved and can change the endpoint microbeam under different voltage conditions. We get the deflection of the beam’s end point as shown in Figure 6. The voltage load vector along with the mass and stiffness matrices were extracted as input random parameters.

Figure 6: The endpoint deflection of the beam.

A generalized polynomial chaos expansion is now generated to represent the dependence of the maximum displacement of the membrane as a function of the input parameter space. This includes geometric parameters thickness , length , width , Young’s modulus , and the initial gap between the membrane and the pull-down electrode. The variability in the thickness and gap is estimated from known processing conditions. The parameter limits are selected based on the requirement; the device cannot experience pull-in instability over the time span of interest. The parameter limits and distribution types are presented in Table 2.

Table 2: The GPC response inputs parameters.

Variable
Distribution
Range
Units

1.6–3.6
μm

549–671
μm

27–33
μm

177–216
GPa

4.0–5.0
μm

Then making use of Smolyak algorithm consisting of only 29 grid points in the two-dimensional random parameter spaces for numerical integration, the grid result is shown in Figure 7.

Figure 7: Smolyak grid (total number of points 29).

Then, we compute 29 times of model reduction simulation at Smolyak point; the results were shown in Figure 8. From the figure, we learn when time ≤ 0.7 s. The most simulation results are almost the same. After that, we get kinds of deflection.

Figure 8: Model reduction simulation results.

We also compute that the resulting probability distribution of the beam maximum deflection is presented in Figures 9, 10, and 11. In Figure 9, we just show geometric parameters thickness , length , and width as random parameters, and Figure 10 shows Young’s module as random parameter, and Figure 11 shows the initial gap distance as random parameters.

Figure 9: The geometric .

Figure 10: Young’s modulus.

Figure 11: The initial gap distance.

7. Conclusions

We have developed a new stochastic model order reduction frame to MEMS model uncertainty and application in single fixed micro-beam. In the computing process, firstly, we construct finite model, then extract the matrix of analysis, get the random input matrix, and use KL series expansion method for input parameters matrix; then, we express the random output variables as standard random variable polynomials by PC method; while obtaining the polynomials of output variables, we use Smolyak method to compute coefficients of PC equation, then take the reduced matrix to PRSIM model reduction method, and compute the statistics such as the mean and standard deviation of the desired system response and deflection of microbeam. We can get that all the computational effort is focused on the estimation of a small set of PC coefficients in the proposed methodology. This leads to a considerable gain compared to the techniques based on Monte Carlo simulation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The paper was supported by the National Natural Science Foundation of China “MEMS Multifields Uniform Simulation Model Construction and Optimism Research” (61100101); the paper was also supported by Nature Science Foundation of Zhejiang Province “MEMS united constraints model order reduction method and solution research (LY14E050026)” and Zhejiang Province Higher School Visiting Scholar Teacher Professional Development Project: MEMS Model order reduction technique and application.

References

Y. Gong, H. Chen, H. Liu, Z. Li, and G. Chen, “MEMS multi-fields uniform model simulation and optimization method research and application,” Sensors and Transducers, vol. 156, no. 9, pp. 237–246, 2013.

M. I. Younis, E. M. Abdel-Rahman, and A. Nayfeh, “A reduced-order model for electrically actuated microbeam-based MEMS,” Journal of Microelectromechanical Systems, vol. 12, no. 5, pp. 672–680, 2003.

A. H. Nayfeh, M. I. Younis, and E. M. Abdel-Rahman, “Reduced-order models for MEMS applications,” Nonlinear Dynamics, vol. 41, no. 1-3, pp. 211–236, 2005.

M. Liu, K. Maute, and D. M. Frangopol, “Multi-objective design optimization of electrostatically actuated microbeam resonators with and without parameter uncertainty,” Reliability Engineering and System Safety, vol. 92, no. 10, pp. 1333–1343, 2007.

S. K. De and N. R. Aluru, “Coupling of hierarchical fluid models with electrostatic and mechanical models for the dynamic analysis of MEMS,” Journal of Micromechanics and Microengineering, vol. 16, no. 8, pp. 1705–1719, 2006.

G. Li and N. R. Aluru, “Efficient mixed-domain analysis of electrostatic MEMS,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 22, no. 9, pp. 1228–1242, 2003.

P. S. Sumant, N. R. Aluru, and A. C. Cangellaris, “A methodology for fast finite element modeling of electrostatically actuated MEMS,” International Journal for Numerical Methods in Engineering, vol. 77, no. 13, pp. 1789–1808, 2009.

G. S. Fishman, Monte Carlo: Concepts, Algorithms, and Applications, Springer Series in Operations Research, Springer, Berlin, Germany, 1996.

W.-L. Loh, “On latin hypercube sampling,” The Annals of Statistics, vol. 24, no. 5, pp. 2058–2080, 1996.

M. Stein, “Large sample properties of simulations using Latin hypercube sampling,” Technometrics, vol. 29, no. 2, pp. 143–151, 1987.

W. Gilks, S. Richardson, and D. Spiergelhalter, Markov Chain Monte Carlo in Practice, Chapman & Hall, London, UK, 1995.

J. S. Kong, D. M. Frangopol, M. Raulli et al., “A methodology for analyzing the variability in the performance of a MEMS actuator made from a novel ceramic,” Sensors and Actuators A: Physical, vol. 116, no. 2, pp. 336–344, 2004.

J. S. Han and B. M. Kwak, “Robust optimal design of a vibratory microgyroscope considering fabrication errors,” Journal of Micromechanics and Microengineering, vol. 11, no. 6, pp. 662–671, 2001.

R. Liu, B. Paden, and K. Turner, “MEMS resonators that are robust to process-induced feature width variations,” Journal of Microelectromechanical Systems, vol. 11, no. 5, pp. 505–551, 2002.

M. Allen, M. Raulli, K. Maute, and D. M. Frangopol, “Reliability-based analysis and design optimization of electrostatically actuated MEMS,” Computers & Structures, vol. 82, no. 13-14, pp. 1007–1020, 2004.

J. Lienemann, E. B. Rudnyi, and J. G. Korvink, “MST MEMS model order reduction: requirements and benchmarks,” Linear Algebra and Its Applications, vol. 415, no. 2-3, pp. 469–498, 2006.

M. Liu, K. Maute, and D. M. Frangopol, “Multi-objective design optimization of electrostatically actuated microbeam resonators with and without parameter uncertainty,” Reliability Engineering and System Safety, vol. 92, no. 10, pp. 1333–1343, 2007.

P. Heydari and M. Pedram, “Model reduction of variable-geometry interconnects using variational spectrally-weighted balanced truncation,” in Proceedings of the IEEE/ACM International Conference on Computer-Aided Design (ICCAD '01), pp. 586–591, San Jose, Calif, USA, November 2001.

Y. Liu, L. T. Pileggi, and A. J. Strojwas, “Model order-reduction of RCL interconnect including variational analysis,” in Proceedings of the ACM/IEEE 36th Annual Design Automation Conference (DAC '99), pp. 201–206, June 1999.

X. M. Zhang, F. S. Chau, C. Quan, Y. L. Lam, and A. Q. Liu, “A study of the static characteristics of a torsional micromirror,” Sensors and Actuators, A: Physical, vol. 90, no. 1-2, pp. 73–81, 2001.

L. Daniel, O. C. Siong, L. S. Chay, K. H. Lee, and J. White, “A multiparameter moment-matching model-reduction approach for generating geometrically parameterized interconnect performance models,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 23, no. 5, pp. 678–693, 2004.

P. Li, F. Liu, X. Li, L. T. Pileggi, and S. R. Nassif, “Modeling interconnect variability using efficient parametric model order reduction,” in Proceedings of the Design, Automation and Test in Europe (DATE '05), pp. 958–963, Munich, Germany, March 2005.

Y. T. Li, Z. Bai, Y. Su, and X. Zeng, “Parameterized model order reduction via a two-directional arnoldi process,” in Proceedings of the IEEE/ACM International Conference on Computer-Aided Design (ICCAD '07), pp. 868–873, November 2007.

B. Bond and L. Daniel, “Parameterized model order reduction of nonlinear dynamical systems,” in Proceedings of the IEEE/ACM International Conference on Computer-Aided Design (ICCAD '05), pp. 486–493, November 2005.

B. N. Bond and L. Daniel, “A piecewise-linear moment-matching approach to parameterized model-order reduction for highly nonlinear systems,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 26, no. 12, pp. 2116–2129, 2007.