Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 5014235, 20 pages
https://doi.org/10.1155/2017/5014235

Research Article

A Two-Step Global Maximum Error Controller-Based TPWL MOR with POD Basis Vectors and Its Applications to MEMS

Ying Liu¹ and Xiaodong Wang2,3

¹School of Arts and Sciences, Shaanxi University of Science & Technology, Xi’an 710021, China
²School of Mathematical Sciences, Peking University, Beijing 100871, China
³Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, China

Correspondence should be addressed to Ying Liu; liu_ying@sust.edu.cn and Xiaodong Wang; xiaodongwang@nwpu.edu.cn

Received 25 May 2017; Accepted 14 September 2017; Published 22 October 2017

Academic Editor: Alessandro Lo Schiavo

Copyright © 2017 Ying Liu and Xiaodong Wang. 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.

In our previous study, we have proposed a linearization point (LP) selection method based on a global maximum error controller for the trajectory piecewise-linear (TPWL) method. It has been demonstrated that this method has many advantages over other existing methods. In this paper, a more efficient version of this method is presented, which introduces a preliminary LP selection procedure and constructs projection matrix by the proper orthogonal decomposition (POD) method. Compared with the original method, the improved method takes much less time for extracting a reduced-order model (ROM) of similar quality and gets some other benefits (such as being easier to implement, having lower memory requirement, and enhanced flexibility). The effectiveness of the new method is fully demonstrated by a diode transmission line RLC circuit. And then, the method is applied to three more complicated microelectromechanical systems (MEMS) devices, which are a micromachined switch, an electrostatic micropump diaphragm, and a thermomechanical in-plane microactuator.

1. Introduction

Originally, model order reduction (MOR) was developed in the area of systems and control theory, which studies properties of dynamical systems by reducing the complexity and preserving the input-output behaviors as much as possible [1]. Over the past decades, the idea of MOR has been widely concerned in both the academic and engineering fields. It has been a hot area of research in the fields including but not limited to integrated circuit design, systems and control, fluid dynamics, computational biology, microelectromechanical systems (MEMS) [2, 3], and computational mathematics.

After extensive study, the theories and applications about linear MOR methods have been well-developed. By contrast, the nonlinear MOR methods are still in their development stage [4–6], although a great deal of effort has been devoted. Owing to good properties of linear problems and maturity of linear MOR methods, the most obvious approach in dealing with nonlinear problems is linearization. However, the linearized model, in general, has good approximation only at the states which are very close to the linearization point (LP). In order to obtain an overall good approximation, a very promising nonlinear MOR method, called trajectory piecewise-linear (TPWL) method, was devised by Rewienski and White in 2001 [7]. The idea of TPWL method is selecting multiple LPs along a training trajectory and approximating the nonlinear model piecewise-linearly [7, 8]. The advantage of this method is the ability to reduce both dimension of the full nonlinear system and cost for evaluating it.

At present, the TPWL method has received much attention due to its good performance in generating high quality reduced-order models (ROMs) [7–21]. But the development of this method is imperfect, and there are still several issues to be resolved. The most serious issues include how to select high quality LPs, which method to choose in projection, and how to construct better weighting functions. Because a reasonable LP selection method has decisive influence on the quality and computational complexity of a final ROM, much effort has been put on this topic [7, 10, 16–20]. By studying the existing LP selection methods, we have proposed another LP selection method in our previous study [14]. The new method is based on a global maximum error controller (GMEC) and the corresponding TPWL method is known as TPWL-GMEC method. Compared with other existing methods, the TPWL-GMEC method can select LPs of higher quality. The obtained ROM has smaller size and higher accuracy than those generated by using other LP selection methods. Besides, this method can improve the model expandability by generating one TPWL model for multiple training inputs. To some extent, the LPs selected by TPWL-GMEC method can be considered to be optimal, since the corresponding ROM always minimizes the maximum error. But as we have pointed out [14], a shortcoming of the TPWL-GMEC method is that it takes relatively long time to extract the ROMs.

In this paper, a more efficient version of the TPWL-GMEC method is presented by introducing two improvements. On the one hand, for obtaining the maximum error state, the TPWL-GMEC method needs to compute the errors at all the states included in the state-set. So, if we can reduce the size of state-set before LPs selection without losing accuracy, the computational cost will be reduced. On the other hand, for updating the ROM after a new LP is selected, the TPWL-GMEC method needs to rebuild the global projection matrix. So, if a predefined global projection matrix is used, the basis vector updates will be avoided. For the two reasons, the improved algorithm introduces a preliminary LP selection procedure and constructs projection matrix by the proper orthogonal decomposition (POD) method [22]. Compared with the original method, the improved method takes less time to extract a ROM of similar quality and gets some other benefits. Numerical test verifies the effectiveness of the new method. And after the test, the method is applied to the MOR of several MEMS devices.

The rest of this paper is organized as follows. In the next section, the foundation of the TPWL method will be briefly reviewed. In Section 3, the basic idea and procedure of our improved algorithm is given. Section 4 shows some simulation results for demonstrating the effectiveness of the new method and Section 5 gives the applications of the method to the MOR of some MEMS devices. Finally, we draw our conclusions in Section 6.

2. Foundation of the TPWL-GMEC Method

The basic idea behind the TPWL method is to select multiple LPs and linearize the nonlinear model near each LP. Then, each local linearized model is reduced by projecting it into a global subspace, and a ROM is constructed as a weighted sum of all local linearized ROMs [8]. In order to describe it clearly, consider the nonlinear system of the following form:

where represents the state vector, is a nonlinear vector function, distributes the input , and maps to the system response .

