Cherif Snoun

Baptiste Bergeot

Sébastien Berger

Robust optimization of nonlinear energy sinks used for mitigation of friction-induced limit cycle oscillations

Keywords: Friction-induced vibration, Nonlinear Energy Sink, Uncertainty, Robust optimisation, Polynomial Chaos

This paper aims at proposing robust methods to optimize nonlinear energy sinks (NES) used for the mitigation of friction-induced vibrations due to mode coupling instabilities. The study is based on a mechanical system composed of two NES coupled to the well-known two-degrees-of-freedom Hultèn's model. In such an unstable system coupled with NES, it is usual to observe a discontinuity in the steady-state amplitude profile which separates the parameter space into two parts corresponding respectively to the mitigated and unmitigated regimes. The discontinuity is predicted by a methodology previously developed by the authors and based on Multi-Element generalized Polynomial Chaos. The method allows to determine the Propensity of the system to undergo a Harmless Steady-State Regime (PHSSR). The objective of the present work is therefore to maximize the value of the PHSSR to obtain a robust optimal design of the NES. To this end, several stochastic optimization problems are presented that take into account the dispersion of the uncertain parameters using two approaches; in the first one, the parameters of the NES are considered as deterministic, and in the second one they are also supposed uncertain but with a known probability law.

Introduction

Friction-induced vibrations rank among the major issues that friction systems may face. This phenomenon can result from the generation of Limit Cycle Oscillations (LCO) induced by dry friction, which are explained by the coupling of two modes of the system [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]; [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF]; [START_REF] Hervé | Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects[END_REF]). Given the parametric uncertainties present in friction systems, it is very difficult to design them so that they always have stable behaviors. It is therefore necessary to use strategies to attenuate vibrations when they appear. As a way to model this phenomenon, the well-known two-degrees-of-freedom Hultèn's model [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF][START_REF] Hultén | Brake squeal -a self-exciting mechanism with constant friction[END_REF]) has been widely used in the past.

Nonlinear Energy Sinks (NES) are nonlinear passive vibration absorbers which consist in a spring mass damper with a strong nonlinear stiffness (usually cubic as is the case in this work). The NES can adapt themselves to the Primary Structure (PS) without being tuned to a specific frequency. Their operation is based on the concept of Targeted Energy Transfer (TET) which has become an important passive control technique for reducing or eliminating unwanted vibrations [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF]; [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]. The design and optimization of NES constitute a very dynamic field of research. Recently, a lot of work has for example been produced concerning the study of new configurations and technologies allowing/enhancing TET [START_REF] Tian | Nonlinear flutter suppression and performance evaluation of periodically embedded nonlinear vibration absorbers in a supersonic FGM plate[END_REF]; [START_REF] Xiong | Effect of Electromechanical Coupling on Dynamic Characteristics of a Piezoelectric Nonlinear Energy Sink System[END_REF]), the analysis of the underlying dynamics of the coupled system [START_REF] Bergeot | Scaling law for the slow flow of an unstable mechanical system coupled to a nonlinear energy sink[END_REF]; [START_REF] Habib | Tracking modal interactions in nonlinear energy sink dynamics via high-dimensional invariant manifold[END_REF]; [START_REF] Bergeot | Dynamic behavior analysis of a mechanical system with two unstable modes coupled to a single nonlinear energy sink[END_REF]), the deterministic optimization of NES [START_REF] Geng | Nonlinear energy sink with limited vibration amplitude[END_REF]; [START_REF] Khazaee | A comparative study on optimization of multiple essentially nonlinear isolators attached to a pipe conveying fluid[END_REF]), the stochastic optimization of NES Pidaparthi andMissoum (2018, 2019), the coupling of the NES to an energy harvester [START_REF] Fang | Tuned bistable nonlinear energy sink for simultaneously improved vibration suppression and energy harvesting[END_REF]; [START_REF] Karama | Energy harvesting in a nonlinear energy sink absorber using delayed resonators[END_REF]; [START_REF] Karličić | Periodic response of a nonlinear axially moving beam with a nonlinear energy sink and piezoelectric attachment[END_REF]) or the use of a NES as an energy harvester [START_REF] Zhou | Essentially nonlinear piezoelectric shunt circuits applied to mistuned bladed disks[END_REF]; [START_REF] Li | Potential of a vibro-impact nonlinear energy sink for energy harvesting[END_REF]), etc. Among the many research subfields about NES, the present work takes place within the framework of robust optimization of NES when the latter have uncertain parameters and are used to mitigate Limit Cycle Oscillations (LCO).

The mitigation of LCO by NES has also been widely studied in the past. The first work reported in this framework concerns the LCO mitigation of the Van der Pol oscillator [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive non-linear energy sinks[END_REF]). In the field of friction-induced vibrations mitigation, [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] analyzed the behavior of a Hultèn's model coupled to two ungrounded NES; the same model in considered in this work.

Lately, special attention has been paid to the optimal design of NES in a deterministic context. [START_REF] Nguyen | Design criteria for optimally tuned nonlinear energy sinks-part 1: transient regime[END_REF] have studied the effectiveness of a NES to attenuate free oscillations of a single degree-of-freedom oscillator. The optimization procedure consisted in studying the influence of the NES damping on the efficiency of the energy transfer from the primary system to the NES. [START_REF] Wang | Multiobjective design and performance investigation of a high-rise building with track nonlinear energy sinks[END_REF] presented multi-objective designs of a track NES and a single-sided vibro-impact track NES used as January 13, 2022 effective control devices to mitigate the seismic response of high-rise buildings. In [START_REF] Qiu | Design criteria for optimally tuned vibro-impact nonlinear energy sink[END_REF] a design criteria for optimally tuned Vibro-Impact NES to control vibration under periodic and transient excitation has been proposed. [START_REF] Oliva | Optimal design of nonlinear energy sinks for {SDOF} structures subjected to white noise base excitations[END_REF] developed an approximate design approach based on the use of the Statistical Linearization Technique, and an accurate empirical formulation linking the NES optimal parameters to the characteristics of the main structure and the random excitation.

More recently, the optimization of NES under uncertainties has also been considered. [START_REF] Boroson | Optimization under uncertainty of nonlinear energy sinks[END_REF] has introduced the "discontinuity in the effectiveness" of the NES, i.e., the fact that the effectiveness can change from a high value to a low value following a small perturbation of the design parameters or of the initial conditions. In [START_REF] Boroson | Optimization under uncertainty of parallel nonlinear energy sinks[END_REF], the optimal design of several parallel NES configurations for maximum mean efficiency has been investigated. It has been shown again that the variation of some system parameters or of the initial conditions generates a discontinuity on the energy shape of the system. An optimization algorithm based on the Support Vector Machine (SVM) technique and a kriging meta-model have been used to predict this discontinuity according to the dispersion of some uncertain parameters. [START_REF] Khazaee | A comparative study on optimization of multiple essentially nonlinear isolators attached to a pipe conveying fluid[END_REF] studied the stochastic optimization of multiple NES (configured in parallel and in series) attached to a simply supported pipe conveying fluid. The results show that the optima obtained from stochastic optimization are much more robust than those obtained from deterministic optimization. Pidaparthi andMissoum (2018, 2019) performed optimization of NES under uncertainties in the context of attenuating self-sustained vibrations of an aircraft wing created by aeroelastic instability (which is also due to mode coupling). Depending on the nonlinear properties of the system, the latter may undergo LCO that are either supercritical or subcritical. In this case the bifurcation diagram of the system (which represents the amplitude of the LCO as a function of the airflow velocity over the wing) may exhibit two types of discontinuities, respectively due to the operation of the NES, and to the generation of subcritical LCO. Using again the SVM technique and a kriging meta-model, several stochastic optimization problems have been presented in order to maximize the average reduction of the LCO amplitude.

