Sensitivity analysis of nonlinear frequency response of defected structures

The computation of the steady-stateresponse of large finite element discretized systems subject to periodic excitations is unfeasible because of excessive run time and memory requirements. One could in principle resort to reduced order models stemming from the high fidelity counterparts, which typically require a solution time orders of magnitude smaller. However, when many simulations are required, as in the case of parametric studies, the overall effort could be still significant and the analysis process could be severely hindered. In this work, we propose a sensitivity approach to assess the influence of model parameters on the nonlinear dynamic response. As opposed to the costly evaluation of reduced order solutions over a range of excitation frequencies and model parameters, the sensitivities of a nominal response allow one to approximate the dynamic response by a simple evaluation of an expansion in the directions spanning the parameter space. Special care must be taken on the closure equation that needs to be appended to the system of equations stemming from the harmonic balance method. We discuss the limitations of the current constant frequency approach and propose an improvement. We demonstrate the merits of the proposed approach on a micro-electro-mechanical system affected by parameterized manufacturing defects. Leveraging from a previous contribution, the nonlinear response and the sensitivities are obtained from a reduced order model which is analytical in the defect parameters. Our procedure is able to deliver accurate probability density functions of quantities of interest (e.g. nonlinear resonance peaks, triple solution bandwidth, etc) against statistical distributions of manufacturing defects at negligible computational cost.


Introduction
Numerical methods for steady-state solutions of harmonically forced nonlinear systems are well established in the literature. A popular method is the Harmonic Balance (HB) [1], a frequency-based approach that relies on the expansion of the system steady-state response in truncated Fourier series. This results in a set of nonlinear algebraic equations for the unknown Fourier coefficients, which can be solved to obtain the steady-state response. Continuation algorithms are employed together with HB to compute the Frequency Response Curve (FRC), that is the set of steady-state solutions for varying excitation frequency. For many years, HB has been successfully applied to study deterministic systems.
In the last decade the field of nonlinear dynamics has made some efforts towards the investigation of parameter uncertainties on the FRCs. In order to model random processes and/or to take into account the lack of knowledge on the real system (epistemic uncertainty), mathematical models may feature uncertain input parameters. Common sources of uncertainties in mechanical systems are: variable geometry introduced by manufacturing tolerances, non uniform material properties, nonuniform contact surface parameters, non-ideal boundary conditions, or in case of epistemic uncertainties, inaccurate measured system properties. Uncertainties in the parameters result in non-deterministic frequency responses that can be investigated with Uncertainty Quantification (UQ) techniques.
In this article we focus the attention on uncertainties introduced by small shape imperfections in geometrically nonlinear structures. Defects in the geometry, stemming for example from manufacturing process, may result in non-negligible modifications of the nonlinear frequency response, as in the case of MEMS devices [2,3].
The Monte Carlo Method (MC) [4] is one of the first tools developed in the literature for UQ and is still nowadays considered the benchmark for solution accuracy. This classical approach is based on repeated simulations of the nonlinear system for a set of deterministic parameters, randomly sampled from the stochastic parameter distribution. When the quantity of interest is the stochastic FRC, this translates in a repeated application of HB. For large systems, however, solving the nonlinear algebraic HB equations could be computationally unfeasible even for one simulation. It is for instance the case of MEMS, for which the equations of motions (EoM) are derived from a FE discretization that requires detailed meshing for a thorough description of complex geometry.
A possible option to overcome this issue is to apply HB to Reduced Order Models (ROM) of the original model, namely the Full Order Model (FOM). Projection-based ROMs can be constructed by selecting a suitable and reasonably small collection of vectors (reduced order basis) to approximate the full solution [5]. In this way, the number of degrees of freedom (dof) of the original system is significantly reduced, making HB computationally feasible.
Unfortunately, rate of converge of MC algorithm is slow, thus the number of responses required to obtain relevant statistical quantities is usually large. For this reason, even if computationally feasible, performing a MC analysis based on ROM to propagate uncertainties for nonlinear structures with stochastic imperfections, results in high computational cost. In fact, this approach sees the construction of a new ROM for each of the randomly sampled geometries. The cost associated to a single run for this operation is usually significant and scales with the size of the FOM, strongly impacting on overall time for MC. A possible solution to this problem is presented in [6,7] and relies on the construction of a Parametric Reduced Order Model in the uncertain structural defects (DpROM). This operation is done only once upstream of MC and represents a fixed cost, independent of the number of parameter samples. With this latter approach, the computational cost associated to a single simulation run is considerably cut down, so that the overall computational burden for MC is also greatly reduced. Still, as it was described in [6,7], DpROM requires additional generalized coordinates to properly describe the contribution of imperfections. For this reason, if on the one hand the explicit parametrization eliminates the cost associated to the ROM construction for every new parameter realization, on the other hand the increased size of the ROM impacts on the simulation times.
All in all, it seems more sensible and efficient to perform MC studies relying on methods which do not involve repeated simulations. Polynomial Chaos Expansion (PCE), Taylor series expansions and interval arithmetic are some methods falling into this category. PCE-based methods consist in expanding the HB Fourier coefficients of the response in a truncated series of orthogonal polynomials of the uncertain parameters. Coefficients of chaos polynomials can be determined following two different approaches: intrusive approaches based on a recast of HB equations to account for stochastic parameters or non-intrusive approaches that rely on nullification of the error at selected collocation points. Once coefficients of expansion are known, any large set of responses corresponding to different uncertain parameters could be obtained at almost null computational cost by multiple polynomial evaluations. In [8,9], an intrusive PCE in combination with Multi Harmonic Balance was successfully applied to quantify the effects of uncertainties on the FRF of a linear rotor system. By the same authors, the PCE-based method was modified in order to propagate uncertainties in systems with geometric nonlinearities [10,11]. In that context, the excitation frequency was considered non-deterministic but parameter dependent and thus was expanded in the polynomial chaos basis. This modification of the method was necessary to assure precise approximations, even in correspondence of returning points of the FRC, where multiple solutions for same excitation frequency may exist. In the computation of chaos polynomial coefficients, the author compared points on different FRCs satisfying a prescribed phase condition. In [12], a non-intrusive PCE-based method is used with the Asymptotic Numerical Method -Harmonic Balance (ANM-HB) [13]. In this case, spurious oscillations in the set of responses provided by PCE were successfully avoided by comparing points on different FRCs with the same "arclength ratio", a measure defined by the author. In [14][15][16] instead, stochastic amplitudes and excitation frequency were referred to the phase of first harmonic, in the non-intrusive computation of PCE coefficients. This approach is restricted to the case in which the phase is a monotonic function of the frequency for all the uncertain parameter values. In [17], the authors applied PCE to investigate the effect of uncertainties on the steady-state response of a Jeffcott rotor subjected to rub-impact and parameter uncertainties.
In [18], a surrogate model of the system is created exploiting an expansion of the response in Chebyshev Polynomials to which follows an interval problem to find the bounds of the FRC corresponding to prescribed parameter intervals. Coefficients of polynomial approximations are found by comparing points on different FRCs at the same polar angle in the frequencyamplitude plane, in the "Polar Angle Interpolation" approach as defined by the author. In [19], Chebyshev Polynomials have been exploited to model uncertainties in the steady state response of a system with backlash nonlinearity, in the framework of Interval Harmonic Balance Method (IHBM).
Taylor series-based UQ methods are far less popular than PCE methods in the field on stochastic nonlinear FRCs, probably for limited applicability to small parameter uncertainties. In [20], the author expanded the HB Fourier coefficients of nonlinear FRCs in a second-order Taylor series of the uncertain parameters. Coefficients of the Taylor series, namely Parameter Sensitivities, allowed easy prediction of frequency response bounds corresponding to small parameter intervals. Instead in [21], parameter sensitivities were used to asses statistical distributions of nonlinear FRCs of a bladed disk, for stochastic distribution of friction parameters, showing good agreement with MC results.
Both in [21] and [20], sensitivity-based approximations proved to accurately describe FRCs for small parameter variations from the nominal. However, even if in a nonlinear setting, in the test cases considered by the author, FRCs did not exhibit multiple solutions for the same excitation frequency. Moreover, sensitivity were performed at constant frequency. As it will be demonstrated in the present work, this approach to sensitivities leads to poor results in case of multiple solutions for same frequency. To deal with this problem, a new approach to sensitivity-based approximations of FRCs is here presented, in which frequency is allowed to vary. Once sensitivities are available, any FRC corresponding to any small perturbation in the uncertain parameters can be obtained by updating the nominal solution with Taylor series approximations. This "update" operation is practically inexpensive and can be exploited in a sensitivity-based MC approach to compute the FRCs.
In this paper we present a Sensitivity-Based MC (SB-MC) method for UQ of structures with geometrical imperfections, leveraging the availability of the parametric ROM (DpROM) presented in [7] and the small variance of flaws in real engineering applications.
The work is organized as follows: Section 2 recovers the fundamentals of HB and introduces the notation used throughout the paper; Section 3 reviews the sensitivity approach for nonlinear FRCs at constant frequency and presents our new approach, showcasing the main results and differences on a simple 1-dof problem. Section 4 sketches the main features of the DpROM which is then used in Section 5 to study a MEMS resonator. Finally, in Section 6 conclusions are drawn.

