A Systematic Design Optimization Approach for Multiphysics MEMS Devices Based on Combined Computer Experiments and Gaussian Process Modelling

This paper presents a systematic and efficient design approach for the two degree-of-freedom (2-DoF) capacitive microelectromechanical systems (MEMS) accelerometer by using combined design and analysis of computer experiments (DACE) and Gaussian process (GP) modelling. Multiple output responses of the MEMS accelerometer including natural frequency, proof mass displacement, pull-in voltage, capacitance change, and Brownian noise equivalent acceleration (BNEA) are optimized simultaneously with respect to the geometric design parameters, environmental conditions, and microfabrication process constraints. The sampling design space is created using DACE based Latin hypercube sampling (LHS) technique and corresponding output responses are obtained using multiphysics coupled field electro–thermal–structural interaction based finite element method (FEM) simulations. The metamodels for the individual output responses are obtained using statistical GP analysis. The developed metamodels not only allowed to analyze the effect of individual design parameters on an output response, but to also study the interaction of the design parameters. An objective function, considering the performance requirements of the MEMS accelerometer, is defined and simultaneous multi-objective optimization of the output responses, with respect to the design parameters, is carried out by using a combined gradient descent algorithm and desirability function approach. The accuracy of the optimization prediction is validated using FEM simulations. The behavioral model of the final optimized MEMS accelerometer design is integrated with the readout electronics in the simulation environment and voltage sensitivity is obtained. The results show that the combined DACE and GP based design methodology can be an efficient technique for the design space exploration and optimization of multiphysics MEMS devices at the design phase of their development cycle.


Introduction
MEMS accelerometers are miniaturized sensors for measuring constant, time varying, and quasi-static accelerations and have wide applications in the field of automotive industry [1], machine condition monitoring [2,3], shock sensing [4], precision navigation [5], and consumer electronics [6]. For these applications, one basic requirement is that the MEMS accelerometer must be able to determine the position of a body in space by sensing its acceleration in three axes. This is generally achieved by designing a single axis MEMS accelerometer and mounting orthogonal to each other on a body for acceleration measurement. This results not only increased device footprint area and increased packaging cost, but also measurement error due to misalignment [7]. To resolve these issues, one of the solutions reported in the literature is the monolithic integration of three proof masses, each sensing in a specific axis, in a single chip [8][9][10]. Though, this approach minimizes the packaging cost, but achieving small device size, high sensitivity and small noise floor values remain the main limitations. In comparison to multiple proof masses integrated, MEMS accelerometers based on single proof mass for sensing multiple axis acceleration have proved to be an efficient solution for achieving a small device footprint, low cost, and improved performance [11][12][13].
MEMS are multiphysics devices and generally require the optimization of multiple output responses for a given set of design parameters and microfabrication process constraints. The optimization of MEMS is traditionally carried out by varying one factor at a time and analyzing its effect on an output response either by developing analytical models, FEM models, or topology optimization [14][15][16][17][18][19]. These methods are efficient for MEMS devices with simple geometric configuration but for MEMS with relatively complex geometry and with the requirement of multiphysics design space exploration and multiple output responses to be optimized simultaneously, these methods become inefficient due to high computational costs and modelling complexity.
The optimization of a capacitive MEMS accelerometer design is a multiphysics problem involving electro-thermal-structural interactions. The main performance responses of a capacitive MEMS accelerometer that must be optimized include natural frequency, proof mass displacement, pull-in voltage, capacitance change, and thermomechanical noise. The natural frequency and proof mass displacement are strongly dependent on the stiffness of mechanical suspension beams. For high sensitivity, the natural frequency value should be low which requires low stiffness mechanical suspension beams or large proof mass. However, low mechanical stiffness of suspension beams leads to a low value of pull-in voltage which further limits the maximum value of initial bias voltage that can be applied to the sensing capacitive parallel plates and maximum input acceleration range. Similarly, for high resolution, a large value of capacitance change is desired in a MEMS accelerometer which is dependent on the initial air gap between the sensing parallel plates and their overlap area. A large value of overlap area and small air gap results in large capacitance change but decreases the pull-in voltage threshold. Moreover, the air damping between the sensing parallel plates increases by increasing the overlap area of plates and decreasing their initial gap. Thus, decreasing the quality factor and hence increasing the thermomechanical noise. Thus, it becomes important that for the geometric design optimization of MEMS accelerometer, all the output responses be considered simultaneously instead of the traditional approach of varying one factor at a time and analyzing its effect on a single output response presented in the literature [20][21][22].
The design of experiments (DoE) is a systematic statistical approach for optimization problems and has been a widely used optimization technique in many different fields, for example in manufacturing processes [23], sensor optimization [21,24], and precision agriculture [25]. The traditional DoE approach explores the design space of an optimization problem with a minimum number of physical experiments and the metamodels for output responses are developed by estimating the effect of design parameters using least square regression which is based on the randomness due to experimental variations [26]. Previously, the authors have presented the application of traditional DoE based optimization methodology for MEMS devices [27,28]. The application of the traditional DoE technique based on least square regression for deterministic computer simulations has been a matter of discussion in the literature since it lacks random error and the method of least squares residuals has no statistical meaning [29]. However, Simpson et al. [30] have argued that the application of traditional DoE for computer simulations is a trade-off of appropriateness vs. practicality. Sacks et al. [31] pointed out that the selection of design points for deterministic computer simulations is a statistical experimental design problem and also proposed a Gaussian process (GP) model by considering the deterministic output responses as random stochastic processes.
The design and analysis of computer experiments (DACE) based on GP models are now a popular choice for the development of metamodels, using deterministic computer simulations, in different fields and allow the design space exploration with minimum computational cost [32][33][34][35]. For MEMS devices, the DACE based design optimization, combined with the FEM based MEMS simulation tools and latest high computing machines, can greatly minimize the MEMS development cycle time and reduce the costs involved in the traditional iterative microfabrication runs for the realization of a fully functional MEMS device.
In this paper, we present a DACE based systematic and efficient design methodology for MEMS in general and MEMS accelerometer in particular by using Latin hypercube sampling (LHS) technique to create a design space with different combinations of geometric design parameters and GP based metamodelling for the multi-response optimization. The optimization for the MEMS accelerometer is carried out by considering the microfabrication process constraints and environmental operating conditions. Figure 1 shows the schematic of the proposed two degree-of-freedom (2-DoF) MEMS accelerometer design conceptualized considering the constraints imposed by the commercially available SOIMUMPs microfabrication process offered by MEMSCAP Inc. USA [36]. The main constraints of the SOIMUMPs process are (a) no bottom electrode is available for the out-of-plane proof mass displacement sensing and (b) the design cannot have anchors in the middle of the proof mass, thus limiting the position of the mechanical springs to the outer side of the mass. For sensing the proof mass displacement, corresponding to an input acceleration, parallel plate combs are attached to each of the four sides of the proof mass. The stator combs are located between the two rotor combs, attached to the proof mass, in a gap-antigap configuration. To obtain maximum capacitance change, the ratio of the larger gap to the smaller gap between combs is kept at 3 with gap and antigap values of 2.5 µm and 7.5 µm, respectively. The suspension beams are T-shaped and are designed symmetrically on the four corners of the proof mass to minimize the cross-axis coupling and increase stability. To minimize the effect of sudden high acceleration values, end stoppers are designed on the proof mass corners. The capacitance sensing mechanism is based on the gap closing mechanism and both stator and rotor combs are attached in differential sensing configuration to obtain the maximum capacitance change. The thickness of the structural silicon-on-insulator (SOI) layer in the whole structure is 25 µm.