Various probabilistic approaches may be used to study the propagation of uncertainties in the deterministic model (DM) of a mechanical system. A review of the numerical methods dedicated to the stochastic prediction can be found in [START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF]. Efficient methods such as the generalized Polynomial Chaos (gPC) or the Multi-Element generalized Polynomial Chaos (ME-gPC) have been developed as a less costly alternative to the common Monte Carlo approach. These methods have been for example used to perform a stochastic analysis of the dynamic behavior of friction systems [START_REF] Nechak | A polynomial chaos approach to the robust analysis of the dynamic behaviour of friction systems[END_REF][START_REF] Nechak | Prediction of random self friction-induced vibrations in uncertain dry friction systems using a multi-element generalized polynomial chaos approach[END_REF][START_REF] Nechak | Non-intrusive generalized polynomial chaos for the robust stability analysis of uncertain nonlinear dynamic friction systems[END_REF][START_REF] Nechak | Robustness of stochastic expansions for the stability of uncertain nonlinear dynamical systems -application to brake squeal[END_REF]; [START_REF] Sarrouy | Stochastic analysis of the eigenvalue problem for mechanical systems using polynomial chaos expansion-application to a finite element rotor[END_REF][START_REF] Sarrouy | Stochastic study of a non-linear self-excited system with friction[END_REF]; [START_REF] Trinh | Stability analysis of a clutch system with multi-element generalized polynomial chaos[END_REF]; [START_REF] Snoun | Prediction of the dynamic behavior of an uncertain friction system coupled to nonlinear energy sinks using a multi-element generalized polynomial chaos approach[END_REF]).

In this paper, two original optimization methods under uncertainties, based on polynomial chaos theory, are developed for the robust design of two ungrounded NES used to mitigate LCO due to mode coupling instability in the Hultèn's model, which is the primary structure (PS). [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] classified the steady-state response regimes in two main categories related to the dispersion of some uncertain parameters: the mitigated regimes (harmless situations) and the unmitigated regimes (harmful situations). In the context of instability mitigation by means of NES, the LCO amplitude profile as a function of a given bifurcation parameter usually presents a discontinuity between these two regimes, which makes the NES potentially very sensitive to uncertainties. In [START_REF] Snoun | Prediction of the dynamic behavior of an uncertain friction system coupled to nonlinear energy sinks using a multi-element generalized polynomial chaos approach[END_REF], this discontinuity is predicted by means of a ME-gPC based method. This allows to determine the Propensity of the system to undergo a Harmless Steady-State Regime (PHSSR), which is defined as the probability for the system to undergo a mitigated regime. The aim of this paper is therefore to maximize the PHSSR value to obtain a robust optimal design of the NES with a low computational cost. For this purpose, several stochastic optimization problems are presented for a PS having uncertain parameters, in which the NES parameters are either considered as deterministic, or are also supposed uncertain but with a known probability law.

The article is organized as follows. In Section 2, the two degrees-of-freedom Hultèn's model coupled to two NES is presented. In Section 3, the possible steady-state regimes of the coupled system are studied, and the PHSSR is defined. Section 4 describes the polynomial chaos theory and recalls the stochastic optimization algorithm used to compute the PHSSR. The optimization formulation under uncertainties is described in Section 5, and the results are presented in Section 6. Finally, conclusions are given in Section 8.

The friction system under study

The studied mechanical system is similar to that previously used by the authors in [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF]; [START_REF] Snoun | Prediction of the dynamic behavior of an uncertain friction system coupled to nonlinear energy sinks using a multi-element generalized polynomial chaos approach[END_REF]. It is recalled hereafter for the sake of clarity.

The primary system

The primary system (PS) consists in the two degree-offreedom (DOF) Hultèn's phenomenological model [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF][START_REF] Hultén | Brake squeal -a self-exciting mechanism with constant friction[END_REF]). This simple model can reproduce the modecoupling phenomenon and is composed of a mass m held against a moving strip. The contact between the mass and the strip is modeled by two plates supported by two different springs, with linear components k 1 and k 2 and cubic components k N L 1 and k N L 2 , and two dampers with damping coefficients c 1 and c 2 (see Fig. 1). We assume that the friction coefficient is constant and that the strip moves at a constant velocity. The relative velocity between the strip velocity and dx 1 /dt or dx 2 /dt is supposed to be positive, and the friction contact is modeled by the Coulomb's law, i.e. F T = µF N where F T is the tangential friction force, F N is the normal force and µ is the friction coefficient. Using Newton's second law the equations of motion, that constitute the primary system in this work, are derived as follows:

• µ • µ c 1 c 2 k 2 , k N L 2 k 1 , k N L 1 x 2 x 1 m NES 1 NES 2 (a) NES 1 c h k h , k N L h h 1 (b) NES 2 c h k h , k N L h h 2 (c)
mẍ 1 + c 1 ẋ1 + k 1 x 1 -µk 2 x 2 + k N L 1 x 3 1 -µk N L 2 x 3 2 = 0 (1a) mẍ 2 + c 2 ẋ2 + k 2 x 2 + µk 1 x 1 + µk N L 1 x 3 1 + k N L 2 x 3 2 = 0. (1b)

Mechanical model with Nonlinear Energy Sinks

Two strongly cubic and identical ungrounded NES with masses m h , damping coefficients c h and cubic stiffnesses k N L h , are attached on the primary system (see Fig. 1). Introducing the NES displacements h 1 (t) and h 2 (t) in Eq. (1) leads to the following equations of motion of the coupled system

ẍ1 + η 1 ω 1 ẋ1 + ω 2 1 x 1 -µω 2 2 x 2 + φ 1 x 3 1 -µφ 2 x 3 2 + η h ω 1 ( ẋ1 -ḣ1 ) + ξ h (x 1 -h 1 ) + φ h (x 1 -h 1 ) 3 = 0 (2a) ϵ ḧ1 + η h ω 1 ( ḣ1 -ẋ1 ) + ξ h (h 1 -x 1 ) + φ h (h 1 -x 1 ) 3 = 0 (2b) ẍ2 + η 2 ω 2 ẋ2 + ω 2 2 x 2 + µω 2 1 x 1 + µφ 1 x 3 1 + φ 2 x 3 2 + η h ω 1 ( ẋ2 -ḣ2 ) + ξ h (x 2 -h 2 ) + φ h (x 2 -h 2 ) 3 = 0 (2c) ϵ ḧ2 + η h ω 1 ( ḣ2 -ẋ2 ) + ξ h (h 2 -x 2 ) + φ h (h 2 -x 2 ) 3 = 0, (2d)
where

η i = c i / √ mk i (with i = 1, 2) are the relative damping coefficients, ω i = √ k i /m are the natural angular frequen- cies, φ i = k N L i
/m; as for the NES parameters, ϵ = m h /m is the mass ratio between the NES and the primary structure,

ξ h = k h /m, η h = c h / √ mk 1 and φ h = k N L h /m.
Because strongly cubic NES are considered, one has ξ h ≪ φ h . In general, the mass ratio ϵ is in the range 0.01 -0.1 as in the present paper.