Harmonic balance for mechanical systems
We here briefly describe the HB method and introduce notation; the interested reader is referred to [1] for a complete dissertation. Let us consider the EoM of a ndof nonlinear mechanical system, stemming from the FE discretization of the continuous problem: where q ∈ R n is the vector of generalized coordinates, M, D, K ∈ R n×n are respectively the mass, damp-ing and stiffness matrices, f is the vector of nonlinear forces and f ω is an external periodic forcing with angular frequency ω, defined as f ω (t) = f ω,c cos(ωt) + f ω,s sin(ωt). (2) At steady state the response is assumed to be periodic, with the same period of excitation T = 2π ω and is expanded in a truncated Fourier series with H harmonic components as [1, cos(ωt), sin(ωt), . . . , cos(H ωt), sin(H ωt)] , where ψ ∈ R 1×(2H +1) is the basis of the Fourier expansion, I n is the identity matrix of dimension n, Q ∈ R n(2H +1)×1 is the unknown vector of harmonic amplitudes (with Q 0 , Q c,k , Q s,k , ∈ R n×1 , for k = 1, . . . , H ) and ⊗ denotes the Kronecker product. By taking time derivatives of the ansatz, expressions for velocities and accelerations writė where the derivative operator ∇ is defined as From Eq. (1), the residual in time domain is defined as Time dependence is eliminated by applying a weighted residual approach in which the weight functions are the components of the Fourier basis ψ. In this way we define the HB residual as . . .
where R( Q, ω) ∈ R n(2H +1) and Notice that the so defined HB residual is nothing but the collection of the first H+1 Fourier coefficients (using the sine-cosine representation) of the time domain residual r(t).
By plugging Eqs. (3a), (4) and (5) into Eq. (6) we obtain in which the expression of the linear residual R lin was derived by applying the mixed-product property 1 of the Kronecker product. Considering for example the residual of mass related terms we have that where we also have exploited orthogonality properties of Fourier basis functions. Moreover, the residual vector of forcing F ω in Eq. (7) can be developed as Nullifying R in Eq. (7) leads to the HB residual equation, that is a set of nonlinear algebraic equations in the unknown vector Q:

Numerical solution
Usually the analyst's interest is not in a single frequency response, but in the FRC, that is a branch of solutions of the HB equations in a frequency span. For this reason, ω is taken not as a constant parameter, but as an additional unknown, so that the FRC can be obtained using a Newton-Raphson method along with a pseudo arc-length continuation algorithm. According to this scheme, discrete points X i = [ Q i T , ω i ] T on the FRC are obtained, at the i th step of the continuation algorithm, by solving: where p (X i , X i−1 , ds) ∈ R is the additional closure equation depending on the adopted continuation algorithm, X i−1 is the solution of the previous continuation step and ds is the user provided value of the pseudo arc length.
At each continuation step the iterative Newton-Raphson method is employed to compute the solutions of the nonlinear algebraic system of equations: where the subscript k denotes the iteration number and J is the Jacobian matrix. While the linear part of the residual and of the Jacobian can be trivially evaluated, usually their nonlinear counterpart must be computed numerically, although in some special cases analytic solutions are available (e.g. for polynomials). Nonlinear residual is numerically computed by means of the Alternating Frequency-Time (AFT) algorithm, while we refer the reader to Appendix C for details on Jacobian evaluation. According to the AFT scheme, displacement vector q(t) and velocity vectorq(t) are evaluated at equidistant time snapshots in the period from the value of the Ansatz Coefficients Q. Upon substitution of these quantities into the nonlinear force relation f (q,q), time snapshots of the nonlinear force vector are retrieved. Eventually the integral in Eq. (7) is numerically approximated by applying the Discrete Fourier Transform (DFT) to the snapshots of nonlinear forces, efficiently implemented with the Fast Fourier Transform (FFT) algorithm. It is important to remark that the number of time samples N in the AFT algorithm impacts on accuracy of the numerically computed residual. If N is not sufficiently large, aliasing phenomena may occur leading to inexact residual evaluation. For the derivatives, because of nonlinearities, the number of required time samples is usually larger. For more details on the choice of N the reader is again referred to [1].

Sensitivity analysis
The nonlinear mechanical system in Eq. (1) can depend on a set of parameters, here collected in vector p ∈ R m . In this case the EoM read: in which the mass, damping and stiffness matrices along with the nonlinear force vector and forcing are assumed to be dependent on p.
By applying the HB method, we obtain a system of nonlinear algebraic equations parametrized by p: In many structural applications p can be used to model a source of uncertainty acting on the system. The response can be computed for a nominal value of parameter p * and variations in the response can be investigated by considering its sensitivity to p. In this context we define a nominal solution point of the parametric HB equations as a point { Q * , ω * , p * } satisfying: In the following, • * denotes a quantity referred to the nominal conditions, although often the superscript will be omitted to ease notation.

Sensitivity at constant frequency
An approach commonly followed in the literature [20][21][22][23] is to compute the sensitivity of the response at constant excitation frequency ω. For each nominal solution point, first-order sensitivities are obtained by total differentiation of the parametric HB equations in Eq. (14) with respect to the parameter vector p, at fixed nominal angular frequency ω * : where all partial derivatives are implicitly assumed to be evaluated at the considered nominal solution point { Q * , ω * , p * } and under the hypothesis of non-singular Jacobian. Through further differentiation we get where we denoted with "·" the contraction operation between last dimension of the first (multidimensional) matrix and the first dimension of the second matrix, while similarly "· i j " denotes the contraction between the i-th dimension of the first matrix and the j-th dimension of the second one 2 . Moreover, the order of elements in derivatives corresponds to the order of derivation (see Appendix C). Again, all quantities are assumed to be evaluated at the considered nominal solution point, while the first-order sensitivity matrix is obtained from Eq. (16). From Eq. (17), second-order sensitivities can be retrieved by solving m 2 systems of linear equations in n(2H + 1) unknowns: where I, J ∈ {1, . . . , m} and the sign ":" indicates all elements in the corresponding matrix dimension. The computational burden can be reduced exploiting Schwartz theorem so that the number of systems of equations to be solved for is reduced to m 2 2 + m 2 . Notice also that the solution of Eqs. (16) and (18) involves the inversion of the same Jacobian, thus enabling a faster solution process. The expressions for the partial derivatives of the residual R with respect to Q and p appearing in the definition of the sensitivities are given in full in Appendix C.
Finally, it is possible to obtain Taylor series approximations of the FRCs corresponding to small perturbations of the parameter vector from the nominal as where Q * is the nominal coefficient vector and δ p ∈ R m is the perturbation of the parameter vector from its nominal value p * .

