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

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.


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][5][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][8][9][10][11][12][13][14][15][16][17][18][19][20][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][17][18][19][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.

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 x() ∈   represents the state vector, f :   →   is a nonlinear vector function, B ∈  × distributes the input u() ∈   , and C ∈  × maps x to the system response y() ∈   .
Let {x 0 , . . ., x  } be the state-set associated with a certain training input u() for  ∈ [ 0 ,   ], where x 0 , . . ., x  are approximate values of x() at time instants  0 <  1 < ⋅ ⋅ ⋅ <   =   .The LPs {x lin 0 , . . ., x lin −1 } ⊂ {x 0 , . . ., x  } are assumed having been selected out according to some rules, where  is the number of LPs.Through linearizing f(x) at those LPs, a set of the local linearized models give where A  is the Jacobian of f(x) evaluated at x = x lin  .Next, the local projection matrices V  for each local linearized model are generated by using a linear MOR method, such as Krylovsubspace methods [7,19] and truncated balanced realization (TBR) method [23].The global projection matrix V ∈  × is constructed by merging and orthogonalizing V  .According to V, the local linearized models are reduced to size ; namely, where where w (z) is the weighting function.For any z, it satisfies the conditions w (z) ∈ [0, 1] and ∑ −1 =0 w (z) = 1.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  = 25.
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 (1) Given a training input u(); (2) Simulate the nonlinear system (1) and store the states at all time steps x 0 , . . ., x  ; (3) Choose the initial state x 0 as the first LP x lin 0 , and set  = 1; (4) Update V and built the current TPWL model; (5) Simulate the current TPWL model to get the responses x0 , . . ., x ; (6) Find out the maximum error state x  : Set  max = ‖x 0 − x 0 ‖ RMSE and  = 0, where the root mean square error (RMSE) is defined by end for (7) if  max >  and x  has not been selected as a LP then select the ( + 1)th LP x lin  = x  , set  ←  + 1; (8) if  is increased in step (7) and  <  ( is the maximum allowable LP number) then go to step (4).( 9) Return the selected LPs {x lin 0 , x lin 1 , . . ., x lin −1 } and end the loop.
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.

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.

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, x 0 and x  are, respectively, the initial and final state vectors.Generally, for a given nonlinear system, a higher accuracy requirement usually needs smaller  (i.e., bigger ).

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.
(preliminary LP selection) (1) Given a training input u() and a threshold  > 0; (2) Choose the initial state as the first LP candidate x 0 , and set  = 1; (3) Simulate the nonlinear system (1)  In addition to improving efficiency, combining with POD method can also bring other benefits.First, when the Krylovsubspace 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.

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.
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.

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  = 1, inductance  = 10, capacitance  = 1, and the constitutive equation of diodes   (V) = exp(40V) − 1.The current source entering into node 1 is the input () = (), and the voltage at node 1 is the output response () = V 1 ().The system equations of the nonlinear transmission line RLC circuit can be represented as 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  = 500,  = 10, and Δ = 0.001, respectively.Other userdefined 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  = ‖x  − x 0 ‖/10.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.
For a more detailed comparison, Tables 1 and 2 show the order of ROMs , the LP number , the actual maximum error  max , the time cost for extracting the ROMs  extracting , and the simulation time of the ROMs  ROM of four TPWL-GMEC methods with different  for the single-and dualsinusoidal 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  ROM 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.
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 ‖x  − x 0 ‖/5.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 ‖x  − x 0 ‖/100 and ‖x  − x 0 ‖/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  = ‖x  − x 0 ‖/10 according to the above general rule of choosing .Here below are the details of applications.

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.
As shown in Figure 6, the micromachined switch consists of a polysilicon beam suspended over a silicon substrate.The beam has length  = 610 m, width  = 40 m, and height ℎ = 2.2 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,  elec 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. elec is approximated by with  0 the permittivity of vacuum and V the applied voltage.
Figure 7 shows the results obtained by TPWL-GMEC2 (POD) method for different .Two inputs considered in the simulations are V() = 7() (() = 0 for  < 0 and () = 1 for  ≥ 0) and V() = 7 cos() with  = 4 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  = 1 − 3, 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 V() = 7 cos(), 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.
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 V() = 7() when  = 1 − 4, 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 V() = 7() has good approximation to the periodic cosine input V() = 7 cos().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 V() = 5() and V() = 6(), small deviation for the input V() = 8(), and large deviation for the input V() = 9().Notice that the pull-in phenomenon occurs when V() = 9() and does not when V() = 7().Therefore, the ROM for V() = 7() does not fully contain the response information of V() = 9().When we use this ROM to simulate V() = 9(), the simulation result certainly becomes unsatisfactory.Nevertheless, the results in Figure 8 still show the good expandability of the TPWL-GMEC2 (POD) method.
Since the input V() = 9() 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.,  = ‖x  − x 0 ‖/10) is used, a satisfactory approximation is not obtained even when  = 1 − 4 (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  = 1 − 4. So, if we want to get good approximation,  should be decreased.The result in Figure 9(b) shows great improvement when  = ‖x  − x 0 ‖/30 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.

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.
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  0 .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,  = ℎ 3 /12(1 − ] 2 ) is the plate stiffness where  is Young's modulus and ] is Poisson's ratio of the material.And Δ 2 is the double Laplacian operator in polar space coordinates; that is, The electrostatic pressure under an applied potential V 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  0 the natural frequency of the structure.By defining dimensionless variables û = / 0 , r = /, and t =  0 , ( 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 0 ≤  ≤ 1 and  > 0.
In simulations, ( 19) is discretized by standard finite difference scheme with  = 500 equidistantly distributed grid points, and the corresponding compact form of the spatial discretized system can be written as where x = ( 1 ,  2 , . . .,   )  is the state vector, b(x) is the nonlinear input vector, M, E, and K are mass matrix, damping matrix, and stiffness matrix, respectively.The time step size for solving system ( 22) is set as Δ = 1.0×10 −3 , and the center point deflection of the diaphragm () = (0, ) is selected as the output.Among three dimensionless model parameters, the mass coefficient is fixed at  = 1.0, 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 Krylovsubspace method or transform second-order system to firstorder 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 V is the POD projection matrix and A  is the Jacobi of b(x) at x  .
Figure 11 shows the results obtained by TPWL-GMEC2 (POD) method for different .When  = 1 − 2, obvious difference can be found between the results of the ROMs and full nonlinear model.But when  is decreased to 1 − 3, 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  = 1−4, 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  = 1 − 3.
In Figures 11(a) and 11(b), two distinct response behaviors of the diaphragm can be observed.For  = 10, 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  = 14, 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.
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.

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 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 legpair.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.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.E and J 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 n is the outward unit normal vector to the boundary.For the temperature equation, initial condition is set as (, , , 0) =  0 with  0 being the ambient temperature.Boundary conditions on the device surface (except the lower surface) are where  conv = ℎ( 0 − ) and  rad = ( 4 0 −  4 ) 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.,  =  0 ), then boundary condition on the lower surface of the device can be written as where  cond =  eff ( 0 − )/ 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,  eff is the effective thermal conductivity from the device to the substrate.For the contact segment and the overhead segment, we have  eff = /(ℎ  /  + ℎ  /  ) and  = ℎ  + ℎ  and  eff = /(ℎ air / air () + ℎ  /  +ℎ  /  ) and  = ℎ air +ℎ  +ℎ  , respectively.Here  air (),   , 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  = 5582.The time step size for solving the full nonlinear system is set as Δ = 1.0 × 10 −6 , and the current on the input electrode is chosen as the input (i.e., () =  0 ).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  0 = 5mA, 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].
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  0 = 5 mA.
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  0 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.
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  = 1 − 4, the maximum difference is less than 0.1%.This demonstrates the effectiveness of the TPWL-GMEC2 (POD) method from another perspective.

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.
and C  = V  C. Finally, a TPWL ROM can be represented as a weighted combination of all local linearized ROMs z ()  = −1

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

Table 5 :
Geometric parameters of a single leg-pair.Structure parameter Value Length of shuttle   = 50 m Width of shuttle   = 30 m Thickness of shuttle ℎ  = 3.5 m Length of leg   = 250 m Width of leg   = 3.0 m Thickness of leg ℎ  = 3.5 m Amount of inclination  = 3 m Length of pad   = 100 m Width of pad   = 30 m Thickness of pad ℎ  = 5.5 m Height of air ℎ air = 2 m Height of Si3Ni4 ℎ  = 0.6 m Height of substrate ℎ  = 500 m and hence causes a linear motion of the shuttle in the desired direction.

Figure 13 :
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.

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

Figure 16 :
Figure 16: Results obtained by TPWL-GMEC2 (POD) method for  0 = 5 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 3 :
Results obtained by the original TPWL method proposed by Rewienski and White.

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

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

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