Python-Based Open-Source Electro-Mechanical Co-Optimization System for MEMS Inertial Sensors

The surge in fabrication techniques for micro- and nanodevices gave room to rapid growth in these technologies and a never-ending range of possible applications emerged. These new products significantly improve human life, however, the evolution in the design, simulation and optimization process of said products did not observe a similarly rapid growth. It became thus clear that the performance of micro- and nanodevices would benefit from significant improvements in this area. This work presents a novel methodology for electro-mechanical co-optimization of micro-electromechanical systems (MEMS) inertial sensors. The developed software tool comprises geometry design, finite element method (FEM) analysis, damping calculation, electronic domain simulation, and a genetic algorithm (GA) optimization process. It allows for a facilitated system-level MEMS design flow, in which electrical and mechanical domains communicate with each other to achieve an optimized system performance. To demonstrate the efficacy of the methodology, an open-loop capacitive MEMS accelerometer and an open-loop Coriolis vibratory MEMS gyroscope were simulated and optimized—these devices saw a sensitivity improvement of 193.77% and 420.9%, respectively, in comparison to their original state.


Introduction
The design, simulation, and optimization process of a MEMS device include a sequence of designing the initial geometry, mechanical parameter simulation and optimization, designing of the electrical interface, and simulation of the complete system [1].
To design, simulate, and optimize MEMS inertial sensors, engineers typically separate the process in two very distinct-yet symbiotic-domains: mechanical and electrical domains. This workflow is often strictly divided and a great deal of simplification is applied to one of the domains in order to allocate computational resources to achieve a complete simulation and optimization of the other [2].
MEMS mechanical structures are designed with a computer-assisted design (CAD) software and commonly comprise thousands of degrees of freedom (DoF) which lead to a heavy computational cost when simulating mechanical behavior. To bypass this obstacle, engineers have used reduced-order modelling methods to build system-level models [3][4][5]-bringing the thousands of DoFs down to a few, frequently used when designing closed-loop control systems [6]. Younis et al. [7] developed a reduced-order model to study the behavior of electrostatically actuated microbeams-based MEMS, using a macro model to reduce the computational time. The aforementioned method can be helpful when designing the sensor electrical interface, but it fails to take into consideration the full complex mechanical structure, and consequently the interaction between electrical and mechanical domains.
Moreover, the typical optimization methodology combines multiphysics software and a programming language interpreter program. Wang et al. [8,9] presented a MEMS mechan-

Finite Element Method (FEM)
The finite element method was implemented in the proposed co-optimization system to simulate complex mechanical geometries.
The FEM approaches any simulation problem by subdividing a continuous entity into finite smaller parts, solving each one individually, and reassembling them. Many physical processes, and especially most solid mechanics ones, can be described by partial differential equations (PDEs). For a computer to solve PDEs, it applies the FEM to divide a complex system into smaller subparts-finite elements. This division process is called space discretization, and the generation of a mesh. This is a way of transcribing a 2D or 3D object into a series of mathematical points that can be analyzed. For static solid mechanics studies, the main category of PDEs are elliptic, which can be solved using a variational FEM method [13]. A variational method has its basis on the principle of energy minimization: when a boundary condition is applied, the configuration where the total energy is minimum is the one that prevails. This process starts with the multiplication of PDEs by a test function, then integrate the resulting equation over the domain, and finally perform integration by parts with second-order derivatives [13].

Python Language and Libraries
The Python programming language was used to develop the co-optimization system in this study. It is a high-level language and well suited for scientific and engineering environments. Its highly modular nature and clean syntax provide a simple and direct code writing suitable in many scientific applications [14]. The main advantage of Python language lies in the countless number of library modules provided by either official Python or the global developer's community, which offers numerous potential combinations and applications.