Let be the state-set associated with a certain training input for , where are approximate values of at time instants . The LPs are assumed having been selected out according to some rules, where is the number of LPs. Through linearizing at those LPs, a set of the local linearized models give

where is the Jacobian of evaluated at . Next, the local projection matrices for each local linearized model are generated by using a linear MOR method, such as Krylov-subspace methods [7, 19] and truncated balanced realization (TBR) method [23]. The global projection matrix is constructed by merging and orthogonalizing . According to , the local linearized models are reduced to size ; namely,

where , , , , and . Finally, a TPWL ROM can be represented as a weighted combination of all local linearized ROMs

where is the weighting function. For any , it satisfies the conditions and . The most used weighting functions are as follows:

where , , and is a positive constant to control the compactness of weighting functions and usually taken as .

Through the above description, the LP selection is an essential part of the whole MOR process. The simulation time of the TPWL ROM is very short, so the responses of the current ROM model and the full nonlinear model at each time step can be easily compared. Because of this fact, in our previous study [14], the state with the global maximum error is selected as a new LP in each cycle. The corresponding TPWL-GMEC method is shown in Algorithm 1. Numerical tests have demonstrated that Algorithm 1 can select high quality LPs. The corresponding TPWL models usually have higher accuracy and smaller sizes than those generated by other methods. Besides, this method can improve the model expandability by generating one TPWL model which contains the information from multiple training trajectories. However, it needs relatively long time to extract the ROMs, so a more efficient version needs to be presented.

(1) Given a training input ;

(2) Simulate the nonlinear system (1) and store the states at all time steps

;

(3) Choose the initial state as the first LP , and set ;

(4) Update and built the current TPWL model;

(5) Simulate the current TPWL model to get the responses ;

(6) Find out the maximum error state :

Set and , where the root mean square error

(RMSE) is defined by

for  

Compute ;

if    then  , ;

end for

(7) if   and has not been selected as a LP

then  select the LP , set ;

(8) if   is increased in step () and ( is the maximum allowable LP number)

then go to step ().

(9) Return the selected LPs and end the loop.

Algorithm 1: Original TPWL-GMEC method.

3. A More Efficient LP Selection Algorithm

In this section, a more efficient version of the TPWL-GMEC method is given by introducing the following improvements.

3.1. Two-Step LP Selection Strategy

Recall that the original TPWL-GMEC method selects LPs from the whole state-set. The errors at all the states included in the state-set need to be computed for selecting a new LP, and this process brings much computational cost. Therefore, the first strategy to improve the efficiency of the TPWL-GMEC method is through reducing the size of state-set before LPs selection. That is, a preliminary LP selection procedure is added to the original algorithm.

First, according to comparing the distance between states, a set of LP candidates meeting the accuracy requirement would be selected. Specifically, it estimates the states along a training trajectory one by one and selects the current state as a new LP candidate once the distance between the current state and the previous LP candidates exceeds a predefined threshold . In fact, this approach is just the LP selection method used by the founders of TPWL method [8], and the only difference is that we use smaller to meet higher accuracy requirement in this work. It should be noted that the LP candidates are not suitable for directly generating a ROM, because their quantity is often large and many redundant states are included. Hence next, the LP candidates are optimized to obtain the final LPs by using the GMEC based selection algorithm. Because the preliminary selection procedure is very efficient and the LP candidates are much less than the whole states, the two-step LP selection strategy is expected to greatly improve the efficiency of the TPWL-GMEC method.

In the above process, the threshold needs to be set in advance. In theory, it can be set as any positive real number. But an improper may bring undesirable effects on the efficiency of TPWL-GMEC method and the accuracy of the finial ROM. In applications, it is usually enough to calculate reasonable by using the information about the initial and final states. That is,

where is a positive number, and are, respectively, the initial and final state vectors. Generally, for a given nonlinear system, a higher accuracy requirement usually needs smaller (i.e., bigger ).

3.2. Combining with POD Method

In the original TPWL-GMEC method, the local projection matrices for each local linearized model are generated by using the Krylov-subspace method, and the global projection matrix is constructed by merging and orthogonalizing the local projection matrices. Once a new LP is selected and a new local linearized model is added, the global projection matrix must be updated to rebuild the ROM. This process brings much computational cost. By contrast, when a set of states is used to directly generate the global projection matrix by POD method [25], the global projection matrix is ready before LP selection and does not need to be repeatedly updated. Therefore, the second strategy to improve the efficiency of the TPWL-GMEC method is through combining with the POD method.

In addition to improving efficiency, combining with POD method can also bring other benefits. First, when the Krylov-subspace method is used, only after all the LPs are selected, the final global projection matrix can be constructed. So the matrices associated with the LPs, including mass matrices and Jacobi matrices, cannot be reduced immediately after their generation. Consequently, a lot of memory is allocated to temporarily store these matrices such that the simulation scale of the TPWL-GMEC method is limited. However, after combining with POD method, the global projection matrix is constructed before LPs selection. So, the matrices can be reduced immediately once they are assembled, and the memory requirement will be greatly reduced. Second, because the POD bases are merely related to state vectors and it is not constrained by the specific form of the system equations, the application range of the corresponding TPWL-GMEC method would be greatly increased.

3.3. The Improved Algorithm

By using any of the above improvements, the efficiency of the TPWL-GMEC method would be expected to be improved. In order to increase the efficiency significantly, an improved algorithm containing two improvements is presented. Algorithm 2 shows the detailed procedure of the improved TPWL-GMEC method.

(preliminary LP selection)

(1) Given a training input and a threshold ;

(2) Choose the initial state as the first LP candidate , and set ;