A practical case study
In order to asses quality of sensitivity-based approximations, let us consider the practical case study of a Duffing oscillator subjected to Coulomb friction. The equation of motion for the harmonically excited system is where is the friction force regularized with the hyperbolic tangent smoothening function, f lim is the limit friction force and is the velocity tolerance [24]. A sketch of this system and of the regularized friction law is provided in figure 1. Sensitivities of the system's response are computed at constant frequency with respect to the parameter vector p = [k, k 3 , f lim , F] T , taking as nominal parameters k = 1, k 3 = 1, f lim = 0.05 and F = 0.08, with m = 1, c = 0.06 and = 0.02.
In figure 2, we show the sensitivity-based approximations for a −5% reduction of f lim (Fig. 2a), for a +5% increase of the forcing F (Fig. 2b), for a −4% reduction in both the linear and cubic stiffness coefficients k and k 3 (Fig. 2c) and for all these variations at the same time (Fig. 2d), that is for δ p % = [−4, −4, −5, +5] T . This latter parameter perturbation results in a significant increase in the resonance peak amplitude, since external forcing increases while the stiffness and the limit friction force decrease, hence it is a good test case for evaluating sensitivity-based approximations.
Both first-and second-order sensitivity-based approximations well capture the variation in the response resulting from the perturbations of the parameters. As expected, second-order sensitivity-based approximations are more accurate, especially close to the resonance peak. However, if the nominal forcing is increased from F = 0.08 to F = 0.12, considering once again the percentage parameter variation δ p % = [−4, −4, −5, +5] T , some points on the perturbed FRC are not anymore approximated correctly by the sensitivities, as shown in Fig. 3. This is due to the fact that sensitivities at constant frequency do not allow the points on the nominal FRC to move to neighboring points on the perturbed FRC, as they are at different frequency.

Sensitivity with the normal-direction approach
To overcome the limits of approximations based on sensitivities computed at constant excitation frequency, we introduce a new approach in which the frequency ω is allowed to vary with the parameter vector p. First, let us collect the harmonic amplitude vector Q and excitation frequency ω in a unique column vector as We now seek a unique mapping from each nominal solution point to a point on the FRC corresponding to a slightly perturbed parameter vector p = p * + δ p. More precisely, we implicitly define X as a function of p only as where U and V are open sets in the neighbourhood of p * and X * , respectively, and g(X( p), p) ∈ R is an additional closure equation. Notice that all the points belonging to a branch of the FRC in V , corresponding to any parameter value p ∈ U , satisfy the first condition R = 0 in Eq. (24), therefore the additional closure equation is essential to identify a unique point on the FRC, as qualitatively depicted in Fig. 4a.
In other words, this scalar equation ensures that the number of equations equals the number of unknowns, that is dim(R) = dim(X) = n(2H + 1) + 1. If R (X * , p * ) = 0 and if det ∂R ∂X * = 0, the Implicit Function Theorem ensures the existence of X( p), U and V as defined in Eq. (24).
Differentiating the implicit function defined in Eq. (24) it is then possible to compute the total derivatives of X with respect to p, and thus first-and second-order sensitivities as and for I, J ∈ {1, . . . , m}. Again, refer to Appendix C for the expressions of all partial derivatives of R with respect to X and p.
Eventually, sensitivity-based approximations are obtained as Notice that in the constant frequency approach, discussed in Section 3.1, the implicit function definition with g = ω( p) − ω * is tacitly assumed before taking the total derivatives with respect to p. A qualitative sketch of this parametrization is provided in Fig. 4b. The choice of the closure equation g has a strong impact on the accuracy of sensitivity-based approximations. As already mentioned, while many points satisfy the HB residual equation (blue dots in Fig. 4a), not all of them are at the same distance from the starting nominal point (black dot). As the sensitivity-based solution uses a Taylor expansion, the more the mapping from X * to X( p * + δ p) is close to linear in the expansion direction, the greater the accuracy we can expect. To this purpose, let us consider a nominal solution point {X * , p * } and a first-order expansion of the HB residual in (14): where δX = X − X * and δ p = p − p * . For a perturbation δ p of the parameter vector, we obtain an approximated FRC by nullifying the first variation of the HB residual δ R and solving for δX the resulting one time underdetermined system: The associated homogeneous system writes where τ * is one of the ∞ 1 solution vectors and is tangent to the nominal FRC at the nominal solution point. Thus, the linear system in Eq. (29) admits solutions of the form in which δX p is a particular solution of Eq. (29), c ∈ R an arbitrary constant and τ * a particular non-null solution of Eq. (30). Therefore, in a neighbourhood V of {X * , p * }, the linear approximation of the FRC corresponding to p = p * + δ p is given by the set of points {X : X = X * + δX p + cτ * , c ∈ R}, while the linear approximation of the nominal FRC, corresponding to p = p * , is given by the set of points {X : X = X * + cτ * , c ∈ R}. Thus, in the space of harmonic coefficients Q and frequency ω, the linear approximations of the two FRCs are two parallel lines, as sketched in figure 5a.
Since the Taylor series approximations in Eq.
In order to find a solution for Eq. (30), let us partition the tangent vector as τ * = [τ Q * T , τ * ω ] T and solve the linear system of equations with arbitrary non-null τ *   ω (a, b, c) or, neglecting the frequency dimension, in the coefficients space only (d, e, f). In (a, d), for the two mappings, modulus of the response is computed for the nominal solution (k * , ), and, for δk% = +10%, with second-order sensitivities ( ) and by re-running HB ( ). In (b-e), the two mappings (X (k), ) in the modulus-frequency plane. In (c,f), the two mappings in the Nyquist plot