Python Simulation and Optimization Software
A complete MEMS co-simulation and co-optimization program comprises different essential blocks. This type of software needs a geometry designer and processor with meshing abilities, a FEM simulation block powerful enough to process different mesh sizes with varying degrees of complexity, a dedicated electrical domain script capable of interpreting the mechanical results of each MEMS design, and a GA that takes into consideration both mechanical and electrical performance parameters.
A software covering all the aforementioned blocks was developed with a general structure depicted in Figure 1. The software tool also comprised a viscous damping calculation. In this way, the optimization considers an additional important performance parameter. interpreting the mechanical results of each MEMS design, and a GA that takes into consideration both mechanical and electrical performance parameters. A software covering all the aforementioned blocks was developed with a general structure depicted in Figure 1. The software tool also comprised a viscous damping calculation. In this way, the optimization considers an additional important performance parameter.

MEMS Geometry Design in Python
A competent geometry design and meshing tool require several fundamental characteristics: the ability to create different shapes and perform Boolean operations on them, the capacity to generate a customizable mesh to be assigned to the created geometry, and the possibility to import and export files. Two Python libraries were chosen-pygmsh [15,16] and meshio [17]. Pygmsh is used for geometry building and mesh generation, the meshio library is applied to generate a file that can be read by the other software blocks.

FEM Simulation for Displacement and Modal Analysis
The FEM is employed in this tool to solve PDEs problems involving linear elasticity equations: the first is to calculate the displacement resulting from an applied force, and the second is to perform a modal analysis-obtaining the MEMS eigenfrequency modes.
When a force is applied to a body Ω, the equations describing the related deformations of the geometries with isotropic elastic conditions are described in [18].
To study MEMS inertial sensors, the knowledge of natural frequencies and the corresponding mode shapes are essential. Modal analysis solves an eigensystem-a group of all eigenvectors belonging to a linear transformation matched with the corresponding eigenvalue. In an eigensystem, the eigenvector represents the mode shape and the eigenvalue represents the frequency.

Electronic Domain Simulation
The electronic domain simulation block is the third. In the readout circuit, the number of comb fingers and overlapping area varies with the geometry, and the actuation scheme changes for each design. Two MEMS inertial sensors were simulated, i.e., a capacitive accelerometer and a linear vibratory gyroscope. The MEMS accelerometer has a proof-mass between two electrodes, which results in differential sensing defined by Equation (2). The MEMS vibratory gyroscope has a proof-mass between two sets of comb-fingers which moves along the y-axis, modifying the gap between the fixed electrodes and the moving ones as described by Equation (4). In this group of equations, N stands for the number of comb fingers, t is the thickness, L represents the comb fingers length, and d signifies the gap between the fingers.

MEMS Geometry Design in Python
A competent geometry design and meshing tool require several fundamental characteristics: the ability to create different shapes and perform Boolean operations on them, the capacity to generate a customizable mesh to be assigned to the created geometry, and the possibility to import and export files. Two Python libraries were chosen-pygmsh [15,16] and meshio [17]. Pygmsh is used for geometry building and mesh generation, the meshio library is applied to generate a file that can be read by the other software blocks.

FEM Simulation for Displacement and Modal Analysis
The FEM is employed in this tool to solve PDEs problems involving linear elasticity equations: the first is to calculate the displacement resulting from an applied force, and the second is to perform a modal analysis-obtaining the MEMS eigenfrequency modes.
When a force is applied to a body Ω, the equations describing the related deformations of the geometries with isotropic elastic conditions are described in [18].
To study MEMS inertial sensors, the knowledge of natural frequencies and the corresponding mode shapes are essential. Modal analysis solves an eigensystem-a group of all eigenvectors belonging to a linear transformation matched with the corresponding eigenvalue. In an eigensystem, the eigenvector represents the mode shape and the eigenvalue represents the frequency.