(3) Simulate the nonlinear system (1) to obtain the state-set and select LP

candidates:

for  

Simulate the nonlinear system (1) to get the state and record it;

if  ,  then select as the th LP candidate , record

the corresponding time instant , and set ;

end for

(secondary LP selection)

(4) Given the maximum error limit ;

(5) Select a set of states from the state-set of the nonlinear system;

(6) Construct the POD projection matrix through applying the SVD algorithm

to the set of states, that is, ;

(7) Choose the first LP candidate as the first LP , and set ;

(8) Built the current TPWL model;

(9) Simulate the current TPWL model to get the responses ;

(10) Find out the maximum error state :

Set and

for  

Compute ;

if    then  , ;

end for

(11) if   and has not been selected as a LP, then  select the

LP , set , and go to ();

(12) if    then return the selected LPs and end the loop.

Algorithm 2: Improved TPWL-GMEC method.

It seems that Algorithm 2 is more complex than the original one (Algorithm 1). But the improved TPWL-GMEC method is actually simpler than Algorithm 1 in program implementation. Particularly, when is set as zero, Algorithm 2 only contains the second improvement.

4. Numerical Validation

For validating the improved TPWL-GMEC method, we present here a diode transmission line RLC circuit. As shown in Figure 1, it consists of resistors, inductors, capacitors, and diodes, where the resistance , inductance , capacitance , and the constitutive equation of diodes . The current source entering into node 1 is the input , and the voltage at node 1 is the output response . The system equations of the nonlinear transmission line RLC circuit can be represented as

Figure 1: Schematic of a nonlinear transmission line RLC circuit.

For this example, we use single- and dual-sinusoidal inputs to validate our method. The order of the full nonlinear model, the simulation time, and the time step size are set as , , and , respectively. Other user-defined parameters except the threshold are the same as we used in [14].

For ease of comparison, we denote the TPWL-GMEC method without any improvement as TPWL-GMEC (Krylov), with the two-step LP selection strategy as TPWL-GMEC2 (Krylov), with the POD basis as TPWL-GMEC (POD), and with two improvements as TPWL-GMEC2 (POD). Figure 2 shows the simulation results of the ROMs generated by the TPWL-GMEC2 (POD) method with . As can be seen, the difference between the simulation results of the ROMs and the full nonlinear model becomes invisible as we use a small enough . This behavior is very similar to that of the TPWL-GMEC (Krylov) method demonstrated in [14]. Further, the error comparisons shown in Figure 3 indicate that the TPWL-GMEC2 (POD) method can inherit the good performance of the TPWL-GMEC (Krylov) method.

(a)

(b)

Figure 2: Transient results of the TPWL-GMEC2 (POD) models for the diode transmission line RLC circuit with different : (a) the single-sinusoidal input ; (b) the dual-sinusoidal input .

(a)

(b)

Figure 3: Root mean square error comparisons of TPWL-GMEC2 (POD) ROMs and TPWL-GMEC (Krylov) ROMs for the diode transmission line RLC circuit: (a) the single-sinusoidal input ; (b) the dual-sinusoidal input .

For a more detailed comparison, Tables 1 and 2 show the order of ROMs , the LP number , the actual maximum error , the time cost for extracting the ROMs , and the simulation time of the ROMs of four TPWL-GMEC methods with different for the single- and dual-sinusoidal inputs, respectively. From two tables, one easily finds that four methods undoubtedly meet the accuracy requirement for all the cases, and the LP numbers and the actual maximum errors of the four methods for every case do not appear to be much different from each other. As to the order and simulation time of the ROMs, four methods show two distinctive performances. That is, the methods with POD basis have the order of ROMs independent of and often possess smaller and less for greater . Most importantly, the TPWL-GMEC2 (POD) method shows the highest efficiency for almost all cases, which is followed in descending order by TPWL-GMEC2 (Krylov) method, TPWL-GMEC (POD) method, and TPWL-GMEC (Krylov) method. For small , the efficiency of the TPWL-GMEC2 (POD) method is even several times higher than the original TPWL-GMEC (Krylov) method. Table 3 shows the results obtained by the standard TPWL method proposed by the Rewienski and White [7]. By comparison, the ROMs obtained by the TPWL-GMEC2 (POD) method always have better accuracy than those obtained by the standard TPWL method with the same LP number. Furthermore, in some cases the efficiency of the TPWL-GMEC2 (POD) method is even higher. In summary, the improved TPWL-GMEC method can extract ROM with higher quality at low time cost, and thus it makes up for the shortcoming of the original TPWL-GMEC method.

Table 1: Comparisons of the four methods for the input .

Method

TPWL-GMEC (Krylov)
3
1

1.528
0.05

TPWL-GMEC2 (Krylov)
3
1

0.681
0.049

TPWL-GMEC (POD)
10
1

1.894
0.073

TPWL-GMEC2 (POD)
10
1

0.856
0.072

TPWL-GMEC (Krylov)
18
5

9.802
0.548

TPWL-GMEC2 (Krylov)
18
5

4.324
0.551

TPWL-GMEC (POD)
10
5

6.8
0.157

TPWL-GMEC2 (POD)
10
5

1.602
0.158

TPWL-GMEC (Krylov)
25
12

37.266
1.52

TPWL-GMEC2 (Krylov)
26
13

23.883
1.706

TPWL-GMEC (POD)
10
11

23.146
0.811

TPWL-GMEC2 (POD)
10
12

9.291
0.855

Table 2: Comparisons of the four methods for the input .

Method

TPWL-GMEC (Krylov)
3
1

1.529
0.05

TPWL-GMEC2 (Krylov)
3
1

0.75
0.05