Limitations of the closure equation
Let us now consider the nonlinear oscillator presented in Section 3.1.1, this time with m = 1, k = k 3 = 10, c = 0.05, F = 0.1 and f lim = 0. The forced steady-state response for this system can be accurately captured with a single harmonic resulting in a FRC in a three dimensional space (Q s,1 , Q c,1 , ω), allowing for an effective graphical representation of the mapping introduced by the constraint equation g.
In Fig. 6a we plot the sensitivity-based approximation for a +10% increase of the linear stiffness k. As it can be observed, in some frequency ranges the accuracy of the sensitivity-based approximation is low.
In Fig. 6b and 6c, we plot the mapping introduced by closure equation (32) for equidistant interval increase in stiffness k, in the modulus-frequency and coefficients space (Nyquist plot), respectively. The mapping (red lines) has been obtained by exactly solving Eq. (24) for different values of the parameter p = k, starting from different nominal solution points. It can be observed that increasing k the FRC shifts to higher frequencies and the resonance peak gradually decreases, as expected. However, it can also be observed that since we are mapping each nominal solution point to the closest point on the perturbed FRC considering also the distance in the frequency space, we end up associating different "operating points". This is clearly visible around frequency ω ≈ 3.4, where two red lines cross. The reason why this happens is that the proposed approach minimizes the distance between solutions in the Q − ω space. Ultimately, this attempt to minimize also δω makes the mapping highly nonlinear and harder to capture with a Taylor expansion, effectively reducing the parameter validity range of the approximation.

A modified closure equation
Following these considerations, we modify the closure equation in such a way that the distance in the frequency space is no longer weighted in the mapping. The modified closure equation we choose is where τ Q * is defined in Eq. (33). In other words, we map the nominal solution point to the point on the perturbed FRC such that the variation in the HB coefficient vector δ Q = Q( p) − Q * is orthogonal to the tangent restricted to the HB coefficients space. A qualitative geometric sketch of this parametrization is provided in Fig. 5c. As proved in Appendix A, with this closure equation we have that under the (generally met) two conditions det( ∂ R ∂ Q ) = 0 and ∂ R ∂ω = 0, the Implicit Function Theorem hypothesis are satisfied. Moreover, in Appendix B, we also prove that this modified approach minimizes ||δ Q|| in a linear approximation of FRCs.
Going back to the example of the nonlinear oscillator already presented in the previous section, with the same increase in stiffness of 10%, we now see that second-order approximations of the mapping based on closure equation (34) (Fig. 6d) are more accurate than those based on closure equation (32) (Fig. 6a). Moreover, notice that in this second approach the mapping is orthogonal to the FRC in the space of HB coefficients (space of A 1 = Q c,1 and B 1 = Q s,1 , in this example) as shown in Fig. 6f. This translates into a map where the modulus is approximately constant, as depicted in Fig. 6e, since now the frequency ω is free to move and the minimization is restricted to the norm of the HB coefficients variations.
A notable case is the one of pure translation in frequency of the FRC. Approximations based on closure equation (32), indeed, may lead to inaccurate results when the FRC is mainly subjected to a translation along the frequency axis for a given variation of the parameter vector. This is the case of the nonlinear oscillator introduced in section 3.1.1 with nominal parameters m = 1, k = k 3 = 10 3 , c = 0.2, F = 3 and f lim = 0 for an increase in linear stiffness δk = 10%. For this example we provide in Fig. 7 the HB coefficients Q s,1 ,Q c,1 and the modulus ||Q s,1 + i Q c,1 || (i is the imaginary unit) for the nominal and the perturbed FRCs, computed by running HB with a single harmonic. With reference to the current section, we recall that starting from each nominal solution point it is possible to obtain a linear approximation of the FRC corresponding to a perturbation in the parameters, through the linearization of the HB residual provided in Eq. (28). This approximation is given by the set of points {X : X = X * +δX p +cτ * } and is plotted in Fig. 7 with dashed line, for different nominal expansion points X * identified with markers " * ". Notice that each approximation describes well small different sections of the curves. In a sensitivity-based expansion (Eq. (27)) limited to first-order terms, with the choice of closure equation g, we uniquely identify only one of these points, that will become part of the final approximation of the perturbed FRC. With refer-  Fig. 7 this is equivalent to select the points on the linear approximations (dashed lines) marked with "×", if closure equation g weighs variations in frequency (Eq. (32)), or marked with "+", if frequency dimension is neglected (Eq. (34)). As it can be observed, in the first case, since distance in frequency is weighted, the closure equation "forces" the points on the perturbed FRC to be at a similar frequency of the corresponding nominal solution points. In this way, the method displays the same issues presented by the constant frequency approach, resulting in poor approximations. Conversely, in the second case, since closure equation does not weigh variations in frequency, approximation points can freely move in the frequency dimension, allowing for greater accuracy.