Electronic Domain Simulation
The electronic domain simulation block is the third. In the readout circuit, the number of comb fingers and overlapping area varies with the geometry, and the actuation scheme changes for each design. Two MEMS inertial sensors were simulated, i.e., a capacitive accelerometer and a linear vibratory gyroscope. The MEMS accelerometer has a proof-mass between two electrodes, which results in differential sensing defined by Equation (2). The MEMS vibratory gyroscope has a proof-mass between two sets of comb-fingers which moves along the y-axis, modifying the gap between the fixed electrodes and the moving ones as described by Equation (4). In this group of equations, N stands for the number of comb fingers, t is the thickness, L represents the comb fingers length, and d signifies the gap between the fingers.
To drive the gyroscope's proof-mass into resonance, electrostatic actuation is applied. Electrostatic actuators rely on the force between two electrodes when a voltage is applied between them [19]. Parallel-plate actuation electrodes are commonly built to apply a force in a specific direction-aligned with the desired motion direction and DoF of the target mass.
In the designed MEMS gyroscope, a balanced actuation mechanism was applied. The actuating comb-fingers generate the desired actuation force by sliding parallel to each other, as described in Equation (5)-a potential V 1 = V DC + V AC is applied to one set of electrodes and another potential V 2 = V DC − V AC to the opposing set.
The calculation of the sensing and actuating mechanism is included in the software to study the influence of electrical parameters on the device performance through the genetic algorithm and thus achieving electro-mechanical co-optimization.
A capacitance-to-voltage converter was implemented in the code, to calculate a voltage output from the acceleration induced capacitance change. The converter is based on the converting block designed by Utz et al. [20], as shown in Figure 2 with a governing equation stated in Equation (6), in which V DD is considered 5 V and V CM is considered 2.5 V.
, , To drive the gyroscope's proof-mass into resonance, electrostatic actuation is applied. Electrostatic actuators rely on the force between two electrodes when a voltage is applied between them [19]. Parallel-plate actuation electrodes are commonly built to apply a force in a specific direction-aligned with the desired motion direction and DoF of the target mass.
In the designed MEMS gyroscope, a balanced actuation mechanism was applied. The actuating comb-fingers generate the desired actuation force by sliding parallel to each other, as described in Equation (5)-a potential is applied to one set of electrodes and another potential to the opposing set.

(5)
The calculation of the sensing and actuating mechanism is included in the software to study the influence of electrical parameters on the device performance through the genetic algorithm and thus achieving electro-mechanical co-optimization.
A capacitance-to-voltage converter was implemented in the code, to calculate a voltage output from the acceleration induced capacitance change. The converter is based on the converting block designed by Utz et al. [20], as shown in Figure 2 with a governing equation stated in Equation (6), in which VDD is considered 5 V and VCM is considered 2.5 V.

Damping Calculation
In this work, the damping was modeled as viscous damping, which is the primary contributor to overall damping of the system. This type of damping mechanism occurs when the gas surrounding the vibratory structures introduces viscous effects caused by the internal friction of the gas trapped in the middle of vibratory structures such as comb fingers. Here, two viscous film damping effects, i.e., squeeze film damping and slide film damping, were modelled to calculate the quality factor of the drive and sense modes for both the accelerometer and gyroscope. The simulated MEMS inertial sensors were considered to be surrounded by air (εr = 1) at the atmospheric pressure of 1013.25 mbar.

Damping Calculation
In this work, the damping was modeled as viscous damping, which is the primary contributor to overall damping of the system. This type of damping mechanism occurs when the gas surrounding the vibratory structures introduces viscous effects caused by the internal friction of the gas trapped in the middle of vibratory structures such as comb fingers. Here, two viscous film damping effects, i.e., squeeze film damping and slide film damping, were modelled to calculate the quality factor of the drive and sense modes for both the accelerometer and gyroscope. The simulated MEMS inertial sensors were considered to be surrounded by air (ε r = 1) at the atmospheric pressure of 1013.25 mbar.
In Equation (7), K n represents the Knudsen number which is a measure of gas rarefaction effect. λ f is the ratio of the mean free path and d refers to the gap containing the gas. This equation makes the ambient pressure of the sensor's working environment. Equation (8) denotes the effective viscosity of squeeze film damping in which µ stands for the mean viscosity. Equation (9) describes the effective viscosity of slide film damping. Equation (10) refers to the viscous damping coefficient. Finally, the quality factor (Q) is determined by Equation (11), in which ω denotes angular frequency, and M is the moving mass.