Possible steady-state regimes

Four main types of steady-state regimes (SSR) may be generated when one or several NES are attached on a primary system which undergoes an LCO: complete suppression of the instability, mitigation through Periodic Response (PR), mitigation through Strongly Modulated Response (SMR) or no mitigation. Those four regimes can also be observed on the studied system (2) [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF]), and are presented in Fig. 2 which shows for each case the displacements x 1 (t) -with and without NES attachments -with respect to time. The regimes described above are classified into two categories depending on whether the LCO is considered as mitigated (Figs. 2(a), 2(b) and 2(c)) or not (Fig. 2(d)).

In this work the Quantity of Interest (QoI) under consideration is the amplitude A wN ES x1 of the variable x 1 in the system with NES (2) and within a steady-state regime. It is defined as follows

A wNES 1 = max [x SSR 1 (t)] -min [x SSR 1 (t)] 2 , (3)
where x SSR 1 (t) is the times series of the variables x 1 obtained from the numerical integration of the coupled system (2) within the steady-state regime. For comparison purposes, the amplitude A woNES 1 of the system without NES (1) is also computed.

The time integration of system (2) (resp. ( 1)) is performed between t b = 0 and t e = 4 seconds and the amplitude A wNES 1 (resp. A woNES

1

) is computed in the interval [0.9t e t e ] in which we assume that the system has reached its steady-state regime. Using the following set of parameters

ω 1 = 2π100 (rad⋅s -1 ), ω 2 = 2π85 (rad⋅s -1 ), η 1 = 0.02, η 2 = 0.06, φ 1 = 10 5 (N⋅kg -1 ⋅m -3 ), φ 2 = 0 (N⋅kg -1 ⋅m -3 ), ϵ = 0.05, ξ h = 0.001 (N⋅kg -1 ⋅m -1 ), η h = 0.02, φ h = 1.4 ⋅ 10 5 (N⋅kg -1 ⋅m -3 ), (4)
the amplitudes A wNES . For increasing values of µ, this discontinuity corresponds in general to the transition from the SMR to the no suppression regime and separates mitigated regimes from unmitigated ones. The value of µ at the discontinuity is called mitigation limit and is denoted as µ ml .

The Propensity to be in a Harmless Steady-State Regime (PHSSR) [START_REF] Snoun | Prediction of the dynamic behavior of an uncertain friction system coupled to nonlinear energy sinks using a multi-element generalized polynomial chaos approach[END_REF]) can be computed from the mitigation limit µ ml . The PHSSR is the probability of obtaining a mitigated regime given the probability laws governing the dispersion of the primary system parameters. Here, only µ is considered as uncertain and it is assumed to follow a uniform probability law within the interval [0, µ e = 0.4]; the PHSSR is then simply defined as

PHSSR = ∫ µ ml 0 1 µ e dµ = µ ml µ e .
(5)

Multi-Element Polynomial Chaos based algorithm to detect a discontinuity

The ME-gPC method previously developed by the authors in [START_REF] Snoun | Prediction of the dynamic behavior of an uncertain friction system coupled to nonlinear energy sinks using a multi-element generalized polynomial chaos approach[END_REF] to locate the discontinuity in the QoI profile is recalled in Section 4.2. Elements of the Polynomial Chaos theory are first presented in Section 4.1.

Elements of the Polynomial Chaos theory

Generalized Polynomial Chaos (gPC)

Let us consider a dynamical system with r uncertain physical parameters β j (j = 1, . . . , r) assumed to be uniformly distributed within a given space ∏ r i=1 [a i , b i ], and let ξ j (j = 1, . . . , r) be the independent random variables within the space [-1 1] r , linked to β j by

β j (ξ j ) = a j + b j 2 + b j -a j 2 ξ j , (j = 1, . . . , r). (6)
The gPC theory [START_REF] Wiener | The homogeneous chaos[END_REF]; [START_REF] Cameron | The orthogonal development of nonlinear functionals in series of fourier-hermite functionals[END_REF]; [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF]) states that the system output variable X(ξ 1 , ..., ξ r ), called the quantity of Interest (QoI), may be approximated by a truncated orthogonal polynomial function series as

X(ξ 1 , ..., ξ r ) ≈ Np ∑ j=0 xj ϕ j (ξ 1 , ..., ξ r ), (7)
where ϕ j (ξ) are orthogonal polynomials which represent the stochastic part of the process, whereas xj are the gPC coefficients that take into account the deterministic part of the process. In this paper, according to the Askey-scheme [START_REF] Askey | Some basic hypergeometric orthogonal polynomials that generalize Jacobi polynomials[END_REF]), the orthogonal polynomials ϕ j (ξ 1 , ..., ξ r ) are built from Legendre polynomials because the uncertain parameters follow a uniform probabilistic law. From [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF], the number of terms N p + 1 of the expansion is given by

N p + 1 = (p + r)! p!r! , ( 8
)
where r is the number of uncertain parameters and p is the order of the gPC. To obtain the approximated values of the QoI X, the coefficients xj of the truncated series ( 7) must be computed. To that end, the non-intrusive regression method [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF]) is used in this paper.

In this method, the gPC coefficients are built from Q values of the QoI X. The corresponding Q simulations may be performed at points chosen with the Latin Hypercube Samples (LHS) method [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]) that will be referred to as the Numerical Experimental Design (NED) in this paper. A minimum of Q = N p simulations is required but in practice, Q = kN p simulations are used, with k a small integer usually equal to 2, 3 or 4.

4.1.2. Multi-Element Generalized Polynomial Chaos (ME-gPC) When the QoI is nonlinear, high polynomial orders may be required to reduce the approximation error. Thus, when the number of uncertain parameters r is high, the number of simulations required to compute the gPC coefficients may be excessive, leading to a prohibitive computational cost. An alternative is to split the vector ξ = (ξ 1 , ..., ξ r ) into a collection of m non-intersecting elements and to use a low order polynomial approximation on each element. That is the basic idea behind the ME-gPC [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF].

The local physical variables

β k j , in the k th element ∏ r i=1 [a k i , b k i ], are expressed in terms of independent uni- form random variables ξ k j in [-1, 1] r through β k j = b k j + a k j 2 + b k j -a k j 2 ξ k j . (9)
Consequently, Legendre polynomials can be used locally and a gPC expansion developped in each element as follows:

X k ( ξk ) ≈ Np ∑ j=0 xk,j ϕ j ( ξk ), (10)
where X k ( ξk ) is the random process corresponding to the k th element.

From the gPC theory, the approximated local mean Xp,k and variance σ 2 p,k in the k th element can be expressed analytically from the gPC coefficients as

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Xp,k = xk,0 σ 2 p,k = 1 2 r Np ∑ j=1 x2 k,j ⟨ϕ 2 j ⟩. (11)
Several criteria have been proposed in the literature for the convergence of the ME-gPC algorithm. They depend on the aim of the ME-gPC. The criteria chosen in this study are presented in the next section.

Detection of a discontinuity in the QoI with the ME-gPC

This section recalls the method based on the ME-gPC and proposed in [START_REF] Snoun | Prediction of the dynamic behavior of an uncertain friction system coupled to nonlinear energy sinks using a multi-element generalized polynomial chaos approach[END_REF] to locate a discontinuity in a Quantity of Interest (QoI) derived from the DM of a mechanical system. The aim of the method is not necessarily to obtain an accurate representation of the QoI, but to locate precisely the discontinuity, i.e. to know the values of the uncertain parameters for which the discontinuity in the QoI profile appears.