Gaussian Process Modelling
A generalized form of Gaussian distribution is termed as Gaussian process (GP), which is a type of continuous stochastic process. It is used to describe the probability distribution over functions and the model developed by GP is a conditional probabilistic model. When fitting data from a deterministic computer experiment, the Gaussian process model is frequently utilized. Sacks et al. [31] presented these models as models for computer experiments. They are desired because they give an exact fit to the data gathered from computer experiments and only need a single parameter for every design parameter. The optimization design space can be represented as d = ( , y ) of k runs, where i = 1, 2, … . , k; is the D-dimensional vector of input design parameters for each run and y is the scalar value of the output response. The design parameters and the output responses for the k runs can be written in a matrix = , , … . . , and y = y , y , … . . , y , respectively. The Gaussian process model for an output response is defined as [26]: where μ is the mean of the modelled surface and z( ) is the Gaussian process, as a function of design parameters . The value of output response from any collection of design parameters has a multivariate normal distribution, which defines a Gaussian process. This multivariate normal distribution has a mean μ , where being the vector of ones of k length. The multivariate normal distribution also contains a covariance matrix that is proportional to a correlation matrix with a specific structure that allows the points to form a smooth surface. The covariance is: where σ is the variance which is the proportionality constant between the covariance matrix and the correlation matrix. The matrix ( , ) represents the correlation matrix and is a function of hyperparameters and design parameters. The elements of the correlation matrix can be estimated using the relation [37]:

Gaussian Process Modelling
A generalized form of Gaussian distribution is termed as Gaussian process (GP), which is a type of continuous stochastic process. It is used to describe the probability distribution over functions and the model developed by GP is a conditional probabilistic model. When fitting data from a deterministic computer experiment, the Gaussian process model is frequently utilized. Sacks et al. [31] presented these models as models for computer experiments. They are desired because they give an exact fit to the data gathered from computer experiments and only need a single parameter for every design parameter. The optimization design space can be represented as d s = (x i , y i ) of k runs, where i = 1, 2, . . . , k; x i is the D-dimensional vector of input design parameters for each run and y i is the scalar value of the output response. The design parameters and the output responses for the k runs can be written in a matrix X = [x 1 , x 2 , . . . . . . , x k ] T and y = [y 1 , y 2 , . . . . . . , y n ] T , respectively. The Gaussian process model for an output response is defined as [26]: where µ is the mean of the modelled surface and z(x i ) is the Gaussian process, as a function of design parameters x i . The value of output response from any collection of design parameters has a multivariate normal distribution, which defines a Gaussian process. This multivariate normal distribution has a mean µ1 n , where 1 n being the vector of ones of k length. The multivariate normal distribution also contains a covariance matrix that is proportional to a correlation matrix with a specific structure that allows the points to form a smooth surface. The covariance is: where σ 2 is the variance which is the proportionality constant between the covariance matrix and the correlation matrix. The matrix R(X, θ) represents the correlation matrix and is a function of hyperparameters θ and design parameters. The elements of the correlation matrix can be estimated using the relation [37]: In Equation (3), θ t ≥ 0 (t = 1, 2, . . . , D) the term r ij specifies the correlation between the output responses at any two points. The correlation matrix is built so that the distance between any two different locations in the input domain is inversely proportional to the correlation of their outputs. In specifically, as the distance between the two input locations approaches zero, the correlation increases to one, and as the distance approaches infinity, the correlation decreases to zero.
The prediction equation to evaluate the value of an output response at other sets of design parameters is given as follows [38]: where r T = [r(x 1 , x), r(x 2 , x), . . . , r(x n , x),] and the values forθ andμ are used in Equation (4) after estimating the maximum likelihood estimates of parameters µ and θ. Figure 2 shows the DACE based optimization methodology implemented in this paper for the multiresponse optimization of MEMS accelerometer. The first step is to identify the main output responses of the MEMS accelerometer to be optimized, which is followed by the selection of the geometric design parameters and operating conditions with the specification of low and high values for the design parameters. Since in the present study, the number of design parameters is eight, thus a statistical Latin Hypercube Sampling (LHS) based technique for the efficient sampling of the design points for the design parameters is implemented. The output responses for each simulation run in the LHS-based design matrix are obtained through Multiphysics FEM simulations. Based on the FEM simulation results for all the design points in the LHS based design matrix, a Gaussian process regression (GPR) analysis is performed and metamodels for all the output responses are obtained. The analysis of these metamodels for the output responses gives a detailed insight into the effect of the individual design parameters and their interaction on the output responses. The GPR analysis allows to analyze a single output response considering the design parameters but to simultaneously optimize all the output responses with respect to the design parameters, an objective function has been defined and a desirability function based approach has been implemented. Finally, the optimized values of the geometric design parameters and corresponding output responses are verified through FEM simulations. If the predicted values from the optimization results are not within 95% confidence interval, then it is recommended to return to the step of generating sample data space with more data points using LHS to obtain a more accurate response surface.   Table 1 shows the design parameters with respective low and high levels considered for the DACE based optimization of MEMS accelerometer. The design parameters are coded as X for i = 1, 2, … . ,8. The low and high levels of the design parameters are based on the geometric configuration of design, operating environmental conditions, and microfabrication process constraints. Figure 1 shows the dimensions of the T-shape spring in terms of lengths and width and represents the design parameters X , X , and X . For the design parameter X , the minimum width that can be set for a beam having a length greater than 100 µm is 6 µm as per the SOIMUMPs process constraint. The proposed MEMS accelerometer is to be optimized for measuring the input acceleration in the range of ±1 g to ±25 g.   Table 1 shows the design parameters with respective low and high levels considered for the DACE based optimization of MEMS accelerometer. The design parameters are coded as X I for i = 1, 2, . . . , 8. The low and high levels of the design parameters are based on the geometric configuration of design, operating environmental conditions, and microfabrication process constraints. Figure 1 shows the dimensions of the T-shape spring in terms of lengths and width and represents the design parameters X 2 , X 3 , and X 4 . For the design parameter X 4 , the minimum width that can be set for a beam having a length greater than 100 µm is 6 µm as per the SOIMUMPs process constraint. The proposed MEMS accelerometer is to be optimized for measuring the input acceleration in the range of ±1 g to ±25 g. Length of suspension beam 1 400 µm 500 µm X 3 Length of suspension beam 2 400 µm 500 µm X 4 Width of suspension beam 6 µm 8 µm X 5 Input acceleration 1 g 25 g X 6 Operating temperature 233.15 K 373.15 K X 7 Operating pressure 100 Torr 760 Torr X 8 Frequency ratio 0.1 0.5 The MEMS accelerometers for most of the applications, especially in space applications are required to operate in the temperature range of −40 • C (233.15 K) to 100 • C (373.15 K). Therefore, the low and high levels for the design parameter operating temperature (X 6 ) are kept in this range. The range for design parameter operating pressure (X 7 ) is kept 100 Torr (sub-atmospheric pressure) for a low level to 760 Torr (atmospheric pressure) for a high level. This range is selected to analyze the dependency of noise and other output responses of the MEMS accelerometer on the operating pressure. The last design parameter X 8 is the frequency ratio (FR), which is the ratio of any frequency value (lower than natural frequency) to the obtained natural frequency of the specific design run.