Genetic Algorithm Optimization
The genetic algorithm (GA) block conducts the electro-mechanical co-optimization. In this study, the genetic algorithm was developed with the following workflow, as illustrated in Figure  For the individuals of the next generation, both an integral copy and mutated copy of the 25 best individuals in the first generation are included and the remaining 50 devices are randomly generated; 4.
Process 3 is repeated for a number of generations-until half of the population converges to a high-performance FOM and an individual is selected.
mented [21,22]: In Equation (7), Kn represents the Knudsen number which is a measure of gas rarefaction effect. λf is the ratio of the mean free path and d refers to the gap containing the gas. This equation makes the ambient pressure of the sensor's working environment. Equation (8) denotes the effective viscosity of squeeze film damping in which µ stands for the mean viscosity. Equation (9) describes the effective viscosity of slide film damping. Equation (10) refers to the viscous damping coefficient. Finally, the quality factor (Q) is determined by Equation (11), in which ω denotes angular frequency, and M is the moving mass.

Genetic Algorithm Optimization
The genetic algorithm (GA) block conducts the electro-mechanical co-optimization. In this study, the genetic algorithm was developed with the following workflow, as illustrated in Figure

Case Study 1: MEMS Capacitive Accelerometer
To show the effectiveness of the developed software, as the first demonstration, an open-loop capacitive MEMS accelerometer is designed, simulated, and optimized. This device comprises a proof-mass suspended by four beams above the substrate. The accelerometer is designed to measure an acceleration in the z-axis by detecting the displacement of the proof-mass in the z-axis, with a mass-spring-damper model illustrated in Figure 4a. ment of the proof-mass in the z-axis, with a mass-spring-damper model illustrated in Figure 4a.
To detect the capacitance change introduced by the acceleration input, the proofmass is located between two electrodes with an overlap area equal to the proof-mass surface area. The initial distance between the proof-mass and the electrodes changes when the proof mass experiences a displacement caused by an acceleration. The MEMS accelerometer structure and its fundamental resonant frequency is illustrated in Figure 4b, with geometric parameters listed in Table 1. The device comprises four L-shaped beams connected to the proof-mass on one end and fixed on the other, which suspend the proof-mass above the substrate. The simulation of the MEMS capacitive accelerometer is based on the flow chart illustrated in Figure 5.  To detect the capacitance change introduced by the acceleration input, the proof-mass is located between two electrodes with an overlap area equal to the proof-mass surface area. The initial distance between the proof-mass and the electrodes changes when the proof mass experiences a displacement caused by an acceleration.
The MEMS accelerometer structure and its fundamental resonant frequency is illustrated in Figure 4b, with geometric parameters listed in Table 1. The device comprises four L-shaped beams connected to the proof-mass on one end and fixed on the other, which suspend the proof-mass above the substrate.