The underlying principle of the method is that the variance of a QoI that contains a discontinuity is higher than that of a QoI that does not contain any discontinuity. Therefore, if the discontinuity lies in a given element of the stochastic parameter space, the variance of the ME-gPC expansion in that element will be high. From these considerations a ME-gPC algorithm has been built, in which at each step and for each element the local variance σ 2 p,k , computed directly from the gPC coefficients using Eq. ( 11), is compared to a threshold θ 1 . If for a given element k

σ 2 p,k ≥ θ 1 , (12)
it is assumed that the discontinuity lies within this element. The latter is then divided by two in each direction of the stochastic space in order to locate more accurately the discontinuity. On the contrary, if the condition Eq. ( 12) is not satisfied, the element k is supposed not to contain the discontinuity and is removed from the algorithm. The threshold θ 1 is chosen so that the variance of the element which contains the discontinuity always satisfies the condition Eq. ( 12). Therefore, two other criteria are introduced to stop the algorithm. The first one is the minimum element size J min defined as

J min = θ 2 J 0 , (13)
where J 0 is the size of the initial element corresponding to the whole stochastic parameter space and θ 2 is a percentage. At a given iteration i, the retained elements of size J i are actually divided into two equal parts in each direction of the stochastic space if the following condition

J i ≥ J min (14)
holds. It is worth recalling that all the elements at iteration i have the same size J i because at each step the division is performed in each direction of the stochastic. If the condition ( 14) does not hold, the algorithm is stopped and we assume that the discontinuity in the QoI profile lies within the remaining elements.

In Section 6, a comparison is performed with a reference method that needs N simulations of the DM to locate a discontinuity in the QoI profile. Therefore, the proposed method is only effective if the number of numerical simulations N ′ required to obtain the gPC coefficients in all the elements is smaller than N . Consequently, the second criterion is that N ′ should verify

N ′ ≤ N, (15)
otherwise the algorithm is stopped.

The algorithm is summarized in Fig. 4. At the first step, only one element is present (i.e. the whole stochastic parameter space) and for a given iteration the parameter space has been divided into K elements during the previous iterations. At each iteration, the gPC coefficients are Consider the iteration in which the stochastic parameter space is divided into elements, (for the first iteration we have ).

i th K K = 1

-Build the PC expansion for the each elements ( ), ( simulations of the DM are needed to built the PC expansions). -Calculate the variance and the size of the elements.

k k = 1,…, K N′ K σ 2 p,k J i N′ ≤ N Consider the first element ( ). k = 1 σ 2 p,k ≥ θ 1 J i ≥ J min
Each of the elements is divided into 2 in each direction of the stochastic parameters space.

K′

and .

i = i + 1 K = 2 r K′ k = K STOP The discontinuity has been located.
Remove the element from the algorithm.

k th k = k + 1
We assume that the discontinuity in the QoI profile is in remaining elements.

K′ STOP

The discontinuity has not been located in general.

elements remain at the exit of the loop in .

K′ k

Figure 4: Algorithm of the proposed method based on the ME-gPC.

computed for the K elements and we check if the corresponding computational cost is larger than the reference one (see Eq. ( 15)). If that is the case, the algorithm is stopped and in general the discontinuity in the QoI profile has not been located. If Eq. ( 15) holds, we check the variance in each of the K elements (see Eq. ( 13)). If it is smaller than the threshold θ 1 , the element is removed from the algorithm. If it is larger, it is kept and K ′ elements remain at the exit of the loop in k. If Eq. ( 14) does not hold, the algorithm is stopped and we assume that the discontinuity in the QoI profile lies within the K ′ remaining elements. If Eq. ( 14) holds, then each of the K ′ elements is divided into two parts in each direction of the stochastic parameter space, then a new iteration begins with K = 2 r K ′ .

As for the reference method, the detection of the discontinuity can be performed by the so-called Jump Criterion (JC). To be considered as an effective jump, the difference between the amplitude at the right of the jump, denoted as A(µ + ml ), and the amplitude at the left, denoted as A(µ - ml ), must be greater than the half of A(µ - ml ), i.e.

JC = A(µ + ml ) -A(µ - ml ) A(µ - ml ) ⩾ 1 2 . ( 16
)
As an example, Fig. 5 shows the last iteration of the algorithm applied to the system (2) for θ 1 = 0.7⋅10 -3 , θ 2 = 1% and p = 1. The Reference curve is similar to that displayed in Fig. 3. One can see that the last interval obtained with the ME-gPC method contains the discontinuity observed on the reference.

Optimization strategy under uncertainties

This section focuses on the optimization methodology proposed to maximize the PHSSR defined in Section 3. Two NES parameters are retained as the design parameters to optimize: the damping coefficient η h and the nonlinear stiffness φ h . The uncertain parameter is again the friction coefficient µ. Following [START_REF] Boroson | Optimization under uncertainty of parallel nonlinear energy sinks[END_REF]; Pidaparthi and Missoum (2019), the problem is treated in two different ways depending on whether the physical design parameters η h and φ h are deterministic or stochastic.

Deterministic design variables

In this case, the parameters to optimize are directly the physical design parameters which are assumed to be known with certainty. The objective of the method is to find the optimal values of the NES parameters x d that maximize the PHSSR value. For this purpose, the PHSSR is determined in the space of uncertain parameters x u for each value of

x d ∈ [x (min) d x (max) d
]. The maximum of the PHSSR is then detected, which corresponds to the optimal value of x d . The method can be summarized as follows: ], called design intervals.