Latin Hypercube Sampling (LHS) Based Space Filling Design
The selection of design space for computer experiments is based on two basic principles: (a) any combination of design parameters should not appear more than once in the design space since replication is only required in physical experiments to account for random errors and (b) the different combinations of design parameters should cover the whole region of the design space, so that different behaviors of output responses in different areas of the design space can be analyzed [39]. This is also because due to the complex nature of computer simulations, the response behavior may change across the design space. The space-filling designs are generally used for DACE due to the evenly spread of design points throughout the design space. Among various types of space-filling designs, the Latin hypercube sampling (LHS) approach is most widely used. The LHS creates the design space for design parameters by maximizing the minimum distance between any two design points and the design points evenly cover the whole design space. This ensures enough degree of freedom to estimate both the linear and quadratic effects of design parameters and minimize the discrepancy between the observed values. For the LHS and other space-filling design, it is recommended that the number of simulation runs must be at least 10X, where X is the number of design parameters [26]. Table S1 shows the LHS based design matrix for the design parameters considered for the optimization of the MEMS accelerometer.

FEM Modelling of 2-DoF MEMS Accelerometer
The output responses, corresponding to the LHS based design matrix, are obtained through FEM simulations in CoventorWare ®® software. The MEMS accelerometer is meshed using three-dimensional (3D) solid tetrahedral 739,599 elements with a total number of 272,580 elements for the central proof mass and elements for the mechanical suspension beams. A fine mesh is used for the mechanical suspension beams, with multiple elements along the thickness, as shown in Figure 3. The material properties of thin film Silicon are included as input in the FEM analysis with Young's Modulus of 169 GPa, Poisson ratio of 0.29, and density of 2300 kg/m 3 [36]. The natural frequency of the MEMS accelerometer is obtained through modal analysis in the CoventorWare MemMech module. The proof mass displacement, corresponding to an input acceleration, is obtained for each simulation run using harmonic analysis and considering the air damping effect. In the dual-axis MEMS accelerometer design, both the  The natural frequency of the MEMS accelerometer is obtained through modal analysis in the CoventorWare MemMech module. The proof mass displacement, corresponding to an input acceleration, is obtained for each simulation run using harmonic analysis and considering the air damping effect. In the dual-axis MEMS accelerometer design, both the squeeze and slide air film damping result in energy dissipation. For an input acceleration, the squeeze film air damping occurs between the stator and rotor combs in the active sense axis, and a slide film air damping occurs between the two sets of combs in the inactive axis. In the present work, for each simulation run, the value of air damping is first computed by using the CoventorWare DampingMM module and subsequently used in the harmonic analysis. In the FEM analysis, the gas rarefaction effects are considered for each simulation run since the dominant airflow regime is slip flow for an air gap of 2.5 µm between the stator and rotor combs of the MEMS accelerometer. For each simulation run in the design matrix, the computed values of squeeze number and Reynolds number are in the range of 10 −3 thus air compressibility and inertial damping effects are ignored. A detailed description of modelling of thin film air damping effects in capacitive sensing combs is given by the authors in [40,41].
For a bias voltage of 2.25 V to the parallel capacitive sensing plates, coupled field electric-structural analysis is performed to find the value of pull-in voltage. The capacitance change for each simulation run is obtained through the CoventorWare MEMS+ module. Based on the FEM analysis based on natural frequency and quality factor, the value of BNEA for each simulation run, at a given operating temperature, is obtained by using the following relation [42].
where K B is the Boltzmann constant, T is the operating temperature, ω n is the natural frequency, m p is the value of proof mass, and Q is the quality factor.