TPWL-GMEC (POD)
16
1

2.21
0.172

TPWL-GMEC2 (POD)
16
1

0.975
0.172

TPWL-GMEC (Krylov)
21
5

10.331
0.707

TPWL-GMEC2 (Krylov)
20
5

4.818
0.651

TPWL-GMEC (POD)
16
3

5.391
0.267

TPWL-GMEC2 (POD)
16
3

1.729
0.265

TPWL-GMEC (Krylov)
32
15

63.062
2.832

TPWL-GMEC2 (Krylov)
33
14

37.775
2.848

TPWL-GMEC (POD)
16
17

32.852
0.963

TPWL-GMEC2 (POD)
16
13

8.21
0.741

Table 3: Results obtained by the original TPWL method proposed by Rewienski and White.

Input

3
1

0.659
0.051

17
5

3.334
0.544

29
12

11.485
2.187

3
1

0.665
0.053

12
3

3.094
0.246

32
13

14.359
2.610

Finally, for giving advice to choose , Figures 4 and 5 show the effects of on the performance of two methods with the two-step LP selection strategy. When decreases, the actual maximum error, the LP number, and the ROM order of the two methods quickly tend to those of the corresponding methods without the two-step LP selection strategy. And to avoid the loss of precision, should be less than . More importantly, for any , the time cost of extracting ROMs by the TPWL-GMEC2 (POD) method is the least, but as decreases the time cost will be gradually increased. According to many simulations (not only shown here), we find that if is between 5 and 100, then the efficiency of the TPWL-GMEC2 (POD) method would be increased by several times. So, the suggested value of is between and .

(a)

(b)

(c)

(d)

Figure 4: Effects of on the performance of the TPWL-GMEC methods for the input and : (a) maximum error; (b) order of ROM; (c) number of LPs; (d) extracting time.

(a)

(b)

(c)

(d)

Figure 5: Effects of on the performance of the TPWL-GMEC methods for the input and : (a) maximum error; (b) order of ROM; (c) number of LPs; (d) extracting time.

5. Applications to MEMS

In this section, the two-step TPWL-GMEC method with POD basis (i.e., TPWL-GMEC2 (POD) method) is applied to MOR of several MEMS devices to further show its performance. The examples presented here include a micromachined switch, an electrostatic micropump diaphragm, and a thermomechanical in-plane microactuator. Without otherwise specified, the simulations in this section will be done by setting according to the above general rule of choosing . Here below are the details of applications.

5.1. Example  1: Micromachined Switch

The first example is a micromachined switch as illustrated in Figure 6 [19, 26]. Due to its complex structure, strong nonlinear pull-in phenomenon, and large computation amount, the performance of a MOR method can be thoroughly tested. Usually this device is regarded as a very challenging example for MOR methods and helps to find out their potential shortcomings.

Figure 6: Schematic illustration of the micromachined switch [19].

As shown in Figure 6, the micromachined switch consists of a polysilicon beam suspended over a silicon substrate. The beam has length  μm, width  μm, and height  μm. When voltage is applied between the beam and the substrate, the beam deflects toward the silicon substrate. It may be reasonable to assume that the deflection of the beam is uniform across its width since the length is much greater than the width. Under this assumption, the dynamical behavior of this coupled electromechanical-fluid system can be governed by the following 1D Euler’s beam equation and 2D Reynolds squeeze film damping equation [27]:

where is Young’s modulus, the moment of inertia of the beam, the stress coefficient, the electrostatic force, the density, the ambient pressure, the air viscosity, the Knudsen number, the height of the beam above the substrate. and the pressure distribution in the air below the beam. is approximated by

with the permittivity of vacuum and the applied voltage.

In order to facilitate the MOR procedure, (8) are transformed by defining new variables

Using the above variables, (8) may be rewritten as

For (11), we consider the following initial conditions

and the following boundary conditions:

where  μm is the length of the gap between the undeflected beam and the substrate.

In simulations, the parameters  GPa,  MPa,  kg/m³,  Pa,  kg/(m·s), and  μm are adopted. Equations (11) are discretized by a central finite difference scheme on a gird, which leads to a full nonlinear system of order . A compact form of the system can be written as

The time step size for solving system (14) is set as  μs, and the center height of the beam is selected as the output.

Figure 7 shows the results obtained by TPWL-GMEC2 (POD) method for different . Two inputs considered in the simulations are ( for and for ) and with  kHz. It can be seen that the results of the ROMs for both inputs converge to those of the full nonlinear model very fast as decreases. When , the accuracy of the ROMs is satisfactory and the sizes of the ROMs (measured by and ) are quite small. This preliminarily shows the good performance of the TPWL-GMEC2 (POD) method for the micromachined switch. Incidentally, for the case , it is not easy for the original TPWL method to select out a set of LPs with high quality since the system response is periodic and many redundant LPs could be selected. However, such difficulty does not appear in the current TPWL-GMEC2 (POD) method because any redundant LP would be removed.

(a)

(b)

Figure 7: Results obtained by TPWL-GMEC2 (POD) method for different (example 1): (a) input voltage ; (b) input voltage with  kHz.

In order to verify whether the TPWL-GMEC2 (POD) method retains the good expandability of the original TPWL-GMEC method, a ROM is generated for the training input when , and then it is used to simulate other testing inputs. Figure 8 shows the corresponding results. Firstly, from Figure 8(a), we see that the ROM generated for the nonperiodic step input has good approximation to the periodic cosine input . That is to say, the ROM generated for the step input can totally contain the response information of the corresponding periodic inputs. Secondly, in Figure 8(b), the ROM shows very good approximation for the inputs and , small deviation for the input , and large deviation for the input . Notice that the pull-in phenomenon occurs when and does not when . Therefore, the ROM for does not fully contain the response information of . When we use this ROM to simulate , the simulation result certainly becomes unsatisfactory. Nevertheless, the results in Figure 8 still show the good expandability of the TPWL-GMEC2 (POD) method.