Optimization Results
The simulation of the MEMS capacitive accelerometer is based on the flow chart illustrated in Figure 5.
To show the effectiveness of the developed software, as the first demonstration, an open-loop capacitive MEMS accelerometer is designed, simulated, and optimized. This device comprises a proof-mass suspended by four beams above the substrate. The accelerometer is designed to measure an acceleration in the z-axis by detecting the displacement of the proof-mass in the z-axis, with a mass-spring-damper model illustrated in Figure 4a.
To detect the capacitance change introduced by the acceleration input, the proofmass is located between two electrodes with an overlap area equal to the proof-mass surface area. The initial distance between the proof-mass and the electrodes changes when the proof mass experiences a displacement caused by an acceleration.  The MEMS accelerometer structure and its fundamental resonant frequency is illustrated in Figure 4b, with geometric parameters listed in Table 1. The device comprises four L-shaped beams connected to the proof-mass on one end and fixed on the other, which suspend the proof-mass above the substrate. The simulation of the MEMS capacitive accelerometer is based on the flow chart illustrated in Figure 5.  The MEMS accelerometer generates a detectable capacitance change, which is then read by a capacitance-to-voltage circuit, governed by Equation (6). For the simulated accelerometer, V DD is defined as 5 V, V CM = 2.5 V, and C int = 300 fF.
The solution of a PDE is strongly related to the density of the mesh. It is, therefore, necessary to perform a mesh convergence study. In this case, the natural frequency is analyzed with different numbers of meshing elements in the suspension beams. As observed in Figure 6, the meshing elements reached the optimal number at 33,824-after this, for the next 6 data points, the variation in the first frequency mode becomes less than 0.15%. Thus, the remaining simulation and optimization process use a number of elements of 33,824 as the optimized meshing element size.
The solution of a PDE is strongly related to the density of the mesh. It is, therefore, necessary to perform a mesh convergence study. In this case, the natural frequency is analyzed with different numbers of meshing elements in the suspension beams. As observed in Figure 6, the meshing elements reached the optimal number at 33,824-after this, for the next 6 data points, the variation in the first frequency mode becomes less than 0.15%. Thus, the remaining simulation and optimization process use a number of elements of 33,824 as the optimized meshing element size. The genetic algorithm optimization process described in Section 2.3.5 is applied to this device for 6 generations, with a FOM defined by Equation (12). In this equation, Sensitivity stands for the output voltage when an acceleration of 1 g is applied, Frequency is the resonant frequency and Qfactor denotes the quality factor of the sensing mechanism under the atmosphere air pressure. The FOM is comprised by weighted elements-allowing the designers to choose the weight of each performance parameter according to the projected application of the MEMS device.
Within 6 generations (100 individuals in each generation), the GA altered the chosen initial geometric parameters: proof-mass length and suspension beam width and obtained an optimized device. The geometric changes and performance are listed in Table 2. The evolution of the device in Figure 7 illustrates the algorithm tends to reduce the suspension beam width and enlarge the proof-mass size. That increases the sensitivity of the MEMS accelerometer by 193.77% (due to lower stiffness in the suspension system combined with a larger proof-mass) but decreases the resonant frequency by 42.54% and the quality factor by 51.51%. Finally, the FOM is increased by 26.79% after the optimization process.
To verify the accuracy of the FEM analysis performed by the proposed software, a COMSOL simulation of the same device was made. The natural frequency of the MEMS accelerometer obtained by COMSOL was 1897.15 Hz, while the natural frequency of the MEMS accelerometer obtained by the proposed software was 1887.9 Hz. The relative difference of the two simulated natural frequencies was 0.5%, proving the accuracy of the FEM analysis performed by the proposed software.  The genetic algorithm optimization process described in Section 2.3.5 is applied to this device for 6 generations, with a FOM defined by Equation (12).
In this equation, Sensitivity stands for the output voltage when an acceleration of 1 g is applied, Frequency is the resonant frequency and Q factor denotes the quality factor of the sensing mechanism under the atmosphere air pressure. The FOM is comprised by weighted elements-allowing the designers to choose the weight of each performance parameter according to the projected application of the MEMS device.
Within 6 generations (100 individuals in each generation), the GA altered the chosen initial geometric parameters: proof-mass length and suspension beam width and obtained an optimized device. The geometric changes and performance are listed in Table 2. The evolution of the device in Figure 7 illustrates the algorithm tends to reduce the suspension beam width and enlarge the proof-mass size. That increases the sensitivity of the MEMS accelerometer by 193.77% (due to lower stiffness in the suspension system combined with a larger proof-mass) but decreases the resonant frequency by 42.54% and the quality factor by 51.51%. Finally, the FOM is increased by 26.79% after the optimization process.
To verify the accuracy of the FEM analysis performed by the proposed software, a COMSOL simulation of the same device was made. The natural frequency of the MEMS accelerometer obtained by COMSOL was 1897.15 Hz, while the natural frequency of the MEMS accelerometer obtained by the proposed software was 1887.9 Hz. The relative difference of the two simulated natural frequencies was 0.5%, proving the accuracy of the FEM analysis performed by the proposed software.