Maximize PHSSR (A(x d , x u )) Subject to x (min) d ⩽ x d ⩽ x (max) d where x d = {η h , φ h }, x u = {µ} (

Stochastic design variables

The physical design parameters are considered here to be random with a known probability law. This case is closer to a real situation, as it is known in particular that obtaining an exact value of a damping coefficient is almost impossible. The parameters to be optimized are then statistical characteristics of these physical parameters (for instance their means or variances) and are called hyperparameters. In this study the hyperparameters considered are the means ηh and φh of the damping coefficient and of the cubic stiffness of the NES, respectively.

The optimization process consists in dividing the design space into subspaces. In each subspace a draw of N values of the design parameters is performed according to their probability laws (in this work we use uniform laws). For each of these N values the PHSSR is calculated and the optimal value of x d (here of ηh or φh ) corresponds to the subspace where the mean of the N PHSSR values is maximum. This can be summarized as follows:

Maximize E[PHSSR (A(x d , x u ))] Subject to x (min) d ⩽ x d ⩽ x (max) d where x d = {η h , φh }, x u = {µ} (18)
where • E is the value of the mathematical expectation (the mean) of all PHSSR evaluations;

• ηh is the mean of the NES damping coefficient values;

• φh is the mean of the NES cubic stiffness values.

Results for deterministic design variables

As explained before, the first method of NES optimization considers that the NES parameters are deterministic and constitute therefore directly the design variables to optimize.

Optimization of the NES damping coefficient

In this section, we consider first that the design space is a one-dimensional space according to the variable η h , represented by the interval [0, 0.04]. The optimization problem ( 17) is therefore solved with x d = η h ∈ [0, 0.04] with the other parameters given Eq. ( 4). The results are expressed in terms of computational cost, on the one hand, and of precision compared to the reference method, on the other hand. For this purpose, the following quantities are introduced:

• µ max ml is the maximum value of the mitigation limit µ ml , and corresponds to the maximum value of the PHSSR, denoted as PHSSR max . The mitigation limit is here searched within the interval [0, 0.4]. The PHSSR value estimated by the ME-gPC method is determined using the upper bounds of the last intervals found by the algorithm presented in Fig. 5. The precision of the corresponding µ max ml value is therefore given by:

µ max ml ± ∆µ µ max ml . (19)
For the reference method, ∆µ is given by the discretization step of the design space, whereas ∆µ = 0.4θ 2 , whith θ 2 defined by Eq. ( 13), for the ME-gPC method.

• E PHSSR max is the relative error of the maximum PHSSR obtained with the ME-gPC method (PHSSR max ME-gPC ) compared to the maximum PHSSR value obtained with the reference method (PHSSR max ref ), and is defined as follows:

E PHSSR max = |PHSSR max ME-gPC -PHSSR max ref | PHSSR max ref . (20)
• η opt h is the optimal value of the NES damping coefficient, i.e. the value for which µ max ml , and thus PHSSR max , are obtained. Since the design space is discretized with a step ∆η h , the optimal value is expressed as follows

η opt h ± ∆η h η opt h . ( 21
)
• Err is the relative error of the optimal design value η opt h obtained with the ME-gPC method, compared to that obtained with the reference method:

Err = |(η opt h ) ME-gPC -(η opt h ) ref | (η opt h ) ref . (22)

Reference optimization

The reference value of the optimal NES damping coefficient, denoted as (η opt h ) ref , is obtained by solving (17) where the amplitude A(x d , x u ), and thus the mitigation limit, directly result from the numerical simulation of (2). The design space for η h [0, 0.04] is uniformly discretized into 51 values. To reduce the computational cost, the mitigation limit µ ml and the PHSSR are first roughly estimated for each η h value using only 1000 values of the friction coefficient µ ∈ [0, 0.4] and applying the criterion ( 16). Fig. 6 shows the resulting evolution of the PHSSR as a function of η h . This first step enables the locating of the interval -here [0.0208, 0.028] -in which the optimal value of the design parameter (η opt h ) ref lies. In this interval the PHSSR is maximum (PHSSR max = 50.10% ± 0.19%) but constant due to the too coarse discretization in µ. To refine this reference result, the PHSSR is then calculated in the design interval η h ∈ [0.02, 0.03] with the same discretization step for η h as previously, but now with 10 000 values of µ to determine the mitigation limit for each value of η h . The total number of simulations needed to obtain the reference value is therefore equal to 50 × 1000 + 10 × 10 000 = 150 000. Fig. 7 shows the new PHSSR as a function of η h . In this case, (µ max ml ) ref = 0.2004±0.02% which gives (PHSSR max ) ref = 50.1% ± 0.02%. The corresponding optimal design value of the damping coefficient is (η opt h ) ref = 0.0240 ± 3.33%.

Optimization with the ME-gPC method

In the proposed approach, the optimal value of the NES damping coefficient, denoted as (η opt h ) ME-gPC , is still obtained by solving (17), but the amplitude A(x d , x u ) and the mitigation limit ( 17) are now determined by the method based on the ME-gPC.

Three optimization procedures are applied with a gPC order p = 1, a variance threshold θ 1 = 2 ⋅ 10 -3 and three values of the element size threshold: θ 2 = 1%, 0.1% and 0.001%. The results of this ME-gPC based optimization are presented in Tab. 1, along with the reference results for comparison. Fig. 8 displays the variations of the PHSSR as functions of η h respectively obtained with the reference method and the three ME-gPC methods with θ 2 = 1%, 0.1% and 0.001%. For θ 2 = 1% the number of simulations is equal to 4324, the optimal value is (η opt h ) ME-gPC = 0.0241 ± 15.14% and the maximum value of the corresponding PHSSR is 50.78%, implying an error compared to the reference solution E PHSSR max = 1.35%. In the case of θ 2 = 0.1%, although the new value of θ 2 has no influence on η opt h , it changes the maximum value of the PHSSR to PHSSR max ME-gPC = 50.09% which induces a relative error E PHSSR max = 0.02%. This error is lower than the previous one (with θ 2 = 1%) but the number of simulations has increased and is now equal to 5165. For θ 2 = 0.001% the optimal value becomes (η opt h ) ME-gPC = 0.0237 ± 3.33%, implying a relative error Err = 1.25% compared to the reference value, which is higher than the errors in the previous cases. This is not illogical since the error is calculated according to the midpoints of the intervals. Nevertheless, the interval obtained with θ 2 = 0.001% is much narrower (±3.33%) than that obtained with θ 2 = 0.01% (±15.14%), and contains the value resulting from the reference method. The total number of simulations is in this last case equal to 7026, which remains 22 times smaller than the number needed with the reference method.

The plateaus observed in Fig. 8 are due to the fact that for two close values of η h the mitigation limits are found in the same element. They are related to the value of the threshold θ 2 : the lower it is, the more accurate the localization of the mitigation limit is, with in counterpart an increase of the computational cost. The precision is therefore driven by the value of the threshold θ 2 .

Finally, with a computational cost much lower than the reference one, the proposed approach proves efficient to identify the value of the optimal NES damping coefficient while taking into account the dispersion of the friction coefficient.

Optimization of the NES nonlinear stiffness

The aim of this section is to find the optimal value of the cubic stiffness of the NES φ h assumed here to be deterministic. The damping coefficient value η h is here set to 0.024. Similarly to the previous section, the following quantities are defined:

• φ opt h is the optimal value of the NES cubic stiffness for which µ max ml , and thus PHSSR max , are obtained. It may be again expressed with respect to the discretization step ∆φ h of the design space as

φ opt h ± ∆φ h φ opt h . ( 23
)
• Err is now the relative error of the optimal design value φ opt h obtained with the ME-gPC method, compared to that obtained with the reference method:

Err = |(φ opt h ) ME-gPC -(φ opt h ) ref | (φ mid h ) ref .
(24)

Reference optimization

The same strategy has been applied for the reference optimization of the NES cubic stiffness φ h . First, to get an idea of the shape of the PHSSR as a function of φ h , a small number of simulations is performed with a discretization step equal to 10 on the design space φ h ∈ [0, 900 000] N⋅kg -1 ⋅m -3 . For each value of φ h the PHSSR is computed using 1000 values of µ within the interval [0, 0.4]. Fig. 9 shows the resulting evolution of the PHSSR as a function of φ h . The mitigation limit appears to decrease monotonically with the cubic stiffness value. The optimal value φ opt h is therefore the lowest value of φ h for which the jump criterion ( 16) is met and is found equal to 100 000 N⋅kg -1 m -3 . The next step then consists in focusing on a reduced design space composed of the lowest values of φ h in order to understand the dynamic behavior of the system for these values and to refine the optimization result. The retained design interval φ h ∈ [0, 150 000] N⋅kg -1 ⋅m -3 is uniformly discretized into 100 values, and for each one the mitigation limit is directly computed from (2) using again 1000 values of µ within [0, 0.4] interval. The number of simulations is thus equal to 100 × 1000 = 100 000.

Fig. 10 shows the evolutions of the jump criterion ( 16) and of the PHSSR as functions of φ h , respectively. For low values of φ h (φ h ∈ [0 46 970] N⋅kg -1 ⋅m -3 a jump exists but, as mentioned before, we consider that the difference between the amplitudes of the limit cycles before and after the jump is too small for the concerned regimes to be qualified as mitigated. The difficulty of this optimization problem is then to find the lowest value of φ h for which the discontinuity jump is high enough.

It is clear that the solution depends on the value of the threshold for the jump criterion (here 50%). The optimal value here is therefore (φ opt h ) ref = 46 970 ± 3.22% N⋅kg -1 ⋅m -3 , which corresponds to the maximum PHSSR value PHSSR max = 52.25% ± 0.19%. The PHSSR curve exhibits again stationary phases linked to the discretization step for µ. Indeed, for two close values of φ h , the discontinuities in the amplitude profiles as functions of µ appear for very close µ values whose difference is smaller than the discretization step.

Optimization with the ME-gPC method

We denote by (φ opt h ) ME-gPC the optimal value of the NES cubic stiffness coefficient obtained by solving (17) where the amplitude A(x d , x u ) and the mitigation limit are assessed using the ME-gPC based method.

The gPC order is again set to p = 1 and the mitigation limit, as explained before, is approached by the up-per bound of the last element found. Tab. 2 compares the results obtained with this ME-gPC method, for several thresholds values θ 1 and θ 2 , to those arising from the reference procedure (in which, as detailed in the previous section, the values of φ h are taken from the design space [0, 150 000] N⋅kg -1 ⋅m -3 which is uniformly discretized).

Fig. 11 shows the evolution of the PHSSR as a function of φ h for θ 1 = 3.5 ⋅ 10 -3 % and θ 2 = 0.1%.

For a given value of the threshold θ 2 , it can be observed from Tab. 2 that the precision of the results depends on the jump threshold θ 1 . Indeed, when the value of θ 1 increases, the error PHSSR max decreases from 3.15% to 1.68% and the error Err decreases from 29.03% to 3.22%. Another benefit of increasing θ 1 is that it reduces the number of required simulations, as can be seen in Tab. 2. However, for higher values of θ 1 , the ME-gPC method appears unable to detect the discontinuity in the amplitude profile, as the jump in the amplitude values becomes lower than the threshold θ 1 , which makes the research of the opti- Figure 9: PHSSR as a function of φ h for η h = 0.024 and φ h ∈ [0 900 000] N⋅kg -1 ⋅m -3 . The mitigation limits are assessed using the reference method.

mum θ 1 value difficult. Another way to reduce the error levels is to increase the value of the precision threshold θ 2 . The optimum value of the jump detection threshold θ 1 depending on the precision threshold θ 2 , the optimum pair of thresholds is finally found as θ 1 = 3.5 ⋅ 10 -3 and θ 2 = 0.1 %. The corresponding relative errors are 0.19 % for PHSSR max and Err = 0 % for φ opt h , while the number of simulations has increased from 5706 to 7023. It is therefore observed that the precision level depends on the accepted computational cost. In this study, the final computational cost (7 023) with the ME-gPC method remains much weaker than that of the reference method (100 000). Remark: as in section 6.1 for the optimization of the NES JC PHSSR coefficient, identical values of PHSSR max may observed with the ME-gPC method. The widths of the plateaus decrease when the precision threshold θ 1 is reduced.

Results for stochastic design variables

The second optimization technique considers that the NES parameters are also uncertain with a known probability law. Thus, the design variables to be optimized are not directly the NES parameters but one of their statistics called hyperparameters. In this work, the design parameters follow uniform probability distributions and the considered hyperparameter is the mean xd . Regarding the Figure 11: PHSSR as a function of φ h for η h = 0.024 and φ h ∈ [0 150 000] N⋅kg -1 ⋅m -3 . The thresholds are respectively: θ 1 = 3.5 ⋅ 10 -3 and θ 2 = 0.1%. The reference results are obtained as described in Sect. 6.2.2. system parameters, as previously only the friction coefficient µ is assumed to be uncertain.

Optimization of the NES damping coefficient

Reference optimization

In this section, the optimization problem ( 18) is solved with a single design variable (hyperparameter) defined as the mean ηh of the NES damping coefficient η h , which is a stochastic parameter following a uniform law. The NES cubic stiffness is fixed and equal to φ h = 140 000 N⋅kg -1 ⋅m -3 . The design space is divided into N I = 5 intervals as shown in Tab. 3. For each interval the lower and upper bounds are determined as follows: η inf h = ηh -0.1η h and η sup h = ηh + 0.1η h where ηh is the mean. The mean values are chosen so that the upper bound of one interval coincides with the lower bound of the next interval.

In each interval, N d = 10 random samples of damping coefficient values η h are drawn according to the uniform probability law. For each of these 10 samples, the mitigation limit and the PHSSR are determined directly from the numerical simulation of (2) using 10000 values of the friction coefficient µ. The objective function to be maximized is now the mean E[PHSSR] of the 10 PHSSR values obtained in each interval.

Tab. 3 shows the reference results of the ηh optimization. For each interval, the mean of the 10 samples of η h is calculated and compared to the theoretical ηh value in order to validate the number of samples. The mean PHSSR value E[PHSSR] is also indicated, along with the error on E[PHSSR] for a confidence level of 99%. The resulting values (errors lower than 0.15%) show that the number of samples N d = 10 of η h in each interval is sufficient.

Fig. 12 displays the evolutions of E[PHSSR], PHSSR max and PHSSR min as functions of ηh . In this case, the maxi-mum value of the objective function is 50.07%, which corresponds to an optimal theoretical value of ηh equal to 0.0237. The number of simulations is equal to 5 × 10 × 10 000 = 500 000.

Optimization with the ME-gPC method

Tab. 4 shows the results of the optimization of ηh by the ME-gpC based method. The gPC order is again chosen as p = 1, θ 1 is set to 1 ⋅ 10 -3 and two values of θ 2 are tested: θ 2 = 1% and θ 2 = 0.001%. The design space is the same as for the reference optimization and the same quantities are considered.

Tab. 5 compares the results obtained with the three strategies (reference optimization and ME-gPC based optimizations with θ 2 = 1% and θ 2 = 0.001%). The corresponding evolutions of E[PHSSR], PHSSR max and PHSSR min according to ηh are displayed in Fig. 12.

In the case of θ 2 = 1%, the number of simulations is equal to 4613, leading to a computational cost reduction of 98.98%. The maximum value of the objective function is here 50.78%, which corresponds to the same optimal value ηh = 0.0237 as for the reference case. The relative error of the maximum PHSSR mean value with respect to the reference is E E[PHSSR] max = 1.41%.

For θ 2 = 0.001%, the number of simulations is increased to 8006, which still induces a high computational cost reduction of 98.38%. The maximum value of the objective function is equal to that of the reference case, 50.07%, leading to a relative error E E[PHSSR] max = 0%. The corresponding optimal value of ηh is again equal to 0.0237.

As for the optimization with deterministic design variables, it is again shown that a proper choice of the precision threshold θ 2 enables the reduction of the error levels. However, an error reduction leads to an increase of the required simulation number (from 4613 to 8006), and thus the retained precision level depends on the accepted computational cost. In both cases presented above, the computational costs with the ME-gPC optimization remains much weaker than that of the reference method (500 000). Moreover, it can be noticed that the optimal values of the mean of η h , ηh = 0.0237, and of the PHSSR mean, E(PHSSR) = 50, 07%, are the same as those obtained with the deterministic optimization. This result is logical because the curve of the PHSSR as a function of η h varies very little around the optimal value. Thus, the PHSSR mean remains practically constant.

Optimization of the NES nonlinear stiffness

Reference optimization

The last optimization problem proposed in this paper (still represented by ( 18)) considers a single design variable (hyperparameter) φh which is here the mean of the NES nonlinear stiffness φ h ; similarly to the previous optimization problem, φ h is a stochastic parameter supposed to follow a uniform law. The NES damping coefficient η h is set to 0.024. The methodology is similar to that of the previous section with a design space divided into N I = 5 intervals as shown in Tab. 6. The lower and upper bounds of the intervals are again defined as φ inf h = φh -0.1 φh and φ sup h = φh + 0.1 φh , where the mean value φh is chosen so that the upper bound of one interval coincides with the lower bound of the next interval.

As previously, N d = 10 random samples of stiffness values φ h are drawn in each interval according to the uniform probability law, and for each sample the mitigation limit and the PHSSR are directly determined from the numerical simulation of (2) using 10000 values of the friction coefficient µ. The objective function to be maximized is the mean E[PHSSR] of the 10 PHSSR values obtained in each interval.

In section 6.2 it has been shown that the PHSSR is a decreasing function of φ h , which implies that the optimal value is the lowest value of φ h for which the jump criterion (16) is met. It is therefore possible that a part of the interval where this optimal value lies (in concrete terms here, some of the 10 samples in the first interval [38 438, 46 980] N⋅kg -1 ⋅m -3 ) corresponds to φ h values for which the jump criterion ( 16) is not respected. In that case, 10 new samples of φ h values are generated in the part of the interval Tab. 6 shows the reference results of the φh optimization. Similarly to the previous section, the mean of the 10 samples φ h is calculated for each interval and compared to the theoretical φh value. The mean value E[PHSSR] is also indicated, along with the error on E[PHSSR] for a confidence level of 99%. With very low error levels (lower than 0.03%), the number of φ h samples N d = 10 in each interval is again found sufficient. This reference optimization has required a simulation number equal to 5 × 10 × 10 000 = 500 000.

The evolutions of E[PHSSR], PHSSR max and PHSSR min according to φh are displayed in Fig. 13. In this case, the maximum value of the objective function is 52.32% which corresponds to an optimal value of φh equal to 46678. The very low error on the E(PHSSR) for a 99% confidence level obtained in this case is due to the fact that the PHSSR values are close in the interval.

Optimization with the ME-gPC method

Tab. 7 shows the results of the optimization of φh based on the ME-gPC, with again a gPC order p = 1, θ 1 = 4.8 ⋅ 10 -3 and two values of θ 2 : θ 2 = 1% and θ 2 = 0.001%. The design space is discretized as for the reference optimization. For several intervals the error on E[PHSSR] for a confidence level of 99% is null, meaning that the 10 PHSSR values corresponding to the 10 samples are identical. The evolutions of E[PHSSR], PHSSR max and PHSSR min according to φh for the cases θ 2 = 1% and θ 2 = 0.001% are visible in Fig. 13 along with the reference results.

The comparison of the different strategies to optimize the NES nonlinear stiffness (that is, the reference optimization and the two ME-gPC based optimizations with θ 2 = 1% and θ 2 = 0.001%) is presented in Tab. 8. In the case of θ 1 = 4.8 ⋅ 10 -3 and θ 2 = 1%, the number of simulations is equal to 3248, leading to a computational cost reduction of 99.21%. The maximum value of the objective function is 52.73%, which corresponds to an optimal value of φh equal to 46 678 N⋅kg -1 ⋅m -3 , and the relative error on the maximum value of the PHSSR mean is

E E[PHSSR] max = 1.5%.
When θ 1 = 2 ⋅ 10 -3 and θ 2 = 0.001%, the number of simulations is increased up to 6947, inducing a slightly reduced computational cost reduction of 98.30%. The maximum value of the objective function is now 52.33%, the optimal value of φh 46 678 N⋅kg -1 ⋅m -3 , and the relative error on the maximum value of the PHSSR mean E E[PHSSR] max = 0.09%.

The conclusions given for the optimization of the damping coefficient with the ME-gPC based method hold here: the choice of the precision threshold θ 2 enables the control of the error levels but impacts the computational costs, which remain however much lower than that with the reference method (500 000 simulations).

Conclusion

In this paper, the robust optimization of NES used to mitigate limit cycles of a friction system with an uncertain friction coefficient has been studied. The optimization strategy is based on the maximization of the propensity of the system to be in an attenuated regime (PHSSR), which constitutes the objective function of the problem. Two original approaches are proposed, depending on whether the parameters of the NES are considered as deterministic or as uncertain, but with a known probability distribution. In each case, two optimization strategies are developed, that consist in evaluating the objective function by either a deterministic reference optimization or by a method based on the Multi-Element generalized Polynomial Chaos (the so-called ME-gPC method). The comparison between the two methods shows the efficiency of the ME-gPC method to significantly reduce the computational cost while keeping a good accuracy compared to the reference optimization.

These optimization approaches have been applied to determine the optimal values of the NES parameters to maximize their ability to passively attenuate limit cycle of oscillations. Concerning the NES nonlinear stiffness, the optimal value is the lowest value for which the discontinuity jump is considered as high enough. Concerning , PHSSR max and PHSSR min as functions of φh with the reference method, the ME-gPC based method with θ 2 = 1 % and the ME-gPC based method with θ 2 = 0.001 %. The black segments denote the limits of the intervals. The reference results are obtained as described in Sect. 7.2.1. the NES damping coefficient, the curve representing the PHSSR as a function of the damping coefficient is nonlinear and presents a maximum corresponding to the optimal NES damping value. However, the variations of the PHSSR near the optimal value are very small, which makes it more difficult (but also less crucial) to precisely locate this optimal value. The general nature of the proposed strategy suggests that it could be used for the robust optimization of any vibration absorber, providing that the transition from the attenuated state to the non-attenuated state passes through a discontinuity in the amplitude profile of the system as a function of the bifurcation parameter considered.

In future works, several issues should be tackled. First, as one of the difficulties of the proposed approach is to choose relevant values for the thresholds, it could be interesting to develop automatically adaptive thresholds. Another issue lies in the extension of the proposed method to several uncertain parameters. It can also be envisaged to study the efficiency of the method for a model with a higher number of degrees of freedom.

Figure 1 :

 1 Figure 1: The mechanical model. (a) Mechanical system with NES; (b) Zoom on NES 1 ; (c) Zoom on NES 2 .

Figure 2 :

 2 Figure 2: Comparison between time series x 1 (t) resulting from the numerical integration of the friction system with and without NES. (a) Complete suppression, µ = 0.16; (b) Mitigation: PR, µ = 0.18; (c) Mitigation: SMR, µ = 0.2; (d) No mitigation, µ = 0.22. The set of parameters Eq. (4) has been used.

  functions of the friction coefficient µ in Fig. 3. The specific values of µ denoted as µ wo b and µ w b correspond to the Hopf bifurcation points without and with the NES attachments, respectively. The four steady-state regimes described above are visible on the Figure, which highlights a discontinuity (or jump) in the profile of A wNES 1

Figure

  Figure 3: Amplitudes A wNES 1

Figure 5 :

 5 Figure5: Last iteration of the detection algorithm when the friction coefficient µ is uncertain, with θ 1 = 0.7 ⋅ 10 -3 , θ 2 = 1% and p = 1; reference curve obtained from Fig.3.

  x d , x u ) is the amplitude of the displacement of the coupled system (2) defined by Eq. (3) (from now on the amplitude is simply written A instead of A wNES 1 to lighten the notations); it depends on both the uncertain parameters x u and the design variables x d ; and maximum values of the design variables x d respectively; the optimal values of the design parameters are searched within the intervals [x (min) d x (max) d

Figure 6 :

 6 Figure 6: Result of the first step of the reference optimization procedure: PHSSR as a function of η h , for φ h = 140 000 N⋅kg -1 ⋅m -3 .

Figure 7 :

 7 Figure 7: Result of the second step of Reference optimization: PHSSR as a function of η h , for φ h = 140 000 N⋅kg -1 ⋅m -3 .

Figure 8 :

 8 Figure 8: Comparison between the optimization strategies based on the reference method and the ME-gPC for η h ∈ [0 0.04] with φ h = 140 000 N⋅kg -1 ⋅m -3 . (a) PHSSR and (b) Zoom of the PHSSR. The reference results are obtained as described in Sect. 6.1.1.

Figure 10 :

 10 Figure 10: Evolutions of the Jump Criterion (JC) and of the PHSSR as functions of φ h obtained with the reference method.

Figure 12 :

 12 Figure12: Evolutions of PHSSR max and PHSSR min as functions of ηh obtained with the reference method, the ME-gPC based method with θ 2 = 1 %, and the ME-gPC based method with θ 2 = 0.001 %. The black segments denote the limits of the intervals. The reference results are obtained as described in Sect. 7.1.1.

Figure 13 :

 13 Figure13: Evolutions of E[PHSSR], PHSSR max and PHSSR min as functions of φh with the reference method, the ME-gPC based method with θ 2 = 1 % and the ME-gPC based method with θ 2 = 0.001 %. The black segments denote the limits of the intervals. The reference results are obtained as described in Sect. 7.2.1.

Table 1 :

 1 Comparison between the reference optimization and the ME-gPC optimization considering deterministic design variables. The design parameter is η h . The reference results are obtained as described in Sect. 6.1.1.

		Nb of sim-ulations	Computational time	PHSSR max E PHSSR max	η opt h	Err
	Reference	150 000	2.15 days	50.10% ± 0.02%	-	0.0240 ± 3.33%	-
	ME-gPC (θ1 = 2 ⋅ 10 -3 and θ2 = 1%)	4324	1.70 hours	50.78% ± 1.97%	1.35%	0.0241 ± 15.14%	0.41%
	ME-gPC (θ1 = 2 ⋅ 10 -3 and θ2 = 0.1%)	5165	2.06 hours	50.09% ± 0.2%	0.02%	0.0241 ± 15.14%	0.41%
	ME-gPC (θ1 = 2 ⋅ 10 -3 and	7026	3.23 hours	50.09% ± 0.02%	0.02%	0.0237 ± 3.33%	1.25%
	θ2 = 0.001%)

Table 2 :

 2 Comparison between the reference optimization and the Me-gPC based optimization considering deterministic design variables. The design parameter is φ h . The reference results are obtained as described in Sect. 6.2.2.

		Nb of simulations	Computational time	PHSSR max E PHSSR max	φ opt h	Err
	Reference	100 000	1.49 days	52.25% ± 0.19%	-	46 970 ± 3.22%	-
	ME-gPC (θ1 = 3.3 ⋅ 10 -3 and θ2 = 1%)	6 652	2.64 hours	53.90% ± 1.86%	3.15%	33 333 ± 4.50%	29.03%
	ME-gPC (θ1 = 4.2 ⋅ 10 -3 and θ2 = 1%)	5 817	2.31 hours	53.13% ± 1.88%	1.68%	39 394 ± 3.40%	22.51%
	ME-gPC (θ1 = 4.8 ⋅ 10 -3 and θ2 = 1%)	5 706	2.25 hours	53.13% ± 1.88%	1.68%	45 454 ± 3.33%	3.22%
	ME-gPC (θ1 = 3.5 ⋅ 10 -3 and θ2 = 0.1%)	7 023	2.96 hours	52.34% ± 0.19%	0.35%	46 970 ± 3.22%%	0%

Table 3 :

 3 Results of the reference optimization (see Sect. 7.1.1). The design parameter is ηh .

			[η inf h , η sup h ]	Theoretical mean ηh	Mean of the 10 samples of η h	E[PHSSR] (%)	Error on E[PHSSR] for a confidence level = 99% (%)
			[0.0117,0.0143]	0.0130	0.0129	48.36	0.1417
			[0.0143,0.0175]	0.0159	0.0161	49.31	0.1086
			[0.0175,0.0214]	0.0194	0.0193	48.85	0.0220
			[0.0214,0.0261]	0.0237	0.0234	50.07	0.0005
			[0.0261,0.0319]	0.0290	0.0284	49.90	0.0277
	0.010	0.015	0.020	0.025	0.030

Table 4 :

 4 Results of the optimization with the ME-gPC method. The design parameter is ηh .

	θ2 = 1 %	θ2 = 0.001 %

Table 5 :

 5 Comparison between the reference optimization and the ME-gPC based optimization. The design parameter is ηh . The reference results are obtained as described in Sect. 7.1.1.

		Nb of sim-ulations	Computational time	Computational reduction(%) cost	E[PHSSR] max E E[PHSSR] max	ηopt h
	Reference	500 000	7.86 days	-	50.07% ± 0.02%	-	0.0237
	ME-gPC (θ1 = 1 ⋅ 10 -3	4 613	1.92 hours	98.98	50.78% ± 2%	1.41%	0.0237
	and θ2 = 1%)
	ME-gPC (θ1 = 1 ⋅ 10 -3 and	8 006	3.29 hours	98.38	50.07% ± 0.02%	0%	0.0237
	θ2 = 0.001%)
	(here [38 438, 46 980]) for which (16) is satisfied.

Table 6 :

 6 Results of the reference optimization (see Sect. 7.2.1). The design parameter is φh .

					Theoretical mean φh	Mean of 10 samples of φ h	E[PHSSR] (%)	Error on E[PHSSR] for a confidence level = 99% (%)
			[38 438, 46 980]	42 709	46 678	52.32	0.0004
			[46 980, 57 420]	52 200	52 992	51.95	0.0299
			[57 420, 70 180]	63 800	62 746	51.53	0.0296
			[70 180, 85 776]	77 978	78 942	51.07	0.0095
			[85 776, 104 837]	95 306	94 876	50.76	0.0200
	40 000	60 000	80 000	100 000

Table 7 :

 7 Results of the optimization with the ME-gPC method. The design parameter is φh .

	θ2 = 1 %	θ2 = 0.001 %

Table 8 :

 8 Comparison between the reference optimization and the optimization with the ME-gPC method. The design parameter is φh . The reference results are obtained as described in Sect. 7.2.1.E[PHSSR]max E E[PHSSR] max φopt

				Computational
		Nb of sim-	Computational	cost
		ulations	time	reduction			h
				(%)
	Reference	500 000	7.53 days	-	52.32% ± 0.02%	-	46 678
	ME-gPC (θ1 = 4.8 ⋅ 10 -3	3248	1.42 hours	99.21	52.73% ± 1.9%	1.50%	46 678
	and θ2 = 1%)
	ME-gPC (θ1 = 2 ⋅ 10 -3 and	6947	3.06 hours	98.30	52.33% ± 0.019%	0.09%	46 678
	θ2 = 0.001%)