(a)

(b)

Figure 8: Results of a ROM generated for the training input with the testing inputs different from the training input: (a) input voltage with  kHz; (b) input voltages .

Since the input produces distinctive response behavior with other mentioned inputs, next we will try to generate a satisfactory ROM for this input by using the TPWL-GMEC2 (POD) method. If the default value of (i.e., ) is used, a satisfactory approximation is not obtained even when (see Figure 9(a)). In this case, 10 LP candidates are selected by the preliminary selection procedure. But even if all LP candidates are selected, the ROM still cannot meet the accuracy requirement . So, if we want to get good approximation, should be decreased. The result in Figure 9(b) shows great improvement when is used. However, it should be noted that smaller would introduce more LPs and many of them would be gathered at the end of the trajectory. So, it is not economic to further improve the approximation for this case.

(a)

(b)

Figure 9: Results obtained by TPWL-GMEC2 (POD) method for with different : (a) ; (b) .

5.2. Example  2: Electrostatic Micropump Diaphragm

Micropumps have played an important role in biological fluid handling, drug delivery, and thermal management of electronic components [24]. At present, the actuation mechanisms of micropumps include piezoelectric, electrostatic, pneumatic, magnetic, electromagnetic, and thermal expansion. In this subsection, we focus on a diaphragm displacement micropump actuated by an electrostatic mechanism. As illustrated in Figure 10(a), the upper plate is a rigid counter-electrode, the middle is a flexible diaphragm, and the bottom is a pumping chamber. During the operation of the pump, fluid flows into the chamber from the left side and flows out from the right side. Two valves are located at the inlet and outlet, respectively. When a voltage is applied on the pump, the motion of the diaphragm causes pressure difference, and hence the fluid in the reservoir is forced to flow in the microchannel.

(a)

(b)

Figure 10: Electrostatic micropump diaphragm: (a) device layout; (b) schematic illustration of the diaphragm [24].

Because the analysis of the micropump is mainly focused on the dynamic response of the diaphragm under prescribed loads, the micropump in simulations may be simplified as a model shown in Figure 10(b). The upper flexible electrode (i.e., the diaphragm) and the lower ground electrode are both circular with radius . The flexible electrode is clamped all along the boundary and separated from the ground electrode by a uniform gap . The two electrodes are assumed to be homogeneous and perfect conductors, and the space between them is filled with a homogeneous dielectric medium with permittivity .

If the thickness of the diaphragm is far less than the radius such that it can be modeled as a thin plate, the displacement field near the equilibrium state would be calculated by the Germaine-Lagrange equation

coupled with suitable boundary conditions. Here, is the plate stiffness where is Young’s modulus and is Poisson’s ratio of the material. And is the double Laplacian operator in polar space coordinates; that is,

The electrostatic pressure under an applied potential may be written as

Considering inertial effects and damping, (15) can be rewritten as

where is the density, the quality factor (representing the damping), and the natural frequency of the structure. By defining dimensionless variables , , and , (18) can be written in the dimensionless form (the superscript has been omitted)

where the mass coefficient (involving geometric and material parameters) and load parameter (involving electrostatic actuation parameters) are defined by

For (19), we consider the following initial and boundary conditions:

for and .

In simulations, (19) is discretized by standard finite difference scheme with equidistantly distributed grid points, and the corresponding compact form of the spatial discretized system can be written as

where is the state vector, is the nonlinear input vector, , , and are mass matrix, damping matrix, and stiffness matrix, respectively. The time step size for solving system (22) is set as , and the center point deflection of the diaphragm is selected as the output. Among three dimensionless model parameters, the mass coefficient is fixed at , and the dynamic response behavior of the diaphragm is analyzed by changing the values of and .

Obviously, (22) is a second-order ODE system. When the original TPWL method and the original TPWL-GMEC method are applied, we need to adopt second-order Krylov-subspace method or transform second-order system to first-order system [28, 29]. However, for the TPWL-GMEC2 (POD) method, the MOR procedures for second-order system and first-order system are nearly the same. The ROM of (22) generated by the TPWL-GMEC2(POD) method has the following form:

where , , , , , , and is the POD projection matrix and is the Jacobi of at .

Figure 11 shows the results obtained by TPWL-GMEC2 (POD) method for different . When , obvious difference can be found between the results of the ROMs and full nonlinear model. But when is decreased to , the results are almost consistent. It should be mentioned that the accuracy of the ROM is hard to be further improved by only increasing the LP number while maintaining the projection matrix. For the cases considered in Figure 11, when , even if hundreds of LPs are selected, the ROM still cannot meet the accurate requirement. So at this situation, to further improve the accuracy of the ROM, the order of the projection matrix must be increased. By doing this, the size of the ROM would be greatly increased, so we should make a tradeoff between the size and accuracy of the ROM. In the following simulations of this example, the accurate requirement will be all set as .

(a)

(b)

Figure 11: Results obtained by TPWL-GMEC2 (POD) method for different (example 2): (a) and ; (b) and .

In Figures 11(a) and 11(b), two distinct response behaviors of the diaphragm can be observed. For , the diaphragm exhibits periodic fluctuation, and the amplitude is decreasing as the elapse of time. It can be expected that the diaphragm will ultimately reach an equilibrium state after enough long time. But for , the dimensionless deflection of the diaphragm reaches the maximum value 1 very fast. That is, the diaphragm snaps together with the ground electrode. In order to ascertain the critical of the pull-in effect, a lot of simulations are performed and Table 4 gives some results. From the table, we find that the critical shows a slow decline as increases and eventually trends to some constant.