Case Study 2: Linear MEMS Vibratory Gyroscope
The second demonstrator of the proposed software was a MEMS vibratory gyroscope, reproduced from [22]. This device featured a drive frame implemented to nest the

Case Study 2: Linear MEMS Vibratory Gyroscope
The second demonstrator of the proposed software was a MEMS vibratory gyroscope, reproduced from [22]. This device featured a drive frame implemented to nest the proofmass and thus decouples the drive and sense motion. A u-beam suspension system was put together to ensure that both the drive and sense motion only deflect in the correct direction.
The MEMS gyroscope maintains an oscillation in the drive axis and experienced a Coriolis force under an angular-rate input. The Coriolis force will result in an energy transfer from the drive axis to the sense axis which occurs in the form of a proof-mass movement along this axis.
The geometry built by the software is shown in Figure 8b, with initial geometric parameters listed in Table 3. This design comprises eight anchors, illustrated in red: suspension u-beams anchors, stationary electrodes for sensing, and stationary drive electrodes.

Case Study 2: Linear MEMS Vibratory Gyroscope
The second demonstrator of the proposed software was a MEMS vibratory gyroscope, reproduced from [22]. This device featured a drive frame implemented to nest the proof-mass and thus decouples the drive and sense motion. A u-beam suspension system was put together to ensure that both the drive and sense motion only deflect in the correct direction.
The MEMS gyroscope maintains an oscillation in the drive axis and experienced a Coriolis force under an angular-rate input. The Coriolis force will result in an energy transfer from the drive axis to the sense axis which occurs in the form of a proof-mass movement along this axis.
The geometry built by the software is shown in Figure 8b, with initial geometric parameters listed in Table 3. This design comprises eight anchors, illustrated in red: suspension u-beams anchors, stationary electrodes for sensing, and stationary drive electrodes.    The u-beams are designed to perform as a suspension system that keeps the proof-mass above the substrate. The four drive frame beams allow movement of the device along the x-axis (the drive direction), and the four beams that connect the proof-mass to the drive frame facilitate a displacement by the y-axis (the sense direction).
In this gyroscope, drive oscillation along the x-axis was actuated through the lateral electrodes shown in Figure 8b. Differential capacitive sensing was used in the device.

Optimization Results
The simulation and optimization of the vibratory MEMS gyroscope is based on the system-level model illustrated in Figure 9.
The u-beams are designed to perform as a suspension system that keeps the proofmass above the substrate. The four drive frame beams allow movement of the device along the x-axis (the drive direction), and the four beams that connect the proof-mass to the drive frame facilitate a displacement by the y-axis (the sense direction).
In this gyroscope, drive oscillation along the x-axis was actuated through the lateral electrodes shown in Figure 8b. Differential capacitive sensing was used in the device.