A practical case study: monte carlo analysis
Let us consider the example of the nonlinear oscillator with friction presented in Section 3.1.1 and with nominal parameter vector p * = [1, 1, 0.2, 0.12] T . The perturbed FRC for percentage parameter variation δ p % = [−20, −20, +5, +5] T , obtained with the sensitivities in the normal-direction approach formulation, is given in Fig. 8. Notice that the so computed approximation is smooth and solves the problems related to the constant frequency approach. Moreover, greater accuracy is achieved with respect to the case presented in Fig. 3 that featured the same nominal parameters, even if percentage parameters perturbations are increased.
Sensitivity-based approximations can be efficiently used to propagate uncertainty in the system's parameters to the response, following a Monte Carlo (MC) approach. From the knowledge of the stochastic distribution of the parameters, the stochastic distribution of the response is retrieved by computing a set of responses for a sufficiently large number of randomly sampled parameter perturbations.
In contrast to conventional Monte Carlo, instead of running the numerical solution of the nonlinear HB equations, the solution is cheaply obtained with the sensitivities by applying the update Eq. (27), for each parameter variation. This comes at the expense of an additional fixed cost to compute sensitivities, which is easily amortized as the number of MC samples is increased. Finally, stochastic distributions of selected Quantities of Interest (QoI) can be obtained by postprocessing the set of responses.
As a demonstration of accuracy and efficiency of the sensitivity-based MC, we apply the method to the nonlinear oscillator with friction, assuming different stochastic distributions of the parameters. The quantities of interest we consider are the resonance frequency and resonance amplitude. In Fig. 9 we show the Probability Density Function (PDF) for each of these quanti- ties, computed for a Gaussian distribution of the linear and cubic stiffness k and k 3 (Fig. 9a), for a Gaussian distribution of friction limit force f lim (Fig. 9b), for a Gaussian distribution of external forcing F (Fig. 9c) and for the combination of all these distributions (Fig.  9d).
Each of the stochastic distributions is obtained both with conventional MC (i.e. solving the HB equations for the perturbed system) and with the sensitivity-based MC, randomly sampling the parameter distribution 1000 times (same seed). As it can be seen, the distributions computed with the two different approaches are in good agreement. For this 1-dof example, we obtained a speedup factor of 15.

Hardening, softening and stability considerations
Before concluding this section, we report the frequency responses for an oscillator with no friction ( f lim = 0) and with cubic stiffness k 3 = 0 (Fig. 10a) and k 3 = 0.015 (Fig. 10b). The remaining parameters m = 1, k = 1, c = 0.06, F = 0.25 are fixed. In the present case, changing the cubic spring value has the effect of changing the type of nonlinearity from hardening to softening, and vice versa. However, as it can be observed, the accuracy of the approximation is not affected by the type of nonlinearity, but rather by the magnitude of the variation, as expected from Taylor expansion theory. In some cases, the maximum allowed variation can be analytically derived [13], but this usually requires the knowledge and computation of derivatives one order higher than the desired order of the approximation. Lastly, in Fig. 10 we also show with dotted lines the unstable branches of the solutions, obtained computing the eigenvalues of the monodromy matrix [1,25], obtained with the shooting method implemented in NLvib. Although stability information is usually desired, this is not readily available neither for the HB-computed solutions nor for the sensitivity-based approximation. This means that for each point of the solution one has to set a proper analysis to determine stability and, for the monodromy matrix method we used, a forward time integration has to be performed for each solution point. This has been done for the nominal and reference curves. However, using the shooting method, the underlying hypothesis is that the solution point (used to set the initial conditions of the time integration) is actually on a periodic orbit, which is in general not true for the sensitivity-based solution (or even for HB-solutions with an insufficient number of harmonics). For this reason, such a stability analysis proves unreliable for the approximated solutions and is not shown here. As a final remark, we also point out that the computational burden associated to the stability analysis of the FRC would have probably defeated the benefits of the proposed method. Until now we presented the sensitivity analysis method for a generic case, showing some examples with lumped parameter models. In this section, we want to apply the method to FE structures with the aim to study the effect of manufacturing variations on the shape of the structures themselves. This task poses at least two issues. First, even a single solution with the Harmonic Balance method can be too expensive (or even unfeasible) for a FE model with a moderate number of dofs; secondly, shape parametrization techniques often involve a high number of parameters, meaning that a high number of sensitivities would need to be computed. As anticipated in Section 1, we propose to solve these problems applying the sensitivity approach to a Defect-parametric Reduced Order Model (DpROM). The latter, introduced in [6] and extended in [7], is here described in brief, leaving the details on its implementation to the aforementioned works.

DpROM description
The DpROM is a projection-based ROM which can parametrically describe a shape imperfection in a regime of geometric nonlinearity. Shape imperfections, or defects, are simply given by displacement fields which move the nodes of the nominal (i.e. without imperfections) structure into the defected configuration. These displacement fields may come from a previous study, measurements, or they can be artificially constructed. Each displacement field U j can be independently defined, and the final defected configuration will be given by the linear superposition of all the defects, that is: where x d and x 0 are the defected and nominal configuration coordinates, respectively, U = [U 1 , . . . , U m ] is the matrix collecting m defect shapes, and where ξ = [ξ 1 , . . . , ξ m ] T is a column vector collecting each defect shape's amplitude. In this context, ξ takes the role of the parameter vector which governs the final defected shape of the structure under study. The projection basis V ∈ R n×r , being r the number of reduced coordinates η ∈ R r , is built using a set Vibration Modes (VM), their associated (Static) Modal Derivatives (MD) and Defect Sensitivities (DS), collected in the matrices Φ, Θ and Ξ , respectively, so that the displacement vector u can be approximated as This choice, loosely speaking, is an extension of classic modal analysis, where the truncated set of eigenvectors of the system is augmented with MD to take into account nonlinearities and DS to take into account shape parametrization. Retaining a number n vm of VM, we have n vm (n vm + 1)/2 MD and n vm m DS, for a total r = n vm (n vm +2m +3)/2. For more details on this projection method, again, see [7] and references therein.
Operating a Galerkin projection, we can write the governing equations for the DpROM as whereM = V T MV, and the nonlinear elastic forces take on a polynomial shape in the displacements 3 , the 3 Under the hypothesis of linear elastic constitutive law and total Lagrangian formulation, elastic forces feature linear, quadratic stiffness tensors also depending in a polynomial way on ξ as In Eqs. (38), the left-subscripts denote with a number the order of the tensor, with a "n" a tensor computed for the nominal geometry and with a "d" ("dd") a tensor whose last (two) dimensions are of size m (which are otherwise equal to r ). For instance,K 4dd ∈ R r ×r ×m×m . The full expressions for these tensors can be found in [7]. Furthermore, each of the above tensors and the reduced mass matrix writẽ whereP is a generic reduced matrix/tensor,P is evaluated on the nominal volume andP ,i is a corrective term required to approximate the volume change introduced by the i-th defect 4 . For this reason, the mass matrix is also parameter-dependent. In the present work, Rayleigh damping is assumed asD = α MM + β KK and can be introduced in Eq. (37). We will also assume that the external forcing f is approximately constant and does not depend on ξ . Finally, sensitivities can be computed for the DpROM using the derivatives of the stiffness and mass tensors given in Appendix D. All in all, DpROM features a small set of reduced coordinates and an easy way to parametrize shape imperfections which results into polynomials, that can be efficiently used to compute sensitivities. Finally, notice once more that although greatly reducing the computational burden with respect to its corresponding FOM and its non-parametric ROM counterparts, DpROM may still take quite some time if several thousands of runs are required, as in the MC method. In the numerical example, we will show how beneficial the sensitivity analysis is in this sense. and cubic terms, given by the contraction of second, third and fourth-order stiffness tensors over the displacement vector (one, two and three times, respectively).