Table 4: The critical for the pull-in effect of the diaphragm for different .

1
5
10
20
50
100
1000

15.56
13.71
13.38
13.21
13.11
13.07
13.04

Figure 12 shows the influence of and on the dynamic response behavior of the diaphragm. As can be seen, the results of the ROMs show very good agreement with those of the full nonlinear model. From Figure 12(a), it can be found that the diaphragm costs longer time to attain its equilibrium states for greater , and the deflection of the diaphragm increases as increases. From Figure 12(b), greater always results in larger fluctuation amplitude of the diaphragm and longer time for the diaphragm to attain its equilibrium state, where the equilibrium state seems independent of . For quantitative analysis, Figure 13 shows the influence of and on the equilibration time and the equilibration deflection. Here the equilibration time means the time duration between the start and end of the diaphragm fluctuation. From the figures, we see that mainly controls the equilibration time of the diaphragm, while controls both the equilibration time and the equilibration deflection.

(a)

(b)

Figure 12: The influence of and on the dynamic response behavior of the diaphragm: (a) influence of at ; (b) influence of at .

(a)

(b)

(c)

Figure 13: The influence of and on the equilibration time and the equilibration deflection: (a) versus equilibration time; (b) versus equilibration deflection; (c) versus equilibration time.

5.3. Example  3: Thermomechanical In-Plane Microactuator

In this subsection, the TPWL-GMEC2 (POD) method is applied to a thermomechanical in-plane microactuator (TIM) [30, 31]. As shown in Figure 14(a), the TIM consists of a moveable shuttle connected to electrical contact pads anchored on the substrate by slender thermal expansion legs. The legs are connected to the shuttle at a slight inclination angle. As voltage is applied across the contact pads, current flows through the legs and shuttle. The high-current density in the legs causes heating and thermal expansion of the legs and hence causes a linear motion of the shuttle in the desired direction.

(a)

(b)

Figure 14: Thermomechanical in-plane microactuator: (a) the device; (b) a single leg-pair and the cross-section.

Because each leg-pair in Figure 14(a) has the same structure and is applied to the same voltage, modeling the entire actuator may be simplified as modeling a single leg-pair. The single leg-pair model can capture the functionality of the full TIM device and provide a significant reduction in model complexity. A schematic diagram of a single leg-pair is shown in Figure 14(b) and the corresponding geometric parameters are listed in Table 5.

Table 5: Geometric parameters of a single leg-pair.

Structure parameter
Value

Length of shuttle
 μm

Width of shuttle
 μm

Thickness of shuttle
 μm

Length of leg
 μm

Width of leg
 μm

Thickness of leg
 μm

Amount of inclination
 μm

Length of pad
 μm

Width of pad
 μm

Thickness of pad
 μm

Height of air
 μm

Height of Si3Ni4
 μm

Height of substrate
 μm

According to the thermoelectric coupled modeling approach [14], a governed equation for the single leg-pair may be written as

where is the potential, the temperature, the temperature dependent electric conductivity, the density, the specific heat, and the temperature dependent thermal conductivity. and are the electric field intensity and the electric current density, respectively, and can be calculated by

For the electric equation, the following boundary conditions are adopted:

where is the outward unit normal vector to the boundary. For the temperature equation, initial condition is set as with being the ambient temperature. Boundary conditions on the device surface (except the lower surface) are

where and are the heat fluxes generated by heat convection and heat radiation, respectively. Here is the coefficient of heat convection, is the heat radiation rate, and is the Stephan-Boltzmann constant. As to the lower surface of the device, because the air gap is very thin and flow can be ignored, heat may conduct from the device to the substrate. If we assume the temperature of the lower substrate surface is constant (i.e., ), then boundary condition on the lower surface of the device can be written as

where is the heat conduction flux, is the shape factor, is the distance from the lower surface of the device to the upper surface of the substrate, is the effective thermal conductivity from the device to the substrate. For the contact segment and the overhead segment, we have and and and , respectively. Here , , and are the thermal conductivities of air, Si3Ni4 layer, and Si substrate. Finally, may be evaluated as

where and are the corresponding thickness and width of the device.

In simulations, all the material parameters are from [30]. Equation (24) is discretized by a 3D finite element method on a mesh partition with 7,848 tetrahedral grid cells, which leads to a full nonlinear system of order . The time step size for solving the full nonlinear system is set as , and the current on the input electrode is chosen as the input (i.e., ). Details of the discretization process and application of the TPWL method can be found from our previous work [14].

Figure 15 shows the steady state of the single leg-pair for the input current  mA, where the locations −375~−275 μm and 275~375 μm correspond to the pads, the locations −275~−25 μm and 25~275 μm correspond to the legs, and the location −25~25 μm corresponds to the shuttle. As can be seen, the legs are much hotter than the shuttle due to their bigger current density (or greater potential drop from Figure 15(b)) and lower heat dissipating capacity. This is conducive to the motion of the shuttle and agrees well with [30].

(a)

(b)

Figure 15: Results of the full nonlinear model for  mA: (a) steady state temperature; (b) steady state potential.

Table 6 gives some results obtained by the TPWL-GMEC2 (POD) method for different . We find that, for the same accurate requirement , bigger input current always needs more LPs and higher order ROMs. For the fifteen cases listed in Table 6, the highest order of the ROMs is only 30 and no more than 5 LPs are needed. Even though the sizes of the ROMs are so small, the maximum errors of all ROMs are less than 0.01. These results demonstrate once again the good performance of the TPWL-GMEC2 (POD) method. In order to make the results given in Table 6 more intuitive, Figure 16 plots the results for the input current  mA.