Development of GP Based Metamodels for the Output Responses
For the LHS based design matrix, with corresponding output responses obtained through FEM simulations (Table S1), the Gaussian process model is fitted to obtain the metamodels for the output responses. For each output response, the value of hyperparameter θ is obtained for each design parameter by minimizing the negative log-likelihood and the values are reported in Table 2. These values are used in Equation (3) to control the length scale and smoothness of modelled Gaussian surface for each output response. After fitting the Gaussian process model, the predictive values for the output response are obtained based on Equation (4) and interactions between design parameters for each output response are also estimated. Table 3 lists the values for interaction between significant design parameters for each output response.

Significant Design Parameters and Interaction Analysis for Natural Frequency (Y 1 )
From the Gaussian process model analysis, the most significant interaction is between the design parameters X 3 and X 4 for the natural frequency with a value of 0.005, as given in Table 3. Figure 4 shows the 3D response surface plot for the natural frequency of the MEMS accelerometer with respect to design parameters X 3 and X 4 . The response surface plot is obtained using GP based metamodel and the value of all other design parameters are kept at their mean value. The result shows that the natural frequency is more sensitive to the change in X 4 in comparison to X 3 . Moreover, with X 4 = 8 µm the natural frequency is much more sensitive to change in the X 3 value in comparison to when X 4 = 6 µm.

Output Responses Significant Design Parameters Interactio Design Parameters
Interaction Va Natural frequency (Y )

Significant Design Parameters and Interaction Analysis for Natural Frequency ( )
From the Gaussian process model analysis, the most significant interactio tween the design parameters X and X for the natural frequency with a value as given in Table 3. Figure 4 shows the 3D response surface plot for the natural fr of the MEMS accelerometer with respect to design parameters X and X . The r surface plot is obtained using GP based metamodel and the value of all other de rameters are kept at their mean value. The result shows that the natural frequency sensitive to the change in X in comparison to X . Moreover, with X = 8 µm the frequency is much more sensitive to change in the X value in comparison to w 6 µm.

Significant Design Parameters and Interaction Analysis for Proof Mass Displacement (Y 2 )
The maximum interaction value is between the two design parameters X 4 and X 5 for the output response proof mass displacement. Figure 5 shows the 3D surface plot for the proof mass displacement with respect to the design parameters X 4 and X 5 with all other design parameters at their mean value. The plot shows that for input acceleration above 20 g, the proof mass displacement becomes more sensitive to the change in the suspension beam width. Moreover, the proof mass displacement shows a linear increase in the value with the increase in the input acceleration.
for the output response proof mass displacement. Figure 5 shows the 3D surface the proof mass displacement with respect to the design parameters X and X other design parameters at their mean value. The plot shows that for input acce above 20 g, the proof mass displacement becomes more sensitive to the change in pension beam width. Moreover, the proof mass displacement shows a linear inc the value with the increase in the input acceleration.

Significant Design Parameters and Interaction Analysis for Pull-In Voltage ( )
For the output response pull-in voltage, the design parameters interaction highest between X and X . Figure 6 shows the 3D surface plot for the interactio sis between these two design parameters on the pull-in voltage value. The plot sh the pull-in voltage is more sensitive to the change in design parameter X in com to X . Moreover, the effect of change in X value from 6 µm to 8 µm has less the pull-in voltage value when X is equal to 250 µm in comparison to when X 150 µm.

Significant Design Parameters and Interaction Analysis for Pull-In Voltage (Y 3 )
For the output response pull-in voltage, the design parameters interaction value is highest between X 1 and X 4 . Figure 6 shows the 3D surface plot for the interaction analysis between these two design parameters on the pull-in voltage value. The plot shows that the pull-in voltage is more sensitive to the change in design parameter X 4 in comparison to X 1 . Moreover, the effect of change in X 4 value from 6 to 8 µm has less effect on the pull-in voltage value when X 1 is equal to 250 µm in comparison to when X 1 value is 150 µm.