Optimization Results
The simulation and optimization of the vibratory MEMS gyroscope is based on the system-level model illustrated in Figure 9. The software calculates the drive amplitude of the MEMS device with the drive actuation force governed by Equation (5) and damping (in this mechanism, slide-film damping is the most prominent damping factor). For this actuation mechanism, a DC voltage of 8 V and an AC voltage of 4 V are applied.
Due to driving velocity along the x-axis and the angular rate input, the gyroscope experiences a Coriolis force along the y-axis, resulting in a displacement of the proof-mass in y-axis. This displacement in y-axis causes a detectable capacitance change, described by Equations (4) and (6). For the simulated gyroscope, VDD is defined as 5 V, VCM = 2.5 V and Cint = 100 fF.
For this gyroscope, the first frequency mode is analyzed with different numbers of meshing elements. As observed in Figure 10, the meshing elements reached the optimal number of 9737-after this, for the next 6 data points, the variation in the first frequency mode is less than 0.15%. Thus, the remaining simulation and optimization process will use a number of elements of 9737 as the optimized meshing element size. The GA optimization process is applied to the MEMS gyroscope for 6 generations. The FOM is defined by Equation (13). The software calculates the drive amplitude of the MEMS device with the drive actuation force governed by Equation (5) and damping (in this mechanism, slide-film damping is the most prominent damping factor). For this actuation mechanism, a DC voltage of 8 V and an AC voltage of 4 V are applied.
Due to driving velocity along the x-axis and the angular rate input, the gyroscope experiences a Coriolis force along the y-axis, resulting in a displacement of the proof-mass in y-axis. This displacement in y-axis causes a detectable capacitance change, described by Equations (4) and (6). For the simulated gyroscope, V DD is defined as 5 V, V CM = 2.5 V and C int = 100 fF.
For this gyroscope, the first frequency mode is analyzed with different numbers of meshing elements. As observed in Figure 10, the meshing elements reached the optimal number of 9737-after this, for the next 6 data points, the variation in the first frequency mode is less than 0.15%. Thus, the remaining simulation and optimization process will use a number of elements of 9737 as the optimized meshing element size. electrodes shown in Figure 8b. Differential capacitive sensing was used in the devic

Optimization Results
The simulation and optimization of the vibratory MEMS gyroscope is based system-level model illustrated in Figure 9. The software calculates the drive amplitude of the MEMS device with the dri tuation force governed by Equation (5) and damping (in this mechanism, slide-film d ing is the most prominent damping factor). For this actuation mechanism, a DC v of 8 V and an AC voltage of 4 V are applied.
Due to driving velocity along the x-axis and the angular rate input, the gyro experiences a Coriolis force along the y-axis, resulting in a displacement of the proof in y-axis. This displacement in y-axis causes a detectable capacitance change, describ Equations (4) and (6). For the simulated gyroscope, VDD is defined as 5 V, VCM = 2.5 Cint = 100 fF.
For this gyroscope, the first frequency mode is analyzed with different numb meshing elements. As observed in Figure 10, the meshing elements reached the op number of 9737-after this, for the next 6 data points, the variation in the first freq mode is less than 0.15%. Thus, the remaining simulation and optimization proces use a number of elements of 9737 as the optimized meshing element size. The GA optimization process is applied to the MEMS gyroscope for 6 genera The FOM is defined by Equation (13). The GA optimization process is applied to the MEMS gyroscope for 6 generations. The FOM is defined by Equation (13).
In this equation, Sensitivity is the output voltage when the device is subjected to an angular rate of 1 rad·s −1 in the z-axis, ∆f denotes the difference between drive and sense frequency (influences the mechanical gain in the sense mode to the input angular rate, stated in Equation (14)). Lastly, Q sense represents the sense mode quality factor. This FOM was presented in a multiplication fashion without weighted elements, this is useful when the designer is not sure which performance parameters are the most important ones.
During the six generations, the GA changed the geometric parameters to find the optimal device. The geometric parameters and performance are listed in Table 4. The geometric parameters were the suspension beam's width and length, and sense comb fingers' width. The evolution of the optimization is illustrated in Figure 11. In this equation, Sensitivity is the output voltage when the device is subjected to an angular rate of 1 rad·s −1 in the z-axis, Δf denotes the difference between drive and sense frequency (influences the mechanical gain in the sense mode to the input angular rate, stated in Equation (14)). Lastly, Qsense represents the sense mode quality factor. This FOM was presented in a multiplication fashion without weighted elements, this is useful when the designer is not sure which performance parameters are the most important ones.
During the six generations, the GA changed the geometric parameters to find the optimal device. The geometric parameters and performance are listed in Table 4. The geometric parameters were the suspension beam's width and length, and sense comb fingers' width. The evolution of the optimization is illustrated in Figure 11. Figure 11. Evolution of the MEMS gyroscope through the six generations of the GA. During the optimization, the proof-mass became larger while the u-beams, the sense comb fingers, and the proof-mass frame became thinner. The suspension beam width of the optimal design is much smaller than that of the initial design and the proof-mass size of the optimal design is much larger than that of the initial design. These two geometric changes lead to an increase in sensitivity. Moreover, after 6 generations, the drive and sense mode frequencies as well as the frequency split were significantly reduced. A comparison between frequency modes of the initial and optimized design is shown in Figure 12. The suspension beam width of the optimal design is much smaller than that of the initial design and the proof-mass size of the optimal design is much larger than that of the initial design. These two geometric changes lead to an increase in sensitivity. Moreover, after 6 generations, the drive and sense mode frequencies as well as the frequency split were significantly reduced. A comparison between frequency modes of the initial and optimized design is shown in Figure 12. The mode shapes of the MEMS gyro are illustrated in Figure 12. In the optimized design, the undesired y-axis movement in the drive frame is reduced, while the desired proof-mass displacement is slightly enhanced. The performance of the initial and optimized design is listed in Table 4. The sensitivity increased by 420.9% and the frequency The mode shapes of the MEMS gyro are illustrated in Figure 12. In the optimized design, the undesired y-axis movement in the drive frame is reduced, while the desired proof-mass displacement is slightly enhanced. The performance of the initial and optimized design is listed in Table 4. The sensitivity increased by 420.9% and the frequency mismatch decreased by 56.44%. The FOM was improved by 83.17%. However, the quality factor of the sense mode was decreased by 50.31%, this is explained by the decrease in resonant frequency which is inversely correlated with the quality factor.