Table 6: Results obtained by the TPWL-GMEC2 (POD) method for different .

1.0 mA

13
1

13
1

13
2

2.0 mA

16
1

16
2

16
3

3.0 mA

26
1

26
2

26
4

4.0 mA

26
1

26
2

26
4

5.0 mA

30
2

30
3

30
5

(a)

(b)

(c)

(d)

Figure 16: Results obtained by TPWL-GMEC2 (POD) method for  mA: (a) steady state temperature at the central axis; (b) steady state potential at the central axis; (c) transient average temperature of the leg; (d) maximum root mean square errors of the ROMs.

Table 7 shows the efficiency of the TPWL-GMEC2 (POD) method (CPU: AMD Radeon A6-3600, 2.1 GHz; RAM: 4 GB). From the table, the time cost for extracting the ROMs increases as decreases or the input current increases, whereas the maximum extracting time is only 6.443 seconds, which is much improved compared to the original TPWL-GMEC method (27.79 seconds). In addition, the ROMs generated by the TPWL-GMEC2 (POD) method can obtain thousands of times’ speed-up compared to the full nonlinear model. The maximum simulation time of the ROMs shown here is only 0.097 seconds, which benefits the applications to MEMS design.

Table 7: Efficiency of the TPWL-GMEC2 (POD) method.

Full simulation time (s)

Speed-up factor

1.0 mA
185.01

1.759
0.013
14231

1.759
0.013
14231

2.426
0.014
13215

2.0 mA
274.45

2.742
0.022
12475

3.435
0.026
10556

4.153
0.031
8853

3.0 mA
278.07

2.886
0.048
5793

3.777
0.057
4878

5.504
0.072
3862

4.0 mA
253.57

2.685
0.046
5512

3.547
0.054
4696

5.325
0.07
3622

5.0 mA
245.25

3.559
0.067
3660

4.675
0.077
3185

6.443
0.097
2528

At last, Table 8 gives the comparisons of the steady voltage and the peak temperature obtained by the ROMs generated by the TPWL-GMEC2 (POD) method and the full nonlinear model. As decreases, the difference between the results obtained by full nonlinear model and ROMs is reduced. When , the maximum difference is less than 0.1%. This demonstrates the effectiveness of the TPWL-GMEC2 (POD) method from another perspective.

Table 8: Comparisons of the steady voltage and peak temperature.

(mA)
Steady voltage (V)
Peak temperature (°C)

1
2
3
4
5
1
2
3
4
5

Full model
1.74
3.74
6.179
9.225
13.06
47.77
128.6
247.8
407.4
608.2

1.742
3.829
6.631
10.78
13.07
49.33
158.9
378.7
777.6
608.8

1.742
3.74
6.18
9.228
13.07
49.33
128.7
248
407.8
608.8

1.74
3.74
6.179
9.227
13.07
47.8
128.6
247.9
407.7
608.7

6. Conclusion

A more efficient algorithm for selecting high quality LPs in TPWL MOR has been presented by making two improvements to the TPWL-GMEC method. The first improvement is adding a preliminary LP selection procedure to reduce the size of the state-set before LPs selection without losing accuracy, and the second improvement is combining with the POD method to avoid updating the global projection matrix. Simulation results of a nonlinear transmission line circuit demonstrate that the time cost of the improved TPWL-GMEC method is much less than the original TPWL-GMEC method when the obtained ROMs are of the same quality. And then, the applications of the new method to MOR of three MEMS devices, which are a micromachined switch, an electrostatic micropump diaphragm, and a thermomechanical in-plane microactuator, once again fully demonstrate the good performance of the improved TPWL-GMEC method. Particularly, the first example (the micromachined switch) tests the ability of the improved method in dealing with strong nonlinearity; the second example (the electrostatic micropump diaphragm) demonstrates the improved method can be directly applied to second-order systems; and the third example (the thermomechanical in-plane microactuator) shows the performance of the improved method for complex 3D multifield coupling problems. On the whole, the improved TPWL-GMEC method can extract ROMs with high quality at low time cost and make up for the shortcoming of the original TPWL-GMEC method. Moreover, combining with POD method also brings some other benefits, such as lower memory requirement and enhanced flexibility.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the Natural Science Basic Research Plan in Shaanxi Province of China (2017JQ1035), the Scientific Research Program Funded by Shaanxi Provincial Education Department of China (17JK0104), and the China Postdoctoral Science Foundation.

References

W. Schilders, “Introduction to model order reduction,” in Model Order Reduction: Theory, Research Aspects and Applications, vol. 13 of Math. Ind., pp. 3–32, 2008.

M. Kudryavtsev, E. B. Rudnyi, J. G. Korvink, D. Hohlfeld, and T. Bechtold, “Computationally efficient and stable order reduction methods for a large-scale model of MEMS piezoelectric energy harvester,” Microelectronics Reliability, vol. 55, no. 5, pp. 747–757, 2015.

A. Corigliano, M. Dossi, and S. Mariani, “Domain decomposition and model order reduction methods applied to the simulation of multi-physics problems in MEMS,” Computers & Structures, vol. 122, pp. 113–127, 2013.

W. H. A. Schilders, “The need for novel model order reduction techniques in the electronics industry,” Lecture Notes in Electrical Engineering, vol. 74, pp. 3–23, 2011.