Significant Design Parameters and Interaction Analysis for Capacitance Change (
As per the GP model analysis, the interaction between two design param and X for the capacitance change has a maximum magnitude with a value of 0. 3D surface plot in Figure 7 shows that when X is at a high level of 8 µm the change in X on capacitance change is less in comparison to when X is at a low

Significant Design Parameters and Interaction Analysis for Capacitance Change (Y 4 )
As per the GP model analysis, the interaction between two design parameters X 4 and X 5 for the capacitance change has a maximum magnitude with a value of 0.038. The 3D surface plot in Figure 7 shows that when X 4 is at a high level of 8 µm the effect of change in X 5 on capacitance change is less in comparison to when X 4 is at a low level of 6 µm. Moreover, the effect of change in X 4 on the capacitance change is more when X 5 is at a high value of 25 g in comparison to its low value of 1 g. Figure 6. Response surface plot for pull-in voltage as a function of design parameters X

Significant Design Parameters and Interaction Analysis for Capacitance Change (
As per the GP model analysis, the interaction between two design param and X for the capacitance change has a maximum magnitude with a value of 0. 3D surface plot in Figure 7 shows that when X is at a high level of 8 µm the change in X on capacitance change is less in comparison to when X is at a low 6 µm. Moreover, the effect of change in X on the capacitance change is more w is at a high value of 25 g in comparison to its low value of 1 g.

Significant Design Parameters and Interaction Analysis for BNEA ( )
For the output response BNEA, the design parameters that have a significa on the output response BNEA are X , X , and X , respectively. The two design ter interactions X X and X X are comparable to each other with values 0.007 0.007235, respectively, for the output response BNEA. Figure 8a shows the interac of design parameters X and X for the BNEA. The results show that BNEA i linearly with an increase in the X from 150 µm to 250 µm. This can be attribute fact that with the increase in the overlap value the air damping increases whic quently decreases the quality factor. From Equation (5), it is clear that a decrea

Significant Design Parameters and Interaction Analysis for BNEA (Y 5 )
For the output response BNEA, the design parameters that have a significant effect on the output response BNEA are X 1 , X 6 , and X 7 , respectively. The two design parameter interactions X 1 X 7 and X 6 X 7 are comparable to each other with values 0.007285 and 0.007235, respectively, for the output response BNEA. Figure 8a shows the interaction plot of design parameters X 1 and X 7 for the BNEA. The results show that BNEA increases linearly with an increase in the X 1 from 150 µm to 250 µm. This can be attributed to the fact that with the increase in the overlap value the air damping increases which subsequently decreases the quality factor. From Equation (5), it is clear that a decrease in the quality factor leads to an increase in the BNEA value. The effect of an increase in X 7 from a low value of 100 Torr to the atmospheric air pressure of 760 Torr has a highly non-linear effect on the BNEA value for the MEMS accelerometer. Figure 8b shows the interaction plot of the design parameters X 6 and X 7 for the BNEA. The results show that the BNEA increases linearly with an increase in the design parameter X 6 value from 250 K to 375 K. quality factor leads to an increase in the BNEA value. The effect of an increase in X from a low value of 100 Torr to the atmospheric air pressure of 760 Torr has a highly non-linear effect on the BNEA value for the MEMS accelerometer. Figure 8b shows the interaction plot of the design parameters X and X for the BNEA. The results show that the BNEA increases linearly with an increase in the design parameter X value from 250 K to 375 K.

Prediction Accuracy of the Fitted GP Metamodels
The prediction accuracy of the developed GP based metamodels for the output responses is analyzed by estimating the mean absolute error (MAE), root mean square error (RMSE), and correlation coefficient I for each response. The MAE, RMSE, and R for an output response can be calculated as follows [43,44]:

Prediction Accuracy of the Fitted GP Metamodels
The prediction accuracy of the developed GP based metamodels for the output responses is analyzed by estimating the mean absolute error (MAE), root mean square error (RMSE), and correlation coefficient I for each response. The MAE, RMSE, and R for an output response can be calculated as follows [43,44]: Sensors 2021, 21, 7242 where y oi is the actual observed value for an output response for the ith run of the simulation design matrix and y pi is the predicted value of output response for the ith run. The MAE gives an estimation of the mean of the absolute errors for each simulation run, whereas RMSE gives a measure of the standard deviation of the residuals and depicts that how far the predicted points are from the fitted linear regression line. The smaller the MAE and RMSE values are, the better the points fit actual data values. The R-value close to one indicates that the linear relationship between the observed and the predicted values is positive, and the points lie nearly along a fit line. Table 4 lists the values of MAE, RMSE, and R for the output responses. These values suggest that the developed metamodels for the output responses are accurate in predicting the output responses.

Optimization Objective Function
The first step in the multi-response optimization of the proposed MEMS accelerometer design is to define the desired optimization objective function, for a given set of constraints for the design parameters. The objective function for the MEMS accelerometer is given as follows: Minimize-Natural frequency (Y 1 ) The value of the design parameter X 5 is set to a maximum input acceleration of 25 g. Similarly, the design parameters X 6 and X 7 are set to the atmospheric air pressure and room temperature of 760 Torr and 300 K, respectively.

Desirability Function Based Simultaneous Multi-Response Optimization
Various techniques have been reported in the literature to simultaneously optimize multiple output responses, including the distance function approach [45], loss function approach [46], and desirability function approach [47]. The desirability function approach was initially proposed by Harrington [48] in the form of exponential functions and later modified by Derringer and Suich [47] and Del Castillo et al. [49]. In this approach, the estimated value for each output response is transformed to a scale free value (d i (y i (x))) which is termed as desirability. The d i (y i (x )) value is scaled to be between 0 and 1 with 0 and 1 being the least and most optimal values, respectively. The overall desirability function is defined by taking the geometric mean of the individual desirability value for each output response. The desirability functions proposed by Derringer and Suich [47] contain non-differentiable target points and hence only search methods can be used for optimization. Later, Del Castillo et al. [49] proposed alternative piece-wise continuous desirability functions that account for non-differentiable points, and hence more efficient gradientbased algorithms can be applied to obtain the optimal solution. The non-differentiable desirability function for an output response is given as follows: where L, U, and T represent the lower, upper, and target values for an output response, respectively. The term δ = (U − L)/50 defines the small range around the non-differentiable point. The terms a 0 , a 1 , b 0 , and b 1 are constants and a detailed discussion to find the values of these constants are given in Del Castillo et al. [49]. The function f(y i (x)) is the polynomial approximation function that corrects for the non-differentiable points and is given as follows: A, B, C, D, and E are constant parameters. For simultaneous optimization of all the responses, a function named as global desirability is defined by taking the geometric mean of all the individual desirability values for the output responses and is given as follows: (12) where k is the number of output responses to be optimized and z i (0 < z i < 1) reflects the importance of each output response relative to others. The importance value (z i ) for each output response is scaled so that they sum up to one for all the responses.
The optimization of the desirability function is generally carried out using either search or gradient-based algorithms [47,49]. The search algorithm-based optimization is a derivative-free approach and can be applied for the optimization of desirability functions whose derivative does not exist. The search algorithm-based optimization requires the functions to have continuous first derivatives and is a more efficient and widely used method. The multi-response optimization of the desired objective function for the proposed 2-DoF MEMS accelerometer is performed using a gradient descent algorithm. Figure 9 shows the optimal values for the design parameters along with their individual desirability values. Figure 10 shows the individual desirability values for the output responses and corresponding optimized predicted values for each response. The maximum value of overall desirability obtained is 0.688 with a predicted value of 3037.37 Hz for the natural frequency, 0.9029 µm for the proof mass displacement, 6.7618 V for the pull-in voltage, 676.213 fF for the capacitance change, and 0.8061 µg/ √ Hz for the BNEA. value of overall desirability obtained is 0.688 with a predicted value of 3037.37 Hz for the natural frequency, 0.9029 µm for the proof mass displacement, 6.7618 V for the pull-in voltage, 676.213 fF for the capacitance change, and 0.8061 µg/√Hz for the BNEA.

Verification of Predicted Values for the Output Responses
The desirability function based predicted values for the output responses are further verified through FEM simulations to validate both the metamodel and desirability function based optimization approach for the proposed MEMS accelerometer design. The MEMS accelerometer is modelled by using the optimized geometric design parameter values shown in Figure 9.

Natural Frequency Analysis
The predicted optimized value of the natural frequency is verified through FEM based simulation in the CoventorWare MemMech module. Figure 11 shows the natural

Verification of Predicted Values for the Output Responses
The desirability function based predicted values for the output responses are further verified through FEM simulations to validate both the metamodel and desirability function based optimization approach for the proposed MEMS accelerometer design. The MEMS accelerometer is modelled by using the optimized geometric design parameter values shown in Figure 9.

Natural Frequency Analysis
The predicted optimized value of the natural frequency is verified through FEM based simulation in the CoventorWare MemMech module. Figure 11 shows the natural

Verification of Predicted Values for the Output Responses
The desirability function based predicted values for the output responses are further verified through FEM simulations to validate both the metamodel and desirability function based optimization approach for the proposed MEMS accelerometer design. The MEMS accelerometer is modelled by using the optimized geometric design parameter values shown in Figure 9.

Natural Frequency Analysis
The predicted optimized value of the natural frequency is verified through FEM based simulation in the CoventorWare MemMech module. Figure 11 shows the natural frequency of the MEMS accelerometer with the corresponding mode shape. The natural frequency value obtained through simulation is 3038.1 Hz, which lies within the 95% confidence prediction interval i.e., 2980 Hz ≤ Y 1 = 3038.1 Hz ≤ 3090 Hz. frequency of the MEMS accelerometer with the corresponding mode shape. The natural frequency value obtained through simulation is 3038.1 Hz, which lies within the 95% confidence prediction interval i.e., 2980 Hz ≤ Y = 3038.1 Hz ≤ 3090 Hz. Figure 11. Eigen mode shape for MEMS accelerometer having 1st mode along Y-direction at 3038.133 Hz.

Frequency Response Analysis
The predicted value of the output response proof mass displacement is verified by using the harmonic analysis in the CoventorWare MEMS+ module. Figure 12 shows the frequency response curve for the MEMS accelerometer at an input acceleration of 25 g. The results show that at frequency ratio (X ) of 0.5, the displacement amplitude of the proof mass is 0.898 µm. lies within the 95% prediction confidence prediction interval, i.e., 0.891 µm ≤ Y = 0.898 µm ≤ 0.915 µm. The X = 0.5 corresponds to an operational bandwidth of 1519 Hz. However, the MEMS accelerometers are generally designed for operational bandwidth of 0-450 Hz. Figure 12 shows that for X = 0.15 with operational bandwidth of 0-450 Hz, the proof mass displacement amplitude is relatively linear with a maximum value of 0.688 µm for the proof mass displacement.  Figure 13 shows the bias voltage vs. proof mass displacement graph for the final optimized MEMS accelerometer design. The results show that the maximum bias voltage that can be applied to the stator and rotor sensing combs must be less than the maximum pull-in value of 6.79 V. This pull-in voltage value is in close agreement with a predicted value of 6.76 V and lies within the 95% confidence prediction interval, i.e., 6.36 V ≤ Y = 6.79 V ≤ 7.16 V.

Frequency Response Analysis
The predicted value of the output response proof mass displacement is verified by using the harmonic analysis in the CoventorWare MEMS+ module. Figure 12 shows the frequency response curve for the MEMS accelerometer at an input acceleration of 25 g. The results show that at frequency ratio (X 8 ) of 0.5, the displacement amplitude of the proof mass is 0.898 µm. lies within the 95% prediction confidence prediction interval, i.e., 0.891 µm ≤ Y 2 = 0.898 µm ≤ 0.915 µm. The X 8 = 0.5 corresponds to an operational bandwidth of 1519 Hz. However, the MEMS accelerometers are generally designed for operational bandwidth of 0-450 Hz. Figure 12 shows that for X 8 = 0.15 with operational bandwidth of 0-450 Hz, the proof mass displacement amplitude is relatively linear with a maximum value of 0.688 µm for the proof mass displacement. frequency of the MEMS accelerometer with the corresponding mode shape. The natural frequency value obtained through simulation is 3038.1 Hz, which lies within the 95% confidence prediction interval i.e., 2980 Hz ≤ Y = 3038.1 Hz ≤ 3090 Hz. Figure 11. Eigen mode shape for MEMS accelerometer having 1st mode along Y-direction at 3038.133 Hz.

Frequency Response Analysis
The predicted value of the output response proof mass displacement is verified by using the harmonic analysis in the CoventorWare MEMS+ module. Figure 12 shows the frequency response curve for the MEMS accelerometer at an input acceleration of 25 g. The results show that at frequency ratio (X ) of 0.5, the displacement amplitude of the proof mass is 0.898 µm. lies within the 95% prediction confidence prediction interval, i.e., 0.891 µm ≤ Y = 0.898 µm ≤ 0.915 µm. The X = 0.5 corresponds to an operational bandwidth of 1519 Hz. However, the MEMS accelerometers are generally designed for operational bandwidth of 0-450 Hz. Figure 12 shows that for X = 0.15 with operational bandwidth of 0-450 Hz, the proof mass displacement amplitude is relatively linear with a maximum value of 0.688 µm for the proof mass displacement.  Figure 13 shows the bias voltage vs. proof mass displacement graph for the final optimized MEMS accelerometer design. The results show that the maximum bias voltage that can be applied to the stator and rotor sensing combs must be less than the maximum pull-in value of 6.79 V. This pull-in voltage value is in close agreement with a predicted value of 6.76 V and lies within the 95% confidence prediction interval, i.e., 6.36 V ≤ Y = 6.79 V ≤ 7.16 V.  Figure 13 shows the bias voltage vs. proof mass displacement graph for the final optimized MEMS accelerometer design. The results show that the maximum bias voltage that can be applied to the stator and rotor sensing combs must be less than the maximum pull-in value of 6.79 V. This pull-in voltage value is in close agreement with a predicted value of 6.76 V and lies within the 95% confidence prediction interval, i.e., 6.36 V ≤ Y 3 = 6.79 V ≤ 7.16 V.

Capacitance Change Analysis
For the verification of the predicted capacitance change for the final optimized MEMS accelerometer design, the accelerometer model developed in the CoventorWare MEMS+ module is interfaced with MATLAB Simulink. A bias voltage of 2.25 V is applied between the proof mass and stator electrodes, since the proposed accelerometer is designed to be interfaced with the commercial capacitive to voltage conversion readout IC, MS3110 [50]. An input acceleration of 25 g at X = 0.5, i.e., 1519 Hz is given as an input to the MEMS accelerometer, as shown in Figure 14a. Figure 14b,c show the corresponding gap and antigap capacitance. The net change in capacitance is 694 fF which is close to the predicted value of 676.2 fF and lies within the 95% confidence prediction interval, i.e., 618 fF ≤ Y = 694 fF ≤ 735 fF. For X = 0.15, i.e., an input acceleration frequency of 450 Hz, the net change in capacitance for the MEMS accelerometer is 523 fF.

Capacitance Change Analysis
For the verification of the predicted capacitance change for the final optimized MEMS accelerometer design, the accelerometer model developed in the CoventorWare MEMS+ module is interfaced with MATLAB Simulink. A bias voltage of 2.25 V is applied between the proof mass and stator electrodes, since the proposed accelerometer is designed to be interfaced with the commercial capacitive to voltage conversion readout IC, MS3110 [50]. An input acceleration of 25 g at X 8 = 0.5, i.e., 1519 Hz is given as an input to the MEMS accelerometer, as shown in Figure 14a. Figure 14b,c show the corresponding gap and anti-gap capacitance. The net change in capacitance is 694 fF which is close to the predicted value of 676.2 fF and lies within the 95% confidence prediction interval, i.e., 618 fF ≤ Y 4 = 694 fF ≤ 735 fF. For X 8 = 0.15, i.e., an input acceleration frequency of 450 Hz, the net change in capacitance for the MEMS accelerometer is 523 fF.

Capacitance Change Analysis
For the verification of the predicted capacitance change for the final optimized MEMS accelerometer design, the accelerometer model developed in the CoventorWare MEMS+ module is interfaced with MATLAB Simulink. A bias voltage of 2.25 V is applied between the proof mass and stator electrodes, since the proposed accelerometer is designed to be interfaced with the commercial capacitive to voltage conversion readout IC, MS3110 [50]. An input acceleration of 25 g at X = 0.5, i.e., 1519 Hz is given as an input to the MEMS accelerometer, as shown in Figure 14a. Figure 14b,c show the corresponding gap and antigap capacitance. The net change in capacitance is 694 fF which is close to the predicted value of 676.2 fF and lies within the 95% confidence prediction interval, i.e., 618 fF ≤ Y = 694 fF ≤ 735 fF. For X = 0.15, i.e., an input acceleration frequency of 450 Hz, the net change in capacitance for the MEMS accelerometer is 523 fF.

Estimation of Brownian Noise Equivalent Acceleration (BNEA)
For the estimation of BNEA, the MEMS 2-DoF accelerometer with optimal parameters is first subjected to damping analysis in the CoventorWare DampingMM module to obtain squeeze and slide film damping coefficients. The operating air pressure and temperature are set at room conditions, i.e., 300 K and 760 Torr. These coefficients are then used to estimate the value of BNEA based on Equation (5). The value of BNEA is estimated to be 0.805 µg/ √ Hz which is close to the predicted value of 0.8061 µg/ √ Hz and lies within the 95% confidence interval, i.e., 0.797 µg/ √ Hz ≤ Y 5 = 0.805 µg/ √ Hz ≤ 0.815 µg/ √ Hz.

Discussion
The capacitive MEMS accelerometers are multiphysics sensors involving complex interactions between the electrical-mechanical-fluidic domains. Table 5 shows a comparison of different capacitive MEMS accelerometer design and analysis approaches presented in the literature. Most of the MEMS accelerometer design analysis and optimization approaches involve the traditional varying one or two design factors at a time and observing its effect on a specific output response. In most cases, the design analysis is only focused on the geometric design parameters. However, the performance of capacitive MEMS accelerometers is strongly dependent on the sensor operating conditions including air pressure and temperature. The MEMS accelerometer operating temperature and air pressure have a strong influence on the air damping and hence on the dynamic response of these sensors. The design analysis and optimization of multiphysics MEMS devices in general and capacitive MEMS accelerometers, in particular, requires analyzing the output responses related to different physics domains simultaneously with respect to both geometric design parameters and device operating conditions. One of the important performance parameters for MEMS accelerometers is voltage sensitivity which is dependent on the capacitive to voltage readout electronics to estimate the voltage sensitivity of the optimized MEMS accelerometer, the behavioral model of the accelerometer is interfaced with the commercially available capacitance to voltage convertor Universal Readout IC TM MS3110 [50] in MATLAB Simulink environment. This readout circuitry is capable of measuring both single and differential capacitance with a resolution of 4 aF/ √ Hz. Figure 15 illustrates the integration of 2-DoF MEMS accelerometer behavioral model with MS3110. The working principle of the IC for the measurement of differential capacitance is based on charge amplification, followed by sample and hold, a low-pass frequency filter, and a buffer amplifier. A bias voltage of 2.25 V is applied between the proof mass and stator electrodes of the MEMS accelerometer. For an input acceleration along X-axis, the capacitance mismatch between the input capacitors, the difference CS 1 − CS 2 is zero and ∆C will be equal to the difference between the capacitors Cap1 and Cap3 as obtained in the capacitance change analysis section. The feedback capacitor C F value is adjusted to 3.2 pF to keep the output voltage in the range of 0.5 to 4 V. Figure 16 shows the effect of change in input acceleration on the output voltage of the readout interface circuit. The graph shows that the maximum and minimum values for output voltage are 3.93 and 0.573 V respectively with a voltage sensitivity of 65.4 mV/g for the final optimized MEMS accelerometer design. the voltage sensitivity of the optimized MEMS accelerometer, the behavioral model of the accelerometer is interfaced with the commercially available capacitance to voltage convertor Universal Readout IC TM MS3110 [50] in MATLAB Simulink environment. This readout circuitry is capable of measuring both single and differential capacitance with a resolution of 4 aF/√Hz. Figure 15 illustrates the integration of 2-DoF MEMS accelerometer behavioral model with MS3110. The working principle of the IC for the measurement of differential capacitance is based on charge amplification, followed by sample and hold, a low-pass frequency filter, and a buffer amplifier. A bias voltage of 2.25 V is applied between the proof mass and stator electrodes of the MEMS accelerometer. For an input acceleration along X-axis, the capacitance mismatch between the input capacitors, the difference CS − CS is zero and ∆C will be equal to the difference between the capacitors Cap1 and Cap3 as obtained in the capacitance change analysis section. The feedback capacitor C value is adjusted to 3.2 pF to keep the output voltage in the range of 0.5 to 4 V. Figure 16 shows the effect of change in input acceleration on the output voltage of the readout interface circuit. The graph shows that the maximum and minimum values for output voltage are 3.93 and 0.573 V respectively with a voltage sensitivity of 65.4 mV/g for the final optimized MEMS accelerometer design.

Conclusions
In this paper, a simulation-based design optimization approach for the MEMS accelerometer is presented by using DACE based design space sampling and development of metamodels for output responses of interest using GP modelling. The DACE based LHS sampling allows to cover the whole design space of the MEMS accelerometer by using only 10X simulation runs, where X = 8 is the number of design parameters. Moreover, the voltage sensitivity of the optimized MEMS accelerometer, the behavioral model of the accelerometer is interfaced with the commercially available capacitance to voltage convertor Universal Readout IC TM MS3110 [50] in MATLAB Simulink environment. This readout circuitry is capable of measuring both single and differential capacitance with a resolution of 4 aF/√Hz. Figure 15 illustrates the integration of 2-DoF MEMS accelerometer behavioral model with MS3110. The working principle of the IC for the measurement of differential capacitance is based on charge amplification, followed by sample and hold, a low-pass frequency filter, and a buffer amplifier. A bias voltage of 2.25 V is applied between the proof mass and stator electrodes of the MEMS accelerometer. For an input acceleration along X-axis, the capacitance mismatch between the input capacitors, the difference CS − CS is zero and ∆C will be equal to the difference between the capacitors Cap1 and Cap3 as obtained in the capacitance change analysis section. The feedback capacitor C value is adjusted to 3.2 pF to keep the output voltage in the range of 0.5 to 4 V. Figure 16 shows the effect of change in input acceleration on the output voltage of the readout interface circuit. The graph shows that the maximum and minimum values for output voltage are 3.93 and 0.573 V respectively with a voltage sensitivity of 65.4 mV/g for the final optimized MEMS accelerometer design.

Conclusions
In this paper, a simulation-based design optimization approach for the MEMS accelerometer is presented by using DACE based design space sampling and development of metamodels for output responses of interest using GP modelling. The DACE based LHS sampling allows to cover the whole design space of the MEMS accelerometer by using only 10X simulation runs, where X = 8 is the number of design parameters. Moreover, Figure 16. The output voltage of the readout electronics for an input acceleration in the range of −25 g to +25 g along the X-axis.

Conclusions
In this paper, a simulation-based design optimization approach for the MEMS accelerometer is presented by using DACE based design space sampling and development of metamodels for output responses of interest using GP modelling. The DACE based LHS sampling allows to cover the whole design space of the MEMS accelerometer by using only 10X simulation runs, where X = 8 is the number of design parameters. Moreover, combined DACE based sampling and GP modelling allows to analyze the interaction between the design parameters and their effect on the output responses simultaneously instead of conventional varying one factor at a time based optimization approach. The output responses of the MEMS accelerometer including natural frequency, proof mass displacement, pull-in voltage, capacitance change, and BNEA are obtained using FEM and behavioral model simulations for each simulation run in the LHS sampling design matrix. The Gaussian process regression is used for developing the metamodel for each output response which allowed to analyze the effect of design parameters and their interaction on the output responses in detail. The accuracy of the developed metamodels for predicting the output responses is proved by estimating the statistical parameters. An optimization objective function for the MEMS accelerometer is defined, and simultaneous optimization of the output responses is carried out using a combined desirability function and gradient-based algorithm. The predicted optimal values of the output responses are verified through FEM simulation and the obtained values showed a close agreement with the prediction. The optimized values of output responses for the MEMS accelerometer are natural frequency of 3036.37 Hz, proof mass displacement of 0.903 µm at ±25 g, pull-in voltage value of 6.762 V, output capacitance change of 676.213 fF at ±25 g, and BNEA of 0.8061 µg/ √ Hz. Based on the optimized values of the proof mass displacement and capacitance change and for an operational bandwidth of 0-450 Hz, the mechanical, capacitance, and voltage sensitivity of the MEMS accelerometer are 0.027 µm/g, 22.9 fF/g, and 65.4 mV/g, respectively.
The DACE based optimization technique for the multi-response optimization of MEMS accelerometer using Gaussian process regression overcomes the limitations of traditional DoE based optimization and can be implemented for multiphysics MEMS devices with complex geometric configurations. In the future, the application of DACE based optimization methodology considering the microfabrication process uncertainties and size effects can be investigated for the reliability-based optimization of the MEMS accelerometers.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/s21217242/s1, Table S1: The sets of design parameters for computer experiments on MEMS 2-DoF accelerometer design optimization.