Numerical test: MEMS oscillator
To demonstrate the potentialities of the proposed method, we now present the case of a single-axis MEMS gyroscope [26]. The structure is fixed in correspondence of the anchor points and is harmonically driven by comb-finger electrodes along the y-direction (see Fig. 11a). Due to this motion, when an external angular rate is present, the sense mass will experience a Coriolis force which will make the mass move in a direction orthogonal to the drive-axis (i.e. x or zaxes); this movement is then sensed by specifically designed electrodes to infer the angular rate [2]. Since the Coriolis-induced motion amplitude is typically 1-2 orders of magnitude smaller than the one along the drive-axis, manufacturing imperfections may have a great impact over the performances of these transducers and must be taken into account in the design phase.

Model and chosen defects
We constructed a FE model of the device using threedimensional continuum elements (Lagrangian hexahedra with quadratic shape functions), counting 451,089 dofs. We used the material properties of Silicon (Young's modulus E = 148 G Pa, Poisson's ration ν = 0.23, density ρ = 0.93 × 2330 kg/m 3 , where 0.93 is the filling fraction to take into account the holes required by the production process), a Rayleigh damping with α M = 71.4 and β K = 0 (corresponding to a quality factor of 2200 on the first vibration mode), and an equivalent concentrated forcing F = 0.5 μN located in the middle of the drive combs. The DpROM 5 is constructed using 2 VM, 2 shape defects, and their associated MD (3) and DS (4), for a total of 9 dofs.
The first shape defect is the so-called overetch, which corresponds to the erosion of the Silicon wafer happening during the MEMS fabrication process [27]. The latter may indeed result in excessive/deficient erosion of the material (over/underetching) with respect to the nominal overetch value, meaning that the final geometry will appear slightly thicker/thinner than the nominal one. In our FE model, this defect is obtained from a static linear analysis, performed in Comsol Multiphysics 6.0, where an inward/outward normal unitary displacement is imposed to all the vertical walls, while the top and bottom faces are fixed along the z-direction by rollers. This way, the ξ 1 parameter will correspond to the effective overetch (for ξ 1 > 0) or underetch (for ξ 1 < 0) variation. The results of the static analysis are shown in Fig. 12(a,b). The second defect shape is the so-called wall-angle, consisting in the fact that the vertical walls will have an angle with respect to the z-axis. The field we consider for our model, shown in Fig. 12c, is given by being u, v, w the displacement fields along x, y, zaxes, respectively, z 0 the elevation of the bottom surface, θ a typical (3σ -)value for the wall angle and ξ 2 ∈ [−1, 1] the amplitude parameter. Notice that with this choice we consider simultaneously an angle with respect to the xz and yz-planes.

Monte Carlo analysis
For each shape defect amplitude, we consider a Gaussian distribution with null mean value and 3σ -value equal toξ 1 = 1 6 andξ 2 = 1. We first preliminary check that the approximated solutions provided by the DpROM and the sensitivity approach are sufficiently accurate with respect to a solution obtained with the non-parametric corresponding ROM 7 . To this purpose, in Fig. 13(a, d) the FRCs measured at forcing location for ξ = [±ξ 1 ,ξ 2 ] are shown, corresponding to the 3σ -corners of the parameter domain (notice that the sign of ξ 2 does not affect the solution). The FRCs show a significant hardening behavior, induced by the geometric nonlinearities of the clamp-clamp suspension beams connected to the drive frame. Under the assumption that the approximation improves for lower defect amplitudes, we deem these results accurate enough, even though a small deviation from the reference (ROM) is present. Notice that this error is to impute to the DpROM approximation, whereas the sensitivity approach reproduces the DpROM solutions accurately. Upon these considerations and for time constraints, the ROM solutions will not be computed anymore, and the comparison will be only between the DpROM solutions and their respective sensitivitybased approximations. All the computations were carried out in YetAn-otherFEcode v1.3.0 [28], with the addition of dedicated routines for the DpROM and sensitivity computation. For computational efficiency, some parts of the DpROM are built using Julia subroutines for tensor operations [29], as it will be indicated later in the computational times (see Table 1). The FRC are computed in NLvib [1], which is also included in YetAnotherFE-code v1.3.0, using MATLAB 2021b on a local machine equipped with an Intel(R) Xeon(R) Silver 4214 CPU @2.20 GHz and 256 GB RAM @2666 MHz. For the HB method, we used H = 5 harmonics, a number N t = 3H + 1 of time samples for the AFT of the nonlinear terms [30] and twice this number of samples for the AFT of the sensitivities.
For the Monte Carlo analysis, we generated 5000 random ξ parameters vectors extracted from the aforementioned Gaussian distributions, then we ran as many simulations with the DpROM and computed as many solutions with the sensitivity approach. However, since Fig. 13 (a, d) Notice that this is a problem only for the DpROM, where a simulation must be performed, but not for the sensitivity approach, where the solution is simply an update of the nominal solution. However, for the sake of comparison, we discarded the failed cases for both the approaches. Figure 13b, c, e, f show the Probability Density Functions for some QoI, namely the frequency at maximum amplitude (b), the span of frequencies where three solutions exist (c), the maximum displacement along y (e) and along z (f), measured at the forcing location. The latter is usually referred to as quadrature error (in that it is a spurious motion which superimposes to the Coriolis reading), and it is nicely captured by the sensitivity method in spite of its small magnitude. As it can be observed, DpROM and sensitivity-based results are in almost perfect agreement.
Finally, in Table 1 computational times are reported, where we distinguish between fixed costs, computed once and for all independently from the total number of analyses, and variable costs that must be sustained for each new analysis. Considering only the variable costs, we find that the sensitivity approach is almost 300 times faster than the DpROM, and more than 3000 times faster than the repeated construction and evaluation of a non-parametric ROM, with a modest increment of the fixed costs. The latter, however, can be largely amortized when a sufficient number of samples is used for the MC analysis.

Conclusions
In this contribution we proposed a sensitivity-based method to alleviate the computation of the frequency responses of nonlinear systems, when the latter have to be repeated for a large number of parameter realizations and for small parameter deviations from the nominal values. The typical case scenario in this sense is the Monte Carlo analysis. After recalling the governing equations for the Harmonic Balance method, we discussed how to compute sensitivities and how appending an appropriate closure equation to the HB residual equations is essential to achieve accurate results. We then introduced a closure equation ensuring that the sought approximated solution lies along a direction which is normal in the HB coefficient space. With the aid of the example of a nonlinear oscillator, accuracy and limitations of the discussed approaches were thoroughly illustrated. Finally, the method is applied to a large FE model of a MEMS gyroscope. As the direct evaluation of the response of such a system through HB is not feasible, we resorted to a parametric Reduced Order Model which we studied in a previous work (DpROM), and where the parameters represent the amplitudes of some imposed shape defect. Using the HB and DpROM, a Monte Carlo analysis was run for 4612 sampled parameter vectors, computing the Frequency Response Curves for both the DpROM (reference solution) and the sensitivity-based approximation. The probability density functions of some quantities of interest for the two batches of results showed excellent agreement. The computational times, however, were greatly reduced, going from the 2.6 days of the re-simulated DpROM to the 45.2 minutes of the sensitivity-based approach, with an effective speedup factor of approximately 300 (if neglecting fixed costs). All in all, the proposed method represents an accurate and cheap solution to evaluate the frequency responses of systems subject to parameter uncertainties and featuring a marked nonlinear behavior, finding potential applications in the MEMS industry as well as in other fields, as for the study of friction related phenomena [31].
Acknowledgements The authors want to express their gratitude to the AMD R&D group at STMicroelectronics, Cornaredo (IT), and to professor Francesco Braghin for their support.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement. Not applicable.
Availability of data and material Not applicable.

Compliance with ethical standards
Conflicts of interest The authors declare that they have no conflict of interest.
Code availability The Matlab FE code and the Julia routines used in this work can be made available upon request to the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Cre-ative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A Non-singularity of Jacobian matrix
In this subsection we demonstrate that, with the implicit function definition in Eq. (24), conditions for the Implicit function theorem are satisfied if: (i) det( ∂ R ∂ Q ) = 0, for closure Eq. (32); (ii) det( ∂ R ∂ Q ) = 0 and ∂ R ∂ω = 0, for closure Eq. (34). The Jacobian of the extended residual is with α = 1 for closure equation (32), or α = 0 for closure equation (34). Since ∂ R ∂ Q is invertible for hypothesis, the determinant of the extended Jacobian can be developed using Eq. (33) as If α = 1, then condition (43) is always satisfied, since the norm of a tangent vector is positive definite. However, for α = 0, condition (43) is satisfied if and only if τ Q * = 0. As it can be seen from the definition of τ Q in Eq. (33), since ∂ R ∂ Q has full rank for hypothesis, τ Q * is null if and only if ∂ R ∂ω = 0. Hence, in this second case, for ∂R ∂X to be invertible, the additional condition

D Derivatives for the DpROM
With reference to EoM of the DpROM (Eq. (37)), the vector of internal reduced forces writes Recalling the definitions for stiffness tensorsK given in Eq. (38) and volume integration formula in Eq. (39), internal forces and mass tensor derivatives, required to compute the sensitivities for the DpROM write in Einstein notation: where we highlighted dimensions over which summation is performed with indices in lowercase. All derivatives are evaluated at the nominal parameter value ξ = 0.