L. Feng, “Review of model order reduction methods for numerical simulation of nonlinear circuits,” Applied Mathematics and Computation, vol. 167, no. 1, pp. 576–591, 2005.

U. Baur, P. Benner, and L. Feng, “Model order reduction for linear and nonlinear systems: a system-theoretic perspective,” Archives of Computational Methods in Engineerin: State-of-the-Art Reviews, vol. 21, no. 4, pp. 331–358, 2014.

M. Rewienski and J. White, “A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices,” in Proceedings of the IEEE/ACM International Conference on Computer Aided Design. ICCAD 2001. IEEE/ACM Digest of Technical Papers, pp. 252–257, San Jose, CA, USA.

M. Rewienski and J. White, “A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 22, no. 2, pp. 155–170, 2003.

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.

N. Dong and J. Roychowdhury, “General-purpose nonlinear model-order reduction using piecewise-polynomial representations,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 27, no. 2, pp. 249–264, 2008.

B. N. Bond and L. Daniel, “Stable reduced models for nonlinear descriptor systems through piecewise-linear approximation and projection,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 28, no. 10, pp. 1467–1480, 2009.

Y.-L. Jiang and H.-B. Chen, “Application of general orthogonal polynomials to fast simulation of nonlinear descriptor systems through piecewise-linear approximation,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 31, no. 5, pp. 804–808, 2012.

S. A. Nahvi, Mashuq-Un-Nabi, and S. Janardhanan, “Trajectory based methods for nonlinear MOR: Review and perspectives,” in Proceedings of the 2012 IEEE International Conference on Signal Processing, Computing and Control (ISPCC '12), March 2012.

Y. Liu, W. Yuan, and H. Chang, “A global maximum error controller-based method for linearization point selection in trajectory piecewise-linear model order reduction,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 33, no. 7, pp. 1100–1104, 2014.

Y. Liu, W. Yuan, H. Chang, and B. Ma, “Compact thermoelectric coupled models of micromachined thermal sensors using trajectory piecewise-linear model order reduction,” Microsystem Technologies, vol. 20, no. 1, pp. 73–82, 2014.

T. Voß, Model reduction for nonlinear differential algebraic equations [Ph.D. thesis], University of Wuppertal, 2005.

T. Voß, R. Pulch, E. J. ter Maten, and A. El Guennouni, “Trajectory Piecewise Linear Approach for Nonlinear Differential-Algebraic Equations in circuit simulation,” in Scientific Computing in Electrical Engineering, vol. 11 of Mathematics in Industry, pp. 167–173, Springer Berlin Heidelberg, Berlin, Heidelberg, 2007.

K. Mohaghegh, Linear and nonlinear model order reduction for numerical simulation of electric circuits, Logos Verlag Berlin GmbH, 2010.

M. J. Rewienski, A trajectory piecewise-linear approach to model order reduction of nonlinear dynamical systems [Ph.D. thesis], Massachusetts In- stitute of Technology, 2003.

J. A. Martínez, Model order reduction of nonlinear dynamic systems using multiple projection bases and optimized state-space sampling [Ph.D. thesis], University of Pittsburgh, 2009.

X. Chen and Z. Wu, “Review on macromodels of MEMS sensors and actuators,” Microsystem Technologies, pp. 1–14, 2017.

C. G. Wu, Y. C. Liang, W. Z. Lin, H. P. Lee, and S. P. Lim, “A note on equivalence of proper orthogonal decomposition methods,” Journal of Sound and Vibration, vol. 265, no. 5, pp. 1103–1110, 2003.

D. Vasilyev, M. Rewienski, and J. White, “A TBR-based trajectory piecewise-linear algorithm for generating accurate low-order models for nonlinear analog circuits and MEMS,” in Proceedings of the 40th Design Automation Conference, pp. 490–495, June 2003.

E. Bertarelli, R. Ardito, A. Corigliano, and R. Contro, “A plate model for the evaluation of pull-in instability occurrence in electrostatic micropump diaphragms,” International Journal of Applied Mechanics, vol. 3, no. 1, pp. 1–19, 2011.

M. N. Albunni, V. Rischmuller, T. Fritzsche, and B. Lohmann, “Model-order reduction of moving nonlinear electromagnetic devices,” IEEE Transactions on Magnetics, vol. 44, no. 7, pp. 1822–1829, 2008.

Q. Wang, Macro-modeling of micro-electrical-mechanical system devices [Ph.D. thesis], Massachusetts Institute of Technology, 1998.

E. Hung, . Yao-Joe Yang, and S. Senturia, “Low-order models for fast dynamical simulation of MEMS microstructures,” in Proceedings of the International Solid State Sensors and Actuators Conference (Transducers '97), pp. 1101–1104, Chicago, IL, USA.

Z. Bai and Y. Su, “Dimension reduction of large-scale second-order dynamical systems via a second-order Arnoldi method,” SIAM Journal on Scientific Computing, vol. 26, no. 5, pp. 1692–1709, 2005.

B. Salimbahrami and B. Lohmann, “Order reduction of large scale second-order systems using Krylov subspace methods,” Linear Algebra and its Applications, vol. 415, no. 2-3, pp. 385–405, 2006.

C. D. Lott, T. W. McLain, J. N. Harb, and L. L. Howell, “Modeling the thermal behavior of a surface-micromachined linear-displacement thermomechanical microactuator,” Sensors and Actuators A: Physical, vol. 101, no. 1-2, pp. 239–250, 2002.

Y.-J. Yang and K.-Y. Shen, “Nonlinear heat-transfer macromodeling for MEMS thermal devices,” Journal of Micromechanics and Microengineering, vol. 15, no. 2, pp. 408–418, 2005.