Discussion
The software can build geometries of relative complexity, with Pygmsh Python library. Although the most complex geometric operations, such as filet, chamfer, Bezier curve, and arrays are not available, most MEMS inertial sensors can be built by the software through coding. When compared with commercial multiphysics tools, the designed software falls short in its ability to create freeform designs and the lack of graphical interface, which makes it less user-friendly.
The geometry building block passes its parameters and meshes to the FEM part of the program, which can conduct modal analysis and displacement. When compared with COMSOL, only a 0.5% difference in the first frequency mode is reported in Section 3.1.1. COMSOL is a multiphysics software, and so it comprises hundreds of different FEM simulations for different physical domains and phenomena. The software described here, in its current state, can only evaluate displacement and eigenfrequency, but users can easily study different physical domains by inserting the correct equations into the script, referring to the FEnICS documentation [13]. The damping script takes into account the geometric parameters, calculating the squeeze film and slide film damping coefficient, and the quality factor. Finally, the implemented genetic algorithm takes the performance of the devices and proceeds to achieve an electro-mechanical co-optimized design in the end. To demonstrate the effectiveness of the proposed software, two MEMS inertial sensors were designed, simulated, and optimized.
Besides, this software offers the possibility of parallel processing, allowing the computer to process long processes faster and more efficiently.

Conclusions
This study presented a novel electro-mechanical co-optimization methodology for MEMS inertial sensors, entirely based on open source Python. A software comprising geometry design, a FEM simulation, damping calculation, electronic domain calculation, and a GA optimization process, was developed and applied to two MEMS inertial sensors as demonstrators. The developed Python program is a powerful tool that provides designers with limitless customization freedom, presents engineers with complete control of all steps in simulation and optimization and allows efficient management of computational resources according to specific research goals.
The next steps for this software will be the implementation of transient and closed-loop system simulation and optimization. In addition, thermoelastic damping will be added to the damping study block of the system as well as a non-linearity calculation of MEMS devices.
Author Contributions: Investigation, conceptualization, software development, writing-original draft, R.A.E.; conceptualization, methodology, writing-review; C.W.; writing-review, funding acquisition, supervision, M.K. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the INTERREG V-A project "Einstein Telescope EMR Site & Technology" (E-TEST, EMR113).

Conflicts of Interest:
The authors declare no conflict of interest.