Virtual twins of nonlinear vibrating multiphysics microstructures: physics-based versus deep learning-based approaches

Micro-Electro-Mechanical-Systems are complex structures, often involving nonlinearites of geometric and multiphysics nature, that are used as sensors and actuators in countless applications. Starting from full-order representations, we apply deep learning techniques to generate accurate, efficient and real-time reduced order models to be used as virtual twin for the simulation and optimization of higher-level complex systems. We extensively test the reliability of the proposed procedures on micromirrors, arches and gyroscopes, also displaying intricate dynamical evolutions like internal resonances. In particular, we discuss the accuracy of the deep learning technique and its ability to replicate and converge to the invariant manifolds predicted using the recently developed direct parametrization approach that allows extracting the nonlinear normal modes of large finite element models. Finally, by addressing an electromechanical gyroscope, we show that the non-intrusive deep learning approach generalizes easily to complex multiphysics problems


Introduction
Even if the fast development of computational resources has enabled the simulation of complex structures, involving multi-scale and multi-physics phenomena, the ever growing computational costs continuously drives the search of efficient, but accurate, reduced order modeling techniques. Well-known strategies are nowadays routinely applied in linear problems arising from vibratory systems and make use of linear normal modes as an optimal projection basis upon which the equations of motion can be reduced. However, their extensions to problems in nonlinear dynamics is still an open issue, and thus object of intensive research. Early ideas to project the nonlinear equations of motion onto a selected subset of the linear modes basis are progressively being abandoned. Indeed, the presence of strong nonlinear coupling terms between low frequency master modes and high-frequency ones makes the aforementioned approach not efficient, and the initial choice of the trial space is a delicate matter [1]. The Proper Orthogonal Decomposition (POD) method offers a potential gain since it optimally orients the subspace bases to better fit the curvatures of the nonlinear data [2,3,4], and enables the construction of such a subspace from data collected by simulating the behavior of the system for a restricted amount of configurations; however, the use of linear subspaces -possibly of local type -might still represent a computational bottleneck [5]. Other approaches, like the implicit condensation and expansion [6] or modal derivatives methods, have been introduced with the aim of taking the amplitude dependence of modes into account. However, they are limited to moderate transformations and apply only in presence of sufficient slow/fast separation between the slave and master coordinates [7].
On the other hand, truly nonlinear methods for nonlinear vibratory systems are, so far, strictly connected to the concepts of Nonlinear Normal Modes (NNMs) and invariant manifolds, as keys to compute accurate Reduced Order Models (ROMs) [8,9]. These start by defining a nonlinear relationship between the original coordinates and those of the reduced dynamics. In the conservative case, the existence and uniqueness of the searched invariant structures are framed by the Lyapunov centre theorem which guarantees, under non-resonance conditions, the existence of a smooth manifold, tangent at the origin to the associated linear eigenspace, which is invariant. This means that a trajectory of an autonomous system initiated on the manifold will develop on the manifold itself. For dissipative systems, the picture gets more complicated, as the whole phase space is foliated by invariant manifolds tangent at the origin to the linear subspace [10]. The application of these methods to large Finite Element Models (FEMs) has remained sporadic until recently, but is currently receiving an impressive boost by contributions [11,12,13,14] in which direct approaches, called Direct Parametrization of Invariant Manifolds (DPIM), bypass the requirement of computing the whole modal basis. Applications to complex structures with millions of degrees of freedom (DOFs) and featuring also internal resonances and parametric excitation have been recently demonstrated in [13,15]. However, their extension to coupled problems and nonlinearities of generic type is still an open issue, and requires dedicated developments. In parallel, exploiting machine learning methods for constructing realtime surrogates, namely virtual twins or surrogates of nonlinear dynamical systems, has become an area of increasing interest for the system dynamics community. A virtual twin can be considered as a particularly accurate and fast reduced order model of the physical asset, that might be used in the simulation and/or optimization of complex systems in which the asset is allocated or as a key component enabling digital twins.
Even if there is still no consolidated view on the definition of digital twins, they are often interpreted [16,17,18,19,20] as consisting of a physical asset, a virtual representation of that product, and the data connections that feed data from the physical to the virtual representation and vice-versa through mirroring or twinning [21,22,23]. While the operational feedback is naturally provided by the proposed virtual twinthrough the ability of computing quantities and outputs of interest -the twinning stage, realized by means of a physical-to-virtual connection, is not considered in this contribution; nevertheless, it could be profitably included in the process.
The idea of leveraging the modeling capability and learning flexibility of Deep Learning (DL) methods to identify at once nonlinear transformations, nonlinear invariant manifolds and modal dynamics of NNMs from the system response data only, is appealing. Among others, great success has been encountered by the Physics Informed Neural Networks(PINNs) [24], which have been applied in multiple contexts, including solid mechanics [25]. Other DL-based ROM techniques have come as an inspiration to handle the complex reduction process of dynamical systems, unveiling low-dimensional features from black-box data streams [26,27,28]. A peculiar perspective has been taken by Fresca et al. [29,30], who proposed the non-intrusive POD-DL-ROM technique to address the fast simulation of parametrized differential problems. POD-DL-ROMs suitably combine a preliminary POD projection of snapshots, a Convolutional Autoencoder (AE) and a Deep Feedforward Neural Network (DFNN) to enable the construction of an efficient ROM, whose dimension is as close as possible to the number of parameters upon which the solution of the differential problem depends. The encoder part performs an operation of feature extraction forcing the high dimensional data to be reduced, at the bottleneck layer, to few reduced variables. We highlight that this approach builds on the idea that the system dynamics develops on a low dimensional invariant manifold, setting a clear parallel with the DPIM. Compared to other surrogate models exploiting machine/deep learning algorithms, a distinguishing feature of DL-ROMs is their capability to compute, at testing time, the whole solution field, for any new parameter instance and time instant, thus enabling the extremely efficient evaluation of any output quantity of interest depending on the solution field, and generating a truly virtual twin of the structure analysed.
In [31], an extension of the DL-ROM technique has been proposed and applied to microstructures, showing very good predictive capabilities. In particular, a first dimensionality reduction of the data is achieved by means of a POD-Galerkin (POD-G) ROM in order to reduce the costly FOM data generation phase, thus defining the POD-G DL-ROM technique. An alternative approach is represented by the SINDy method (Sparse Identification of Dynamical Systems), a regression technique for extracting dynamics from time-series data [32,33,34]. In the past few years, SINDy has been widely applied to identify models of fluid flows, convection phenomena, structural modelling, and many others. The main appeal of the SINDy approach is to generate, by combining a set of pre-defined (analytical) functions collected in a suitable dictionary, an explicit ROM in the form of first-order differential equations, which can be then manipulated with standard numerical tools. Even if applications have been limited so far to rather small problems, very recent works have fostered the use of SINDy in combination with AE neural networks [35] and preliminary POD reduction. The SINDy approach, with these recent extensions, and the aforementioned POD-DL-ROM, are thus converging to a similar framework, the difference being the way in which the reduced dynamics is integrated and approximated. Convolutional or deep AEs thus emerge as distinctive and powerful tools for dimensionality reduction of dynamical systems.
However, a certain sense of distrust still pervades the community of nonlinear dynamics, and computational mechanics in general, towards these emerging DL approaches which are seen as black boxes lacking reliable foundations. Nevertheless, the strict connection with the DPIM emerges as quite evident, as discussed for instance in [36,37]. In the latter contribution, a data-driven identification of NNMs via physics-integrated DL is proposed; in this case, the architecture integrates prior physics knowledge of the NNMs by embedding physics-based constraints. Intrinsic coordinates are required to have modal properties of a dynamical system, including statistical uncorrelations between desired modal displacements and modal velocities, and potential sparsity of dynamics in the intrinsic space. Their ideas are tested on systems having 2 and 3 DOFs.
The POD-DL-ROM can be considered as the data-driven counterpart of the DPIM. Even if the two techniques arise from different perspectives, they share several analogies, thus generating a sort of duality. Indeed, they both rely on the use of a Full Order Model (FOM) -typically represented by the FEM -and they try to extract salient features of the problem being analyzed by reducing the model to an intrinsic space with very limited dimension. The reduced variables govern the evolution of the dynamics on invariant manifolds. This reduction (or encoding) phase is performed offline with dedicated software and hardware as it can be memory and CPU intensive. A second shared phase is the integration of the reduced dynamics, which is fast and can be carried out almost in real-time. While DPIM achieves this latter task by integrating the reduced differential equations, e.g. with continuation techniques, POD-DL-ROMs adopts a DFNN, which directly approximates the parameter-reduced solution map. Finally, the whole field of FEM nodal displacements of the original FOM can be reconstructed (decoded) with great accuracy. For instance, the ROMs generated with both approaches investigated herein can be queried for the stress field at locations that do not need to be specified a priori, provided the query deals with parameter values falling within the range of (or not too far from) the training parameter set.
Some features however differentiate the two approaches. The DPIM is an almost exact procedure which, for a given FOM with fixed parameters and forcing type, generates a ROM which is highly predictive for all initial conditions and forcing levels. However, it does not generalise automatically to other physics and nonlinearities. Conversely, the DL-based ROM accommodates a whole range of parameters selected by the user in the training phase, can be extended almost automatically to coupled physics, and has real-time performances. However, it still suffers from limited predictive ability and does not guarantee the same level of accuracy as the DPIM. Indeed, such a ROM strategy account for the problem physics only through the snapshots data, employed to train the neural network.
As a natural consequence of these remarks, the aim of this contribution is twofold. After a short introduction to the governing equations in Section 2 and of the DL-based ROM method in Section 3, we further improve the POD-DL-ROM technique [38,29] in Section 3.1 through the use of an arc-length abscissa defined on the Frequency Response Functions (FRFs) to enable the analysis of complex interaction problems in real nonlinear microstructures showing, e.g., internal resonances and coupled electromechanical interactions inducing autoparametric effects. In parallel, we explore in-depth the relationship between the DL-ROM and the DPIM approaches in Section 4. We devote particular attention to the convergence of the POD-DL-ROM with respect to the dimension of the reduced space, analyzing it in terms of both the FRF and the Invariant Manifolds predicted. The applications are proposed in Section 5. In the first two examples addressed, which focus on mechanical geometric nonlinearities, the DPIM simulations are used as a reference solution due to their known accuracy, and the convergence of the POD-DL-ROM approximation to the DPIM one is investigated. In the final example, on the contrary, in order to show that the POD-DL-ROM easily generalises to complex multi-physics without losing its non-intrusive nature, we exploit a commercial code to generate snapshots of an electromechanical gyroscope, aiming at simulating its oscillating behavior in real-time.

Problem formulation
We will focus on mechanical structures subjected to periodic forcing that undergo large transformations but still experience only small strains, a condition which is well described by the Saint Venant-Kirchhoff constitutive model. As detailed in [31], the space discretization of the governing equations by means of finite elements yields a system of coupled nonlinear differential equations representing the Full Order Model (FOM): where the vector U(t) ∈ R N h collects the N h unknown displacements nodal values; M ∈ R N h ×N h is the mass matrix; C = (ω 0 /Q)M is the Rayleigh mass-proportional damping matrix with ω 0 reference eigenfrequency and Q quality factor; F(t; µ) ∈ R N h is the nodal force vector which depends on the vector of n µ parameters µ ∈ P ⊂ R nµ and expresses the actuation due, in general, to a multiphysics coupling e.g. with piezolectricity or electrostatics. In particular, the vector F(t; µ) depends on the angular frequency ω of the actuation and is nonlinearly modulated in amplitude by the coefficient β.
The internal force vector has been exactly decomposed in linear, quadratic, and cubic power terms of the displacement: K ∈ R N h ×N h is the stiffness matrix related to the linearized system, while G ∈ R N h and H ∈ R N h are vectors given by monomials of second and third order, respectively. We stress that the components of these vectors can be expressed using an indicial notation as Eqs. (1) represent our high-fidelity FOM which depends on the input parameters µ. Our goal is the efficient numerical approximation of the solution manifold: The numerical solution of Eqs. (1) to compute the steady state response is a challenge by itself for large scale problems. One option is the use of time marching methods to simulate a sufficiently large number N c of cycles, where N c is typically inversely proportional to the damping. This technique resorts to robust algorithms implemented in most of the commercial software, but when damping is very small, as in most MEMS devices, the computational effort may be unaffordable. In other approaches, like the Harmonic Balance (HB) method [39,40], the unknown displacements are expressed as the sum of Fourier components thus automatically respecting periodicity conditions. However, their implementation requires dedicated codes and non-standard computing facilities.

Data-driven Reduced Order Modeling through DL-ROMs
In this section, we briefly review the construction of deep learning-based ROM (DL-ROM) technique and its extension to the POD-enhanced version (POD-DL-ROM).
The DL-ROM technique [29] is extremely efficient and is able to model highly nonlinear problems by identifying the manifold underlying the dynamics of the system in a complete data-driven and black-box, non-intrusive way. On the other hand, its data-driven nature implies that the sampled data provided must span the parameter space of interest and contain all the information necessary to accurately approximate the solution manifold. In this way, the size of the training dataset increases with the number of parameters considered in the FOM.
DL-ROMs aim at approximating the map (t, µ) → U(t, µ) by describing both the trial manifold, approximating (2), and the reduced dynamics through deep neural networks, which are trained on a set of FOM snapshots. In particular: • to map the FOM solutions in a low-dimensional coordinates vector representation (encoding stage), we use the encoder function of a Convolutive AutoEncoder (CAE) where f E n (·; θ E ) : R N h → R n is obtained as the composition of several layers (some of which are convolutional), depending upon a vector θ E of parameters; • to describe the system dynamics (reduced dynamics learning), the intrinsic coordinates of the ROM approximation are defined as where φ DF n (·; ·, θ DF ) : R (nµ+1) → R n is a DFNN, consisting in the subsequent composition of a nonlinear activation function, applied to a linear transformation of the input, multiple times [41]. Here θ DF denotes the vector of parameters of the DFNN, collecting all the weights and biases of each layer of the DFNN; • to model the reduced nonlinear trial manifold S n h ≈ S (decoding stage) we employ the decoder function of a CAE [42,43], that is, where f D (·; θ D ) : R n → R N h depends upon a vector θ D collecting all the corresponding weights and biases.
The DL-ROM approximationŨ(t; µ) ≈ U(t; µ) is then given bỹ and is computing by solving a suitable optimization problem, in the variable θ = (θ E , θ DF , θ D ) (see, e.g., [29] for further details). The POD-DL-ROM technique [30] consists in applying the DL-ROM technique to the intrinsic coordinates of a linear trial manifold generated through randomized singular value decomposition (rSVD) and approximating S; In particular, to reduce the dimensionality of the snapshots and avoid feeding training data of very large dimension N h , POD is first applied -realized through randomized SVD (rSVD) [44] -to the snapshot set; then, a DL-ROM is built to approximate the map between (t, µ) and the POD generalized coordinates u N (t; µ) = V T N U(t; µ). Note that, in the POD-DL-ROM framework, the encoding stage results from the projection of FOM snapshots onto the linear trial manifold defined by V N and the compression performed by the encoder function of the CAE. Similarly, the decoding stage is formed by the application of the decoder function and the multiplication for the POD basis matrix. The architecture of the POD-DL-ROM neural network, employed at training time, is the one shown in Figure 1; note that, at testing time, as in the DL-ROM technique, we can discard the encoder function.
A variant of the POD-DL-ROM approach, called POD-G DL-ROM, has been recently benchmarked on several MEMS structures including beams, arches and mirrors in [5] with the aim of reducing the cost of the training phase. A limited number of FOM analyses are used to identify an optimal linear basis via SVD, and the ROM obtained projecting the governing equation on this basis is then employed to train the neural networks. However, to assess the capacity of the POD-DL-ROM to tackle other, more involved physical problems, the initial formulation is adopted in this work.

Arc-length on the FRF as control parameter
In applications to microstructures, one of the most meaningful outputs of the analysis is the Frequency Response Function (FRF) which defines the steady-state periodic response of the dynamical system. In a FRF, a selected output like, e.g., the maximum midspan deflection of a beam, or the rotation amplitude of a micromirror, is plotted as a function of the actuation intensity and frequency. In linear systems, the FRF is a single-valued function of the forcing frequency ω that is hence utilized as the ordering parameter also for the snapshot matrix. When nonlinearities are present, this property is generally lost. For instance, in systems behaving as simple duffing oscillators with hardening (or softening) properties, the phase of the response with respect to the forcing signal has been exploited in [31] as order parameter.
In more involved applications, even the phase does not ensure uniqueness of the response, and a more general approach is needed. Inspired by the continuation technique used to integrate the ROM in the DPIM approaches, we opt for the use of the arc-length abscissa along the FRFs. Let us consider for instance the typical FRF for a 1:2 internal resonance of a shallow arch, as illustrated in Figure 2a), where the curves for two different forcing levels are plotted. In order to express both the frequency and the amplitude as single valued functions, we first compute the arc-length abscissa along the FRF. However, the total length of the FRFs strongly depends on the forcing levels; moreover, the introduction of this abscissa induces a misalignment of physically relevant phenomena, like e.g. resonances and mode couplings. These are, on the contrary, naturally clustered at very similar values of the frequency and phase with positive fall-outs in the DL-ROM.
As a solution to this issue, we introduce a piecewise normalization of the arc-length abscissa that synchronizes maxima and minima of the FRFs. Four regions are identified, as illustrated in Figure 2, and the arc-length of each region is rescaled from 0 to 1, yielding a total arc-length of 4. Even if this procedure requires some insight into the physics of the response, it is however very general and powerful, thus representing a major improvement with respect to other procedures which would not be applicable in this context [31].
We stress moreover that the arc-length abscissa is not only used as control parameter in the snapshots passed to the DL-ROM; indeed, it also easily allows to retrieve, a posteriori, the FRF frequency-amplitude relationship through interpolation on the training data.

Direct Parametrization of Invariant Manifolds
The parametrization method of invariant manifolds, introduced by Cabré, Fontich and de la Llave in [45,46,47], has been very recently applied to large FEM models of nonlinear vibrating structures [11,12,15].
In particular, in [15] it is shown that lightly damped and mildly forced systems can be parametrized with very good accuracy by considering only the autonomous part of the equations and adding a posteriori the forcing by projecting F in Eq.(1) on the linear master mode. This is the approach retained in the present work. The general idea is to reduce the dynamics to an invariant manifold tangent to the eigenvectors selected as master modes. Let us assume that n master coordinates are selected, with n N h . These master coordinates are linked to their corresponding vibration modes and the searched invariant manifold is the nonlinear continuation of the subspace spanned by the n vibration modes. Working in the phase space in terms of both displacements and velocities, the invariant manifold becomes 2n-dimensional; to describe the reduced dynamics on this manifold, 2n normal coordinates z are introduced. The 2N original coordinates U (nodal displacements) and V (nodal velocities) are then expressed as a function of the new normal coordinates z as where the two nonlinear mapping functions Ψ(z) and Υ(z) are time-independent unknowns to be computed. The reduced dynamics of the autonomous system governs the evolution onto the corresponding invariant manifold. This latter is also unknown at this stage, and can be expressed as: The aim of the method is to compute Ψ, Υ and f . The procedure consists in differentiating Eq. (7) with respect to timeU and substituting in the equations of motion (1) -suitably rewritten as first-order system -to obtain the so called invariance equations These nonlinear equations can be solved using asymptotic expansions in the unknowns (the reduced coordinates z), as proposed in [47,48]. The procedure for large FEM models has been detailed in [11,12], extended to high-order expansions in [13] and to forced systems in [15], where it has been extensively validated against FOM HB solutions. An open-source version of the code has been published in [49] and is now a reference for geometrical nonlinearities with an impressive list of successful applications ranging from internal resonance, autoparametric excitation and folding of the manifolds due to extremely large transformations, just to mention a few. However, the DPIM formulation has not been extended yet to coupled multiphysics. One of the greatest benefits of the DPIM is the generation of a ROM in the form of an explicit, smallsize system of first-order differential equations -see Eq.(8) -which can be solved with different approaches. MEMS structures, and in general periodically excited systems, are characterized by their steady state regime. Since damping is typically very low, time-marching methods are less appealing than Harmonic Balance, Collocation or Shooting techniques that allow computing the steady state response directly. Furthermore, these methods are compatible with continuation approaches that allow reconstructing the whole FRF, with both stable and unstable branches. Several packages suitable for small-scale problems like the one given by DPIM are freely available. Among them we cite Auto07p [50], Manlab [51,52]; Nlvib [53], MATCONT [54], COCO [55] and BifurcationKit [56], an emerging toolkit for Julia language. All the continuation solutions discussed in this work have been obtained using MATCONT.

Duality with the DL-ROM in a simple case
In order to better highlight the striking duality/analogy between the DL-ROM and the DPIM, we now focus on the real formulation of the DPIM proposed in [11], and restrict the attention to a single mode reduction in a mechanical system with asymptotic expansion truncated at second order. The basic ingredients of the procedure can be expressed in terms of the real coordinates R and S, i.e. the reduced normal coordinates, which are the counterpart of z in Eqs. (7) and (8). In this case, the nonlinear change of variables reads: and allows expressing the FOM nodal values for displacements U and velocities V from the normal coordinates. The reduced model can be formulated in terms of two first-order equations as follows: All the vectors and coefficients in Eqs.(11)- (14) are computed in the parametrization procedure while β denotes the forcing level. Let us now underline some contact points between the DL-ROM and the DPIM (see Figure 3): 1. All the coefficients and vectors in Eqs.(13)- (14) are computed offline starting from the FOM, in a preliminary phase, that can be seen as the equivalent of the encoding phase in the DL-ROM. The costly offline training can be performed on dedicated platforms and software. Both approaches hence are based on a FOM, which is typically built exploiting a finite element discretization.
2. The ROM Eqs.(13)- (14) emerges as the counterpart of the DFNN in Eq.(4) and can be integrated online with almost real-time performance when the model is queried with specific values of the forcing parameters.
3. Eqs.(11)-(12) represent the parallel of the decoding phase which reconstructs the global nodal fields starting from the online integration of the ROM. Thus both ROMs can reproduce the same richness in details of the original FOM, since the decoding phase allows to generate a full field information.
It is worth recalling that the SINDy approach [32,33,34] has been recently coupled with autoencoders for order reduction [35] and shares several of the similarities in common between the DPIM and the DL-ROM. However, differently from the latter, SINDy expresses the reduced dynamics in the form of an ODE system whose coefficients are estimated by means of sparse identification procedures. Figure 3: Schematic comparison between DPIM and DL methods. One can identify a pre-processing stage, an encoding stage, a reduced dynamics solution stage and, finally, a decoding stage. For DL methods, in this scheme we distinguish between black box approaches for the reduced order dynamics, as in the DL-ROM [29], and model based approaches, as in SINDy proposed by Brunton et al. [33] 5 Applications In this section, we discuss three different applications: a micromirror, a shallow arch showing 1:2 internal resonance, and an electromechanical gyroscope. In the first two examples addressed, which focus on geometrical nonlinearities, the DPIM simulations are used as reference solutions, and the convergence of the DL-ROM to the DPIM is investigated. A detailed discussion on the performance of the DPIM in these kind of applications has been developed in [12,13] through extensive validations against the HB approach.
In the final example, on the contrary, we exploit a commercial code to generate FOM snapshots of an electromechanical gyroscope in order to show that the construction of a DL-ROM can be easily generalized to the case of complex multi-physic problems, still retaining its non-intrusive peculiar character. In this third example, validation will be performed against the results of the commercial FOM itself.
Reconstruction of the whole response As highlighted in the previous Sections, the methods analysed are not limited to the evaluation of the FRFs of selected quantities like, e.g., the rotation angle or the midspan displacement, but they accurately regenerate, through the decoding phase, an approximationŨ of the whole displacement field of the original FOM. Next,Ũ can be post-processed to generate all the desired outputs, a typical example being the monitoring the evolution of stresses in the structure to check the admissibility of the given design [31]. In our case, we will useŨ to reconstruct the FRFs for the master and the slave modes, in order to investigate the convergence of the DL-ROM to the DPIM results.
Given the estimate of the nodal displacement vectorŨ at any instant t, we define a generalized modal coordinate u i (t) and modal amplitude A i according to: (|u i (t)|) (15) in which M is the mass matrix, Φ i is the i-th eigenmode and T is the period. Modal coordinates are also used to generate error measures for each mode: where u P i denotes the DPIM solution, u DL i the DL-ROM one, and (u P ) 2 =Ũ T MŨ. The proposed errors are essentially time averages of the instantaneous error of the DL-ROM with respect to the DPIM reference. The former one, E r i , is a relative indicator for a given mode, but it does not account for the absolute importance of the mode itself in the global response. The latter one, E i , instead normalizes the error with respect to a global displacement measure.
Computation of manifolds The underlying key idea and assumption of the two reduction techniques addressed in this contribution is the existence of invariant manifolds on which the system response evolves as a function of the (few) master coordinates. As a consequence, in all the examples we will analyse the existence of DL-ROM manifolds and their convergence, when applicable, to the DPIM ones. While plotting the manifold for the DPIM is straightforward as both displacements and velocities are directly accessible, see Eq. (7), a DL-ROM provides, on the contrary, only a reconstruction of displacements. In order to generate the manifold, the modal velocity v DL i (t) can be computed resorting to a Fourier decomposition of the periodic u DL i (t), and then differentiating it with respect to time. Indeed, this approach allows to smooth high-order components through harmonic truncation, and is preferred, e.g., over finite differences that might generate noisy derivatives.
Because of the inertial and geometrical effects triggered by large rotations, micromirrors are intrinsically nonlinear, and the prediction of their dynamic behaviour is essential to guarantee a proper design and control of the device during operations. Recently, the authors have developed a large scale HB-approach in [39] for the analysis of piezo-actuated mirrors. However, the computational cost entailed by this approach limits its applicability during the design phase and for online monitoring, ultimately stimulating a frantic search for efficient ROMs. This example is a tough challenge not only for classical implicit condensation approaches [6] that cannot deal with large rotations, but also for the most advanced and recent techniques, like the DPIM approach. The torsional mode is not the lowest-frequency one and is not well separated from other eigenmodes, entailing the failure of the quadratic formulation of the DPIM utilized in [12], therefore requiring a high-order expansion, as remarked by Vizzaccaro et al. [13]. In this section, we will consider the MEMS micromirror illustrated in Figure 4, fabricated by ST Microelectronics. The mirror is assumed to be made of isotropic polysilicon [60], with density ρ = 2330 Kg/m 3 , Young modulus E = 167 GPa and Poisson coefficient ν = 0.22. The mirror plate, a circle of diameter 1970 µm , is suspended to a gimbal connected with a torsional beam along the rotation axis and two suspension beams on each side. Thanks to symmetry, only half of the device is modelled with the FEM, involving in total 9732 dofs; Dirichlet boundary conditions  Figure 4c) and on the symmetry plane. On the remaining boundaries, zero traction Neumann boundary conditions are instead imposed. The mirror has been actuated with a fictitious body force proportional to the inertia forces of the master mode, F(t) = Mφ 3 β cos(ωt), in order to enable the comparison with the DPIM approach [13]. The first four eigenmodes are illustrated in Figure 5, where also the eigenfrequencies obtained from a linear FEM eigenvalue analysis are listed. The torsional master mode, of interest herein, is the third one and has a frequency of 29271 Hz. A quality factor Q = 1000 is considered in all the analyses.
Since the central plate is stiff, we adopt its angle of rotation θ as reference output for the FRFs in Figure 6. It should be noted that the amplitude-induced shift of the peak frequency is small with respect to the absolute value. However, deviations of few Hertzs in working conditions affect the optical performance of projectors; hence, the first requirement for the model is to be highly predictive even of tiny deviations.
This example, that has already been discussed in [31], is here revisited using the new arc-length abscissa which corrects the artifacts that were visible at the edges of the frequency range, and focusing on the convergence of the slave modes and manifolds to the DPIM reference ones. The DPIM is also used to compute the training data for the DL-ROM method. In the following, we consider 63434 snapshots providing 161 samples on each period; the parameter space spanned is (β, s) in {1, 1.5, 2, 2.5, 3}µN ×[0 : 2.0]. The arc-length abscissa s is rescaled in this example between 0 and 2, and the alignment of the data is forced only at peaks of the FRFs, i.e. at s = 1. One master mode is used in the DPIM analyses, which implies that the dimension of the reduced space is 2 and the reduced variables can be interpreted as the generalised displacement and velocity of the torsional mode. The DL-ROM FRFs for the rotation angle are collected in Figure 6 and are very accurate even for p = 1, i.e, with just one reduced variable. Actually, the difference in FRFs obtained with increasing p can hardly be appreciated. The DL-ROM is thus extremely efficient in extracting the key features of the response of this lightly damped device. Indeed, as a first approximation, the mirror can be considered as a conservative system which is well described by a single displacement-like parameter.
The accuracy of the DL-ROM is further investigated in Figure 7 which presents the error measures, Eqs.(16)- (17), and the FRFs for the first four modal coordinates u i , Eq.(15), as a function of the forcing frequency. Inspecting the error plots it can be appreciated that, while the master mode is unaffected by the additional reduced variables, these latter on the contrary improve the representation of slave modes. In particular, the low frequency ones benefit from the enrichment of the reduced space, while high-order modes are less sensitive. It should be stressed, anyway, that even if the relative error on these modes remains large, their contribution to the global displacement field is negligible, as illustrated in Figures 7b). Moreover, the inclusion of a third reduced variable has a minimal impact on the response, in accordance with the DPIM that perfectly matches the HB FOM with only a single master mode, i.e., two reduced variables [13]. Furthermore, the DL-ROM training process is subjected to a certain degree of stochasticity also because of the stochastic gradient descent method employed during the training of neural networks and consequently the results obtained with p = 2 seem to be, in some conditions, slightly better than the one obtained with p = 3. All these remarks are confirmed by the FRFs of the modal coordinates in Figure 7. Indeed, the FRF in Figure 7e), referred to the master mode, is already at convergence with p = 1, while those in Figure 7c,d,f), plotting the amplitude of the low frequency slave modes, benefit from the introduction of a second reduced variable and are almost insensitive to the third one. In order to further gain insight into the DL-ROM performance, we inspect the invariant manifolds of Figure 8, computed as explained in the introduction to Section 5. Each manifold plots the modal coordinate u i as a function of u 3 , the coordinate of the master torsional mode, and of its time derivative v 3 . Two low frequency slave modes are considered, i.e. modes 1 and 4 in Figure 5. The DPIM invariant manifolds are represented as smooth grey surfaces, while DL-ROM orbits are plotted as continuous lines of different colors. The orbits obtained with p = 1 are plotted separately in the Figures on the left, in order to better stress the fact that they correspond to single concavity surfaces which are, as expected, velocity independent. Indeed, having one single reduced variable, the DL-ROM cannot introduce a dependence on its time derivative. However, it emerges that the DL-ROM selects, among all options, the single curvature manifolds that better interpolates the correct one. Manifolds of slave modes match accurately the DPIM ones already with 2 reduced variables, as clearly illustrated in Figure 7.

Internal resonance in a shallow arch
Internal resonances (IRs) are associated to energy transfer between modes, and are frequently experienced by MEMS structures mainly due to the very low dissipation. Often IRs are strongly linked to the stability of the associated periodic response, and quasi-periodic regimes might arise as a consequence of Neimark-Sacker (NS) bifurcations [61]. The numerical prediction of such phenomena has been tackled only recently with ad-hoc approaches like the Implicit Condensation [6] or the DPIM [62]. Such a task requires a stability analysis which can be run on small ROMs using dedicated continuation tools, and cannot be performed using FOMs, in general.

Figure 10: Shallow arch. Geometry and eigenmodes
A typical structure that is prone to 1:2 IR is the shallow arch illustrated in Figure 10. The layout, inspired by the one proposed for a bistable structure in [63], has been suitably calibrated so as to trigger the desired IR. A MEMS resonator presenting the similar geometry and IR has been studied in [64]. The geometry and the mesh employed for the FOM are illustrated in Figure 10. The arch dimensions are B = 20 µm, H = 5 µm, L = 530 µm, rise=13.4 µm, s = 10 µm. The mesh consists of quadratic wedge elements and contains 1971 nodes. The device is made of polycrystalline silicon with density ρ = 2330 kg/m 3 and a linear elastic Saint-Venant Kirchhoff constitutive model is assumed, with Young modulus E = 167000 MPa and Poisson coefficient ν = 0.22 [65]. The lowest eigenfrequencies of the modelled structure are reported in Figure 10. The quality factor has been set to Q = 500 and the actuation is provided by a body force proportional to the first eigenmode F(t) = Mφ 1 β cos(ωt) with β load multiplier.
The 1:2 IR occurs as a result of the interaction of the first and the sixth eigenmodes having frequency ratio close to 2 i.e. satisfying the necessary condition to induce a 1:2 internal resonance, and leads to qualitative and quantitative changes in the dynamics. An in-depth analysis has been developed in Gobat et al. [61], where 1:2 IR systems are analysed starting their normal form and the existence of the so-called parabolic modes between the coupled oscillators is demonstrated. Such modes exist on the two branches of the FRF (see Figure 11) and are associated to the two backbones that onset from ω 1 and ω 6 /2, respectively. In order to generate a ROM with the DPIM, the two interacting modes must be both selected as master modes and the dimension of the reduced space is 4. A critical analysis of the convergence of the DPIM to the full-order HB simulations has been presented in [62], and on this basis the DPIM results are here considered as an exact reference. Addressing this problem with the DL-ROM is a very ambitious task that brings this Deep Learning-based approach to a new level of complexity, as the autoencoder has to automatically recognise the transition between totally different pictures of the dynamical response. Moreover, this example is used to further investigate the impact of the dimensions of the reduced space for the autoencoder.
Following the same path already considered for the micromirror, in Figure 11 we first plot the FRFs for the midspan displacement d, which is the typical output for this kind of structure. The reference curves have been computed with the DPIM; note that the DL-ROM FRFs obtained with 2 and 3 reduced variables are almost superposed, thus showing its convergence, as opposed to the curve obtained with p = 1. Indeed, due to the presence of an interaction between two modes, a single reduced variable in the DL-ROM cannot describe the full evolution even for the master modes.
Next, Figure 12 presents the errors, Eqs.(16)- (17), and the FRFs for four selected modal amplitudes, Eq.(15), as a function of the forcing frequency for p = 1 : 5. Inspecting the error plots and the FRFs 12c)-d) for the master modal amplitudes A 1 and A 6 , it can be appreciated that these are already at convergence with p = 2, consistently with what remarked for the mirror, where only one master mode was active and p = 1 provided accurate FRFs for the master amplitude.
Also in this case, increasing the number of reduced variables has an impact only on the representation of slave modes. The check is performed on modes 11 and 13 in Figure 12, where p = 1 : 5. Convergence to the DPIM develops up to p = 4, while the inclusion of a fifth reduced variable has a minimal impact on the response, in accordance with the DPIM that perfectly matches the HB FOM with two modes, i.e. four reduced variables [62]. It should be remarked that these slave amplitudes are at least one order of magnitude smaller than the master ones and, as a consequence, larger errors remain.
In Figure 14, we start comparing the 3D trajectories in the space u 6 , u 1 , v 1 , i.e. we plot the modal displacement for the second master mode in terms of the coordinates of the first master. In this space, the interactions between the two modes leads to foldings and intersections of the orbits envelop as the whole space tends to be filled with orbits. To tackle these difficulties we will focus on a specific forcing level generating the FRF of Figure 13a); furthermore, we will consider five different portions of the FRF denoted by different colors. The same colors are used to represent the envelops obtained with the DPIM in Figure  13b). In order to better investigate the agreement between DPIM and DL-ROM, the different portions are considered separately in Figures 14a)-d), where the smooth grey surfaces are the DPIM results, the DL-ROM orbits are traced as lines and the colored dots recall which region of the FRF is being investigated. From these surfaces, we can appreciate the aforementioned transformation of the manifold which is dominated by a quadratic term turning from negative to positive [61]. Despite the complicated relationship between the orbits on each mode, the DL-ROM approach provides an almost perfect agreement. The curvatures are well reproduced both qualitatively and quantitatively. Since we are considering in Figure 14 only master modes, this result confirms our understanding, previously underlined, that the DL-ROM identifies the most critical features of the response and represents them as accurately as possible: two reduced variables are sufficient to represent their evolution of two lightly damped, almost conservative modes.
The inspection of the invariant manifolds for slave modes is indeed a much more involved task than for the mirror case, as here the manifolds should be plotted as a function of the four variables u 1 , v 1 , u 6 , v 6 . We thus limit ourselves to tracing envelops of slave trajectories for the modal displacement u 13 in terms of u 1 , v 1 , illustrated in Figure 15. The comparison is done for the same forcing level as before, and the five different regions of the FRF. The smallness of the u 13 coordinate makes the comparison difficult, as previously recalled for the FRFs; however, the convergence of the DL-ROM trajectories to the DPIM smooth surfaces is evident.

Electromechanical Disk Resonating Gyroscope
In the previous two examples, an actuation through fictitious body forces has been considered and, thanks to this choice, an extensive validation with the exact DPIM approach has been proposed. However, the DPIM has not been applied to multi-physics problems yet, as each extension requires dedicated and costly developments and coding. Therefore, this application represents a challenging test-bed for the DL-ROM technique, that so far has only been applied to a multi-physics problem in [66], where a fluid-structure interaction problem depending on a set of physical parameters has been addressed. We aim at showing the great versatility of the DL-ROM approach, as well as the possibility to easily generalize its construction to face a complex multi-physics problems. Another positive feature, useful in this context, is the non-intrusiveness of DL-ROMs, as snapshots can be readily generated with commercial codes. These capabilities, that are crucial in proposing realistic applications to microstructures, are underlined in this section by addressing an electromechanical Disk Resonanting Gyroscope (DRG) taken from the library of Coventor MEMS+ TM [67], a leading tool for the analysis of MEMS. The original model available in the software, inspired by the device proposed in [68], has been slightly modified in our benchmarks, and all the details are provided in what follows.
The DRG external ring is modelled with 32 Eulero-Bernulli beams, and each arch suspension is made of 2 Eulero-Bernulli beams. The center support is a rigid body constrained to the ground. A series of parallel-plate electrodes are placed around the ring and electrostatic forces are modelled by the software by means of conformal mapping. We highlight that, compared to the element size and geometry, the expected displacements are small, so that we neglect geometric nonlinearities and all the nonlinear effects are induced Figure 17: DRG: first linear eigenmodes. Due to the symmetry properties, many of the eigenmodes come in degenerate couples sharing the same eigenfrequency and differing only by a rotation around the z axis. Modes 7 and 8 are, respectively, the drive and sense modes.
by the electromechanical coupling.
In Figure 16, the gyroscope components are coloured according to their function. This device allows detecting angular velocities around the z axis by exploiting the two degenerate modes 7 and 8 in Figure 17 which are characterized by radial displacements of the outer ring proportional to cos(2θ) and sin(2θ) with θ polar coordinate running on the external ring. During operations, the drive mode, i.e. mode 7, is excited imposing the bias V AC sin ωt to the blue electrodes and the bias −V AC sin ωt to the yellow ones. The V AC values range from 2 mV up to 10 mV guaranteeing an in-plane displacement of approximately 0.54 µm . A constant potential bias V DC = 1 V is imposed to the gyro (red structure). All the remaining electrodes are grounded. When the device is subjected to an external angular velocity rate ω z , the Coriolis effect exerts in-plane forces on the ring and activates the sense mode (mode 8 in Figure 17) that is detected by the two couples of independent parallel-plate electrodes (violet and grey) rotated by 45 degrees with respect to the drive ones.
However, it has been shown (e.g. in [69]) that the response of the sense mode largely exceeds predictions based on linear models. Indeed in the equation of motion for the sense mode, u 8 , a term proportional to DC sin 2 ωt appears [70], i.e. a linear stiffness term modulated in time at twice the frequency of the mode. This is responsible for a dynamical instability, called autoparametric resonance [71,72], which rather abruptly triggers the activation of the sense mode (see Figure 18d). Simulating the onset of the autoparametric resonance is a challenge for any numerical tool as it induces rapid changes in the overall response, migrating from a single master mode evolution towards a two master modes dynamics. For this reason, we assume ω z = 0 in what follows and focus on autoparametric effects.
The DL-ROM has been trained using 149 900 snapshots representing 100 points per period on the parameter space f = [32.620, 32.6498] kHz, V AC = [2,3,4,5,6,7,8,9,10] mV. We highlight that the training set is larger than in the previous examples mainly due to the changes introduced by the parametric resonance. At the onset of autoparametric effects the drive mode experiences a saturation, while the sense mode abruptly reaches comparable amplitudes [73,74]. The results on a testing set are collected in Figure 18 for different choices of number of reduced variables p. The comparison is performed taking the Coventor results as reference.
The FRF of the drive amplitude d 1 , represented by the radial displacement of the node indicated by the red circle in Figures b), is plotted in Figure 18a). The curves below a given V AC denote a simple softening response of the drive mode, due to the electrostatic nonlinearities, while a plateau starts developing when the sense mode gets autoparametrically activated. Inline with the previous examples, in the presence of two master modes one reduced variable, p = 1, cannot represent adequately the main features and p = 2 is required. No measurable benefits can be observed when further increasing p. Similar remarks also hold for the FRF of the spurious amplitude d 2 plotted in Figure 19d)-f).
Also the error estimates 19a)-b) show a convergence trend with the increment of p and confirm that additional reduced variables only bring benefits to slave modes. It is however worth stressing that, in this particular example, slave modes reach only a negligible amplitudes, within the numerical errors of the procedure, which prevents from presenting meaningful convergence analyses.
Finally, in Figure 20 we compare the envelop of orbits in the space u 8 , u 7 , v 7 . As for the shallow arch example, these master mode orbits tend to fill the whole space. Hence, for the sake of clarity, we focus on the single FRFs illustrated in Figure 20a),b) which are partitioned in three different regions generating the Coventor MEMS+ reference envelops of Figure 20c). The comparison is performed separately in Figures 20d)-f) for the three different portions highlighted. Again, p = 1 is omitted, but the match is excellent already with p = 2.

Conclusions
In this work we have focused on two recent computational approaches, namely the Direct Parametrization of Invariant Manifolds (DPIM) and Deep Learning-based ROMs (DL-ROMs), for the efficient numerical simulation of physics-based virtual twins of nonlinear vibrating multi-physics microstructures, providing innovative contributions along two main directions.
On one side, we have brought the applications addressed by POD-DL-ROMs to an unprecedented level by showing that complex dynamical effects can be properly analysed in a very efficient way. In order to address strongly nonlinear phenomena like internal resonances and autoparametric effects involving the interaction between different modes, we have proposed a new arc-length abscissa which serves as ordering parameter for the collection of snapshots. We have also proved that the POD-DL-ROM construction can be easily extended to multi-physics problems exploiting its non-intrusiveness. This latter indeed allows us to generate snapshots with any FOM, including commercial ones, since FOM solutions are the only simulation data required to train neural network architectures.
On the other side, for problems with geometric nonlinearities, we have validated the POD-DL-ROM approach against the very recently developed exact DPIM approach, stressing the striking analogy that exists between the two techniques: namely, both share a nonlinear encoding phase, the generation of a reduced dynamics model, and a likewise nonlinear decoding phase that allows reconstructing very accurately the full displacement field. The POD-DL-ROM accurately reproduces the exact invariant manifolds with a minimal set of reduced variables. The problems addressed are characterized by low damping, a distinctive feature of most MEMS working in near-vacuum and these systems can be classified as nearly conservative. In this context, the POD-DL-ROM approach identifies the major features of the response, i.e., the FRFs together with the trajectories of the N master modes, with only N coordinates. On the contrary, slave modes might be active in the analysis and their correct reproduction generally requires to use up to 2N reduced variables, as predicted by the DPIM.
The insight gained in this work, associated with the real-time performances of deep learning-based reduced order models, paves the way to the application of the approach in the generation of digital twins of the modelled devices, to be used in the analysis and optimization of the complex systems in which the MEMS are employed, with countless applications in the pervasive Internet of Things.  Figure c): comparison between the p = 1 case and the FOM, clearly highlighting that main dynamic features are not represented adequately with one reduced variable. Figure a): comparisons for p = 2, 3. The curves below a threshold V AC correspond to a simple harmonic resonance of the softening drive mode; a plateau starts developing when the sense mode gets autoparametrically activated. Data for the sense mode are collected in Figures d) and f) illustrating the FRFs of the radial displacement of the node indicated by a red circle in Figures e), which is representative of the sense mode. These FRFs confirm the strong and explosive mode interaction according to which the sense mode reaches abruptly amplitudes of the same order as the drive one. The FRFs here reported have been obtained with a testing dataset made of 1 398 000 instances. The inquiry of the DL-ROM is performed in less than 3 s using a Tesla V100 32GB GPU and an implementation in the Tensorlow DL framework.