Evangéline Capiez-Lernout
email: evangeline.capiez-lernout@univ-eiffel.fr

Christian Soize
email: christiansoize@univ-eiffel.fr

Nonlinear stochastic dynamics of detuned bladed-disks with uncertain mistuning and detuning optimization using a probabilistic machine learning tool

Keywords: nonlinear stochastic dynamics, bladed-disk, detuning, mistuning, reduced-order model, uncertainty quantification, combinatorial optimization, machine learning, probabilistic learning, PLoM HFCM from w c, . q identity matrix in M n degree of freedom. pdf: probability density function. ANN: artificial neural network. FFT: Fast Fourier Transform. HFCM: high-fidelity computational model. KDE: kernel density estimation. MCMC: Markov chain Monte Carlo. NL-ROM: nonlinear reduced-order model. NL-SROM: nonlinear stochastic reduced-order model. PCA: principal component analysis. PLoM: probabilistic learning on manifolds. QoI: quantity of interest. ROB: reduced-order diffusion-maps basis. 2. Stochastic computational model formulation, methodology, and algorithms

The paper deals with the nonlinear stochastic dynamics concerning the detuning optimization in presence of random mistuning of bladed-disks with geometrical nonlinearities. We present an efficient computational methodology for reducing the computational cost, an analysis of the detuning, and the detuning optimization, based on the use of a high-fidelity computational model. A deep computational analysis is presented for a 12-bladed-disk structure that is representative of industrial turbomachines in order to understand the role played by the geometrical nonlinearities on the dynamical behavior and to exhibit the consequences on the detuning effects. For the detuning optimization with a very large number of possible detuned configurations, we propose a reformulation of the combinatorial optimization problem in a probabilistic framework, which is adapted to a probabilistic machine learning tool in order to limit the number of evaluations of the cost function with the high-fidelity computational model. The methodology proposed is validated for the 12-bladed-disk structure for which the exact optimal detuned configuration has been identified. A very good prediction is obtained.

Introduction

The vibrational behavior of turbomachines is known to be particularly complex, requiring the construction of predictive computational models that have also to be efficient in terms of numerical costs. One of the issues concerns the mistuning caused by the small variations of the mechanical properties from one sector to another one, induced by the manufacturing tolerances or by the wear and tear of the structure. Such phenomenon can generate strong localization effects yielding larger dynamical amplifications of the forced response with respect to the perfect cyclic symmetry case [START_REF] Whitehead | Effects of mistuning on the vibration of turbomachine blades induced by wakes[END_REF][START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF][START_REF] Wei | Localization phenomena in mistuned assemblies for cyclic symmetry -part ii: forced vibrations[END_REF]. This amplification strongly affects the fatigue life of the blades, which can cause safety problems and impair the proper operation of turbomachines. Many research have been carried out on this subject, requiring the use of probabilistic approaches for modeling the random character of the mistuning combined to the construction of reduced-order models [START_REF] Bladh | Component-mode-based reduced order modeling techniques for mistuned bladed disks. part 1: theoretical models[END_REF][START_REF] Yang | A reduced-order model of mistuning using a subset of nominal modes[END_REF][START_REF] Capiez-Lernout | Nonparametric modeling of random uncertainties for dynamic response of mistuned bladed-disks[END_REF][START_REF] Sinha | Reduced-order model of a bladed rotor with geometric mistuning[END_REF][START_REF] Mbaye | A reduced-order model of detuned cyclic dynamical systems with geometric modifications using a basis of cyclic modes[END_REF][START_REF] Vargiu | A reduced order model based on sector mistuning for the dynamic analysis of mistuned bladed disks[END_REF][START_REF] Madden | Reduced-order modeling approach for blisks with large mass, stiffness, and geometric mistuning[END_REF][START_REF] Schwerdt | A model reduction method for bladed disks with large geometric mistuning using a partially reduced intermediate system model[END_REF] when linear assumption is considered. There also exist other phenomena that can strongly affect the vibrational behavior of the bladed-disk and that have to be the subject of a dedicated modeling. The presence of nonlinearities increases the complexity of the forced response analysis because couplings and energy transfers can strongly modify the position and the nature of the resonances of the bladed-disk. Contact nonlinearities occur for fan-type bladed-disks for which the disk and the blades are manufactured independently and for which the blade assembly is mounted on the disk with "fir-tree" or dovetail fixings. Such local nonlinearities although occur through rotor/stator interaction between the blade tip and the stator casing or when using under-platform dampers yielding blade to blade contact (see e.g. [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF][START_REF] Krack | Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact[END_REF][START_REF] Batailly | Full three-dimensional rotor/stator interaction simulations in aircraft engines with time-dependent angular speed[END_REF][START_REF] Pesaresi | Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications[END_REF][START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF][START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF][START_REF] Quaegebeur | Nonlinear dynamic analysis of three-dimensional bladed-disks with frictional contact interfaces based on cyclic reduction strategies[END_REF]). A research effort has also been made for the construction of nonlinear reduced-order models allowing for the mistuning analysis of bladed-disk in presence of contact nonlinearities [START_REF] Mehrdad Pourkiaee | A reduced order model for nonlinear dynamics of mistuned bladed disks with shroud friction contacts[END_REF][START_REF] Koscso | Sensitivity and forced response analysis of anisotropy-mistuned bladed disks with nonlinear contact interfaces[END_REF][START_REF] Chen | Analysis of nonlinear modal damping due to friction at blade roots in mistuned bladed disks[END_REF][START_REF] Joachim | Numerical investigation of a mistuned academic bladed disk dynamics with blade/casing contact[END_REF][START_REF] Nyssen | Investigation on the robustness of rotor/stator contact interactions with respect to small mistuning[END_REF][START_REF] Quaegebeur | Impact of mistuned underplatform dampers on the nonlinear vibration of bladed disks[END_REF][START_REF] Pourkiaee | Relative cyclic component mode synthesis: A reduced order modeling approach for mistuned bladed disks with friction interfaces[END_REF]. Moreover, due to technological evolutions that involve lighter materials and thinner blades, the nonlinear geometrical effects induced by finite displacements cannot longer be neglected. The main difficulties concerning the modeling of these nonlinear geometrical effects concern the construction of adapted nonlinear reduced-order models that have to be efficient not only in terms of computational costs but also in terms of predictability. These aspects have also been investigated in the more general context of thin and slender structures [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF][START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF][START_REF] Kim | Nonlinear reduced order modeling of isotropic and functionally graded plates[END_REF][START_REF] Wang | Locally enhanced reduced order modeling for the nonlinear geometric response of structures with defects[END_REF][START_REF] Wang | Structural uncertainty modeling for nonlinear geometric response using nonintrusive reduced order models[END_REF][START_REF] Wang | Nonlinear reduced-order modeling of flat cantilevered structures: Identification challenges and remedies[END_REF]. Concerning the context of turbomachines, various research can be found for the deterministic case [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Grolet | Computing multiple periodic solutions of nonlinear vibration problems using the harmonic balance method and groebner bases[END_REF][START_REF] Martin | Dynamic analysis and reduction of a cyclic symmetric system subjected to geometric nonlinearities[END_REF], but also in presence of contact nonlinearities [START_REF] Delhez | Reduced order modeling of blades with geometric nonlinearities and contact interactions[END_REF] or in presence of mistuning [START_REF] Capiez-Lernout | Mistuning analysis and uncertainty quantification of an industrial bladed disk with geometrical nonlinearity[END_REF].

The intentional mistuning, also called detuning, consists in voluntarily breaking the cyclic symmetry of the structure by using partial or alternating patterns of different sector types. The detuning allows for spreading the frequencies of adjacent blades and thus for reducing the interaction between them. As a consequence, the drastic amplification and localization effects caused by the random mistuning can be greatly reduced. Such technology is particularly relevant because it is a way for reducing the response amplification levels that are induced by the unavoidable random mistuning while increasing its robustness. It has been thoroughly studied in the framework of linear dynamics [START_REF] Castanier | Investigation of the combined effects of intentional and random mistuning on the forced response of bladed disks[END_REF][START_REF] Castanier | Using intentional mistuning in the design of turbomachinery rotors[END_REF][START_REF] Choi | Optimization of intentional mistuning patterns for the reduction of the forced response effects of unintentional mistuning: formulation and assessment[END_REF][START_REF] Mbaye | Robust analysis of design in vibration of turbomachines[END_REF][START_REF] Han | Optimization of intentional mistuning patterns for the mitigation of the effects of random mistuning[END_REF][START_REF] Beirow | Vibration analysis of an axial turbine blisk with optimized intentional mistuning pattern[END_REF][START_REF] Lupini | Mistuning identification technique based on blisk detuning[END_REF][START_REF] Nakos | Mistuning and damping of a radial turbine wheel. part 1: fundamental analyses and design of intentional mistuning pattern[END_REF][START_REF] Sanchez-Alvarez | Key action mechanisms of intentional mistuning[END_REF][START_REF] Beirow | Vibration reduction of a steam turbine wheel by means of intentional mistuning[END_REF].

In a context of sustainable development, technological innovations have to be compatible with energy and environmental issues. This gives rise to the development and to the use of lighter materials, that have to be capable of high thermo-mechanical resistance and of long durability, so that it may be possible that the blade displacements respond in its nonlinear vibration range. In this context, a computational methodology and a robust analysis of the geometrical nonlinear effects on several detuned configurations of mistuned bladed-disk have been proposed in [START_REF] Picou | Robust dynamic analysis of detuned-mistuned rotating bladed disks with geometric nonlinearities[END_REF]. The results underline a complex nonlinear dynamical behavior and highlight a sensitivity of the nonlinear response to the detuning in presence of mistuning. Nevertheless, a too few number of detuned configurations is available and is not sufficient to perform a full optimization analysis with respect to the set of all possible detuned configurations of mistuned structures. This latter one has a huge dimension that increases exponentially with respect to the number of blades.

Novelties of the paper. The first novelty is the development of an efficient computational methodology in nonlinear stochastic dynamics for reducing the computational cost using a high-fidelity computational model (HFCM), and a deep computational analysis of a 12-bladed-disk structure that is representative of industrial turbomachines in order (i) to understand the role played by the geometrical nonlinearities on the dynamical behavior and (ii) to construct the solution of the detuning optimization problem and to analyze the optimization results. For the case for which the set of all the possible detuned configurations cannot be computed with the HFCM, we propose a novel probabilistic approach consisting in estimating a small subset (a few units) of detuned configurations that improve the dynamical responses with respect to the configuration without detuning in presence of random mistuning. This approach is based on a probabilistic formulation of the combinatorial optimization problem and on the use of a probabilistic learning tool for estimating the cost function without using the HFCM.

Organization of the paper. Section 2 is devoted to the formulation, the methodology, the computational model, and the algorithms for the nonlinear stochastic dynamics analysis of rotating detuned bladed-disks with geometrical nonlinearity in presence of random mistuning. The chosen formulation is compatible with the use of the nonparametric probabilistic approach of uncertainties for modeling the random character of the mistuning. We describe in details the specific numerical difficulties encountered (i) in the computation of the nonlinear internal forces and its related tangent stiffness matrix, and (ii) in the use of the nonlinear solver. In Section 3, we present a full analysis of a representative computational model of the 12-bladed-disk structure and the analysis of the detuned configurations. All the numerical parameters involved in the computational process are carefully optimized in order to be able to generate a full data basis for constructing the exact solution of the detuning optimization problem. In Section 4, for the general case for which the number of detuned configurations is very large, we present a probabilistic formulation of the combinatorial optimization problem and a probabilistic learning tool for constructing a surrogate model of the cost function. This formulation is able to capture the amplification level of a detuned configuration of the mistuned structure with respect to the configuration without detuning in presence of random mistuning (pure mistuned configuration). The statistics post-processing of the computational results yields a small subset of optimized configurations whose characteristics are analyzed in details. The proposed approach is validated using the computational analysis performed in Section 3 for the 12-bladed-disk structure.

Assumptions

In the present work, the detuned bladed-disk structure has N blades and is assumed (1) to be made up of a linear elastic material, [START_REF] Ewins | The effects of detuning upon the forced vibrations of bladed disks[END_REF] to be submitted to external forces whose amplitude is sufficiently large so that the structure undergoes geometrical nonlinear effects induced by large displacement, (3) to rotate at a constant velocity Ω around the rotation axis defined by (0, Z), [START_REF] Bladh | Component-mode-based reduced order modeling techniques for mistuned bladed disks. part 1: theoretical models[END_REF] to be constituted of different types of blades whose spatial distribution characterizes a detuned configuration. In this work, two different sector types are considered. These two types have an identical geometry but the blades have different material properties.

We have then to distinguish the tuned configuration from the detuned configurations. The tuned structure is related to the conceptual structure that has a perfect N-order cyclic symmetry. In consequence, the geometrical domain, the constitutive equations and the boundary conditions related to the generating sector are invariant under the 2π/N rotation around the symmetry axis (0, Z). A detuned configuration corresponds to a modification of the tuned structure for which there is a spatial distribution of different types of sectors. For instance, if N = 12 and if there are two different types of sectors A or B, a pattern denoted by 3A4B3A2B is constituted of 3 consecutive blades with type A, 4 with type B, 3 with type A and 2 with type B. The tuned configuration can also be viewed as a particular detuning case. It should be also noted that detuned configurations can also be sub-cyclic which means that they present a cyclic symmetry with order n < N, where n is a divisor of N. All these detuned configurations are characterized by the superscript .

The real systems that are related to these detuned structures are modeled by introducing mistuning (unintentional mistuning) that is modeled by using the nonparametric probabilistic approach of uncertainties [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF]. We then can distinguish the pure mistuning situation, from which the conceptual system is modeled with the tuned structure from the detuned structures that are considered in presence of mistuning.

Description of the geometric nonlinear boundary value problem of a detuned configuration

The superscript is omitted in this Section. Let Ω be the three-dimensional bounded domain of the physical space R 3 corresponding to the steady configuration of a detuned structure observed in the rotating frame that is defined as the reference configuration of the boundary value problem [START_REF] Desceliers | Non-linear viscoelastodynamic equations of three-dimensional rotating structures in finite displacement and finite element discretization[END_REF]. The boundary ∂Ω is such that ∂Ω = Γ 0 ∪ Γ E with Γ 0 ∩ Γ E = ∅ and the external unit normal to boundary ∂Ω is denoted by n = (n 1 , n 2 , n 3 ) (see Fig. 1). The boundary part Γ 0 corresponds to the fixed part (in the rotating frame) of the structure whereas the boundary part Γ E is submitted to an external surface force field. A total Lagrangian formulation is chosen. Consequently, the mechanical equations are written with respect to the reference configuration in the rotating frame. Let x = (x 1 , x 2 , x 3 ) be the position of a point belonging to domain Ω. The displacement field expressed with respect to the reference configuration is denoted as x(x, t) = x 1 (x, t), x 2 (x, t), x 3 (x, t) . It should be noted that the surface force field G(x, t)) = G 1 (x, t)), G 2 (x, t)), G 3 (x, t) acting on boundary Γ E and the body force field g(x, t)) = g 1 (x, t)), g 2 (x, t)), g 3 (x, t) acting on domain Ω correspond to the Lagrangian transport into the reference configuration of the physical surface force field and to the physical body force field applied on the deformed configuration. The nonlinear boundary value problem is written,

Γ 0 Γ E Ω n g G 3
ρ ∂ 2 x i ∂t 2 + 2ρ r ik ∂x k ∂t + ρ r ik r k (x + x ) = ∂ F i j S jk ∂x k + g i (t) in Ω , (1)
x i = 0 on Γ 0 , (2)
F i j S jk n k = G i (t) on Γ E , (3)
in which ρ > 0 is the mass density expressed in the reference configuration, where the (3 × 3) matrix [r(Ω)] whose component [r(Ω)] i j is denoted as r i j for the sake of clarity, is such that

[r(Ω)] i j = r i j = -Ω ε i j3 , (4)
where ε i jk is the Levi-Civita symbol such that ε i jk = ±1 for an even or odd permutation and ε i jk = 0 otherwise. In Eq. ( 17), {F i j } i j is the deformation gradient tensor whose components F i j are defined by

F i j = x i, j + δ i j , (5)
in which δ i j is the Kronecker symbol such that δ i j = 1 if i = j and δ i j = 0 otherwise. In Eq. ( 5), x i, j denotes the partial derivative ∂x i /∂x j . For a linear elastic material, the second Piola-Kirchoff symmetric stress tensor {S i j } i j is written as

S i j = σ (g) i j + a i jk E k , (6)
in which {σ (g) i j } i j is the symmetric Cauchy stress tensor acting on the reference configuration observed in the rotating frame. In Eq. ( 6), the fourth-order elasticity tensor {a i jk } i jk satisfies the usual symmetry and positive-definiteness properties. The Green strain tensor {E i j } i j is written as the sum of linear and nonlinear terms such that

E i j = ε i j + η i j , (7)
in which ε i j = (x i, j + x j,i )/2 and where η i j = x s,i x s, j /2. In the following, when no confusion is possible, symbol t denoting time is removed from the equations. Note that the first and the second partial derivatives of the displacement field with respect to t are denoted by ẋ and ẍ. In the rotating frame, we define the admissible space C ad of the displacements as the set of all the sufficiently differentiable functions x → δx(x) from Ω in R 3 such that δx = 0 on Γ 0 (in the rotating frame). For all fixed t, the weak formulation of the geometric nonlinear boundary value problem consists in finding the unknown displacement field x(t) in C ad such that, for all δx in C ad ,

m(ẍ, δx) + c(ẋ, δx) + d(ẋ, δx) + k (1) (x, δx) + k (2) (x, x, δx) + k (3) (x, x, x, δx) = f (δx) , (8)
in which

f (δx) = Ω δx i g i dx + Γ E δx i G i ds , (9)
with ds the surface element and where

m(x, δx) = Ω ρ x i δx i dx , m(x, x) > 0 , (10)
c(x, δx) = Ω 2 ρ r i j x j δx i dx , c(x, x) = 0 . ( 11
)
The bilinear form k (1) is written as

k (1) (x, δx) = k (e) (x, δx) + k (c) (x, δx) + k (g) (x, δx) ,
in which the positive-definite bilinear form k (e) , the negative-definite bilinear form k (c) , and the bilinear form k (g) are written as

k (e) (x, δx) = Ω a jk m ε m (x) ε jk (δx) dx , (12)
k (c) (x, δx) = Ω ρ r ik r k x δx i dx , (13)
k (g) (x, δx) = Ω σ (g) i j x s,i δx s, j dx , (14)
The multilinear forms k (2) and k (3) are such that

k (2) (x, x, δx) = Ω a jk m η m (x)ε jk (δx) dx + Ω a jk m x s, j δx s,k ε m (x)dx , (15)
k (3) (x, x, x, δx) = Ω a jk m x s, j δx s,k η m (x) dx. (16)
In Eq. ( 8), the damping bilinear form is modeled with a Rayleigh model such that

d(x, δx) = α m(x, δx) + β k (1) (x, δx) , (17)
in which α and β are two given positive or zero real numbers, not equal to zero simultaneously. The multilinear form k (2) (x, x, δx) can be written as

k (2) (x, x, δx) = 1 2 ( k (2) (x, x, δx) + k (2) (x, δx, x) + k (2) (δx, x, x)) , (18)
in which the multilinear form k (2) is defined by

k (2) (x, δx, z) = Ω a jk m x s, δx s,m z j,k dx , (19)
and satisfies the following symmetry property, k (2) (x, δx, z) = k (2) (δx, x, z) for all z in C ad . Finally, using the symmetry and the positive-definiteness property of elasticity tensor a allows for rewriting the multilinear form k (3) (x, δx, z, δz) as

k (3) (x, δx, z, δz) = 1 2 Ω a jk m x α, δx α,m z β, j δz β,k dx , (20)
which shows that k (3) (x, δx, x, δx) > 0.

Construction of the NL-ROM related to a detuned configuration

The computational model is constructed by discretizing Eq. ( 8) by the finite element method and is referred as the high-fidelity computational model (HFCM). Such a nonlinear computational model is not explicitly constructed because of the large number of degrees of freedom occurring in industrial finite element models. This Section is devoted to the construction of the nonlinear reduced-order model (NL-ROM) of a detuned configuration without mistuning of the rotating structure. The superscript characterizing the detuned configuration is omitted for clarity.

Let {ϕ 1 , . . . , ϕ M } be a given finite family of algebraically independent vectors in C ad . We introduce the notation A αβ = a(ϕ β , ϕ α ) in which a(x, δx) is one of the bilinear forms defined by Eqs. [START_REF] Madden | Reduced-order modeling approach for blisks with large mass, stiffness, and geometric mistuning[END_REF] to [START_REF] Batailly | Full three-dimensional rotor/stator interaction simulations in aircraft engines with time-dependent angular speed[END_REF]. Similarly, we introduce the notations K (2)

αβγ = k (2) (ϕ γ , ϕ β , ϕ α ), K (2) αβγ = k (2) (ϕ γ , ϕ β , ϕ α ), K (3) αβγδ = k (3) (ϕ δ , ϕ γ , ϕ β , ϕ α )
, and F α = f (ϕ α ). We define the vector u = (u 1 , . . . , u M ) of the generalized coordinates as a new set of unknowns obtained by projecting the reference nonlinear response x(x, t) on the vector space spanned by the finite family {ϕ 1 , . . . , ϕ M },

x(x, t) = M β = 1 ϕ β (x) u β (t) . (21)
Let δx(x) = ϕ α (x) be a test function in C ad . Substituting x(x, t) and δx(x) into Eq. ( 8), and using Eqs. ( 9) to (20) yields the following set of M nonlinear coupled differential equations such that, for α ∈ {1, . . . , M},

M αβ üβ + (C αβ + D αβ ) uβ + (K (e) αβ + K (c) αβ + K (g) αβ ) u β + F NL α (u) = F α , (22)
in which

F NL (u) = (F NL 1 (u), . . . , F NL M (u))
is the R M -vector of the reduced nonlinear internal forces whose component α is written as

F NL α (u) = K (2) αβγ u β u γ + K (3) αβγδ u β u γ u δ , (23)
in which

K (2) αβγ = 1 2 K (2) αβγ + K (2) βγα + K (2) γαβ , (24)
and where F = (F 1 , . . . , F M ) is the R M -vector of the reduced external forces whose component α is F α = f (ϕ α ). The computational procedure concerning the computation of F NL (u) requires to use the finite element method in order (1) to construct the quadratic stiffness contribution K (2) αβγ , to retrieve the quadratic stiffness contribution K (2) αβγ according to Eq.( 24), to construct the cubic stiffness contribution K (3) αβγδ and to compute the vector of the nonlinear reduced internal forces F NL (u). The choice of such computational procedure is required by the choice of modeling the mistuning phenomenon through the nonparametric probabilistic approach [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF]. All the details concerning the computation of the nonlinear internal reduced force vector can be found in [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF].

The projection basis can be constructed as the linear elastic modes of the rotating structure associated with its first eigenfrequencies ordered in ascending order. It should be noted that in the case of a rotating structure, the generalized eigenvalue problem allowing such a computation is a nonsymmetric eigenvalue problem yielding a complex projection basis, due to the presence of the gyroscopic coupling term. Another projection basis can be chosen as the family of linear elastic modes constructed without considering the gyroscopic coupling term. With such a choice, the speed of convergence with respect to M of the projection remains fast enough for bladed-disks whose blades have a reasonable slenderness. When the gyroscopic coupling term is kept, an alternative method avoiding the calculation of complex modes can be used as proposed in [START_REF] Picou | Robust dynamic analysis of detuned-mistuned rotating bladed disks with geometric nonlinearities[END_REF]. It consists in using the proper orthogonal decomposition method (POD) applied to the nonlinear dynamical response of the rotating structure (including all the coupling terms) issued from a first nonlinear reduced-order model obtained with the linear elastic modes calculated without the gyroscopic coupling term.

Let n be the number of dofs of the HFCM. Let [Φ] be the (n × M) matrix whose columns are the vector of the projection basis (in the finite element discretization). At time t, the physical response of the detuned bladed-disk, constructed using the HFCM, is written as

X(t) = [Φ] U(t) , (25)
in which X(t) is the R n -vector of the physical displacements and U(t) is the R M vectors of the generalized coordinates U 1 (t), . . . , U M (t).

2.4. Construction of the NL-SROM related to a detuned configuration in presence of mistuning In this paragraph, we are interested in constructing the nonlinear stochastic reduced-order model (NL-SROM) related to a detuned configuration of the bladed-disk in presence of random mistuning. It is assumed that the mistuning only affects the linear elastic internal force and the nonlinear forces, which are therefore uncertain. Such mistuning is modeled through the nonparametric probabilistic approach [START_REF] Soize | Uncertainty Quantification[END_REF] for geometric nonlinearities [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF], which is briefly summarized below. A global stiffness (P×P) real symmetric matrix [K NL ], with P = M(M + 1), containing all the linear, quadratic, and cubic elastic stiffness contributions, is constructed such that

[K NL ] =       [K (e) ] [ K (2) ] [ K (2) ] T 2 [K (3) ]       , (26)
in which [K (e) ] is the (M × M) matrix defined by [K (e) ] αβ = K (e) αβ and where the block matrices [ K (2) ] and [K (3) ] are the (M × M 2 ) and the (M 2 × M 2 ) real matrices resulting from the following reshaping operation, [ K (2) ]

αJ = K (2) αβγ , J = (β -1)M + γ , (27)
[K (3) ]

I J = K (3) αβγδ , I = (α -1)M + β , J = (γ -1)M + δ . (28)
In [START_REF] Mignolet | Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems[END_REF], it is proven that [K NL ] is a positive-definite matrix that is modeled by a random matrix [K NL ]. In this work, we use the nonparametric probabilistic modeling proposed in [START_REF] Capiez-Lernout | An improvement of the uncertainty quantification in computational structural dynamics with nonlinear geometrical effects[END_REF] for random matrix [K NL ],

[K NL ] = [L K ] T [G(δ K )] [L K ] + [∆K]with [∆K] = [K NL ] -[L K ] T [L K ] . (29)
The matrix [L K ] is a full (Q × P) matrix for which Q

P, defined by [L K ] = [Λ Q ] 1 2 [Ψ Q ] T .
The columns of [Ψ Q ] are the Q eigenvectors of matrix [K NL ] related to its first Q largest eigenvalues that are stored in the diagonal matrix

[Λ Q ]. The full (Q × Q) random matrix [G(δ K )]
is then constructed according to the nonparametric probabilistic approach [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF] for which the uncertainty level is controlled by a scalar hyperparameter δ K . The random matrix [K NL ] is blockdecomposed similarly to Eq. ( 26), which allows for extracting the random stiffness terms K (e) αβ , K

αβγ , and K (3) αβγδ . The random quadratic term K (2) αβγ is then reconstructed according to Eq. ( 18) by using the formula

K (2) αβγ = 1 2 K (2) αβγ + K (2) βγα + K (2)
γαβ .

The NL-SROM of a detuned configuration of the rotating bladed-disk in presence of mistuning is then written (see Eq. ( 25)) as

X (t) = [Φ ] U (t) , (31)
in which X (t) is the R n -valued random vector corresponding to the nonlinear random displacements at time t of the HFCM for the detuned configuration in presence of mistuning. At time t, the R n -valued random vector U (t) of the generalized coordinates is then solution of the random nonlinear coupled differential equations, which are written, for α ∈ {1, . . . , M}, as

M αβ Ü β + (C αβ + D αβ ) U β + (K (e), αβ + K (c), αβ + K (g), αβ ) U β + F NL, α (U ) = F α , (32)
in which the component α of the R M -valued random vector F NL, (U) is

F NL, α (U ) = K (2), αβγ U β U γ + K (3), αβγδ U β U γ U δ . (33)
The observation (quantity of interest) that is adapted to the detuning context in presence of random mistuning and that allows to quantify the dynamical amplification of the detuned configuration with random mistuning with respect to the configuration without detuning in presence of random mistuning, is defined in Section 3 (see Eq. ( 38)).

Numerical aspects related to the algorithms

This section is devoted to the computational strategy used for solving the NL-SROM. In the detuning context, the number n c of distinct detuned configurations rapidly increases with the number N of blades of the bladed-disk as n c

2 N /N. Furthermore, the NL-SROM is solved with the Monte-Carlo numerical simulation using n sim realizations for which the numerical convergence has to be achieved in order to perform statistics post-processing. This means that there is a large number n c × n sim of nonlinear differential equations with order M to be solved. Note that, in general for industrial bladed-disks, n c is very large and it is impossible to solve n c ×n sim times the nonlinear reduced-order computational model, even with an efficient algorithm. Only a small number N d n c of detuned configurations in presence of mistuning can be computed. This is the reason why a surrogate model will be constructed in Section 4 with a probabilistic learning tool in order to solve the probabilistic-based combinatorial optimization problem and which will allow an optimal detuned configuration to be identified. In consequence, the computational algorithms require to be efficient in terms of computational costs. A careful attention has been paid in order to be able to browse the set of all the possible detuned configurations in presence of mistuning using the HFCM, when n c is not to large. This approach is presented and analyzed in Section 3 and the obtained results will be used for validating the solving method of the probabilistic-based combinatorial optimization problem using a probabilistic learning tool, presented in Section 4.

Computational strategy for dealing with geometrical nonlinearities

In a deterministic finite element context, it should be noted that the construction of the linear reduced operators in Eq. ( 22) is usual and can easily be extracted using commercial softwares. From a computational point of view, the vector of the nonlinear reduced internal forces F NL (U) and the tangential matrix, for which the nonlinear stiffness contribution is written as

[K NL T (U)] αβ δU β = ∂F NL α (U) ∂U β δU β ,
can also be extracted. Nevertheless, vector F NL (U) is usually constructed by directly computing the contribution related to each finite element, without distinguishing the quadratic stiffness entries from the cubic stiffness entries. In such a case, the assemblage is quickly performed as a simple summation. In general, tangential matrix [K NL T (U)] is deduced from finite differences techniques. The present methodology requires the knowledge of entries K (2), αβγ and K (3), αβγδ . Its detailed construction can be found in [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF]. In the detuning context in presence of random mistuning, the computational procedure allowing for computing the random quantities F NL (U) and [K NL T (U)] issued from the NL-SROM is summarized hereinafter. 1. Commercial software extraction or direct construction of the reduced linear elastic stiffness K (e) αβ = k (e) (ϕ β , ϕ α ) according to Eq. ( 12). 2. Construction of the quadratic stiffness contribution K (2) αβγ = k (2) (ϕ γ , ϕ β , ϕ α ) and of the cubic stiffness contribution K (3) αβγδ = k (3) (ϕ δ , ϕ γ , ϕ β , ϕ α ) according to Eqs. [START_REF] Mehrdad Pourkiaee | A reduced order model for nonlinear dynamics of mistuned bladed disks with shroud friction contacts[END_REF]

αβγ , and K (3) αβγδ and retrieving random stiffness K (2) αβγ according to Eq. (30) 6. Computational evaluation of random nonlinear internal reduced force vector F NL (U) according to Eq. (33) 7. Depending on the computational algorithms used for solving nonlinear equations, computing the tangential matrix related to the nonlinear contribution

[K NL T (U)] such that [K NL T (U)] αβ = K (2) αβγ U γ + K (2) αγβ U γ + K (3) αβγδ U γ U δ +K (3) αδβγ U γ U δ + K (3) αδγβ U γ U δ .

Computational strategy for the nonlinear solver

The problem investigated concerns the forced regime of the bladed-disk, that is to say the nonlinear forced response that does not depend on initial conditions. The external excitation is thus a square integrable real-valued function on R whose Fourier transform has a compact support (see Section 3.3). For a given detuned configuration with random mistuning, the NL-SROM is solved in the time domain. An implicit unconditionally stable integration scheme corresponding to the Newmark method with an averaging acceleration scheme is chosen. With such integration scheme, for a given sampling time t i , we then have to solve the following set of nonlinear coupled equations such that, for α = 1, . . . , M,

[K eff ] αβ U β (t i ) + F NL α (U(t i )) = F eff α (t i ) . (34)
In Eq. ( 34), the random matrix [K eff ] is the reduced effective stiffness matrix that is independent from sampling time. Note that [K eff ] is positive definite almost surely for a fixed structure and is without signature for a rotating structure due to the antisymmetric gyroscopic coupling term. The random vector F eff (t i ) models the effective external load that has to be updated at each sampling time t i . There are several ways to solve iteratively such nonlinear coupled equation. The simplest method that is used is the fixed point method. The main advantage is that this iterative scheme is computationally expedient because it does not require the evaluation of the tangential matrix. Nevertheless, the algorithm does not necessarily converge. When the fixed-point algorithm does not converge, it can be replaced by a Newton-Raphson algorithm [START_REF] Borst | Non-Linear Finite Element Analysis of Solids and Structures[END_REF]. It should be noted that such algorithm requires the computation of the tangential matrix that is particularly time consuming. Furthermore, for nonlinear computations corresponding to a snap-back or snap-through behavior, the Newton algorithm does not converge beyond the critical limit corresponding to a null tangential stiffness matrix. In this case, the arc-length algorithm appears to be efficient for computing the solution. Nevertheless, such algorithm introduces a new additional scalar parameter ζ(t i ) that multiplies the right-hand side member of the nonlinear equation. In this case, we then introduce the residual random vector R(U

(t i ), ζ(t i )) whose component α is defined by R α (U(t i ), ζ(t i )) = ζ(t i )F eff α (t i ) -[K eff ] αβ U β (t i ) -F NL α (U(t i )) . ( 35
)
The algorithm is iterative. Several numerical approximations can be made. The first one consists in keeping the same tangential matrix [K NL T (U(t i ))] along the successive iterations. The second one concerns the computation of the increment related to parameter ζ(t i ). The pure spherical arc-length method is based on the fact that there always exists a hypersphere that intersects the load curve [START_REF] Crisfield | A fast incremental/iterative solution procedure that handles "snap-through[END_REF]. It then requires to solve a quadratic equation and to select amongst the two solutions the physical root. Improvements have been proposed based on linearization methods, for which the hypersphere is replaced by an hyperplan that is updated at each iteration [START_REF] Ramm | Strategies for tracing the nonlinear response near limit points[END_REF] or not updated [START_REF] Riks | The application of Newton's method to the problem of elastic stability[END_REF]. For a given sampling time t i , such iterative scheme allows for computing a solution (U(t i ), ζ(t i )) but does not allow for controlling the value of ζ(t i ). In the present case, the process is incrementally repeated in order to increase the value of ζ(t i ) that has to be controlled to 1. In order to achieve that, an adaptive arc length, depending on the number of iterations necessary to obtain the convergence of the previous increment is also implemented according to [START_REF] Borst | Non-Linear Finite Element Analysis of Solids and Structures[END_REF]. An unusual procedure is then added to the algorithm concerning parameter ζ(t i ). The state of the algorithm corresponding to the previous increment is stored. As the new incremental value of ζ(t i ) is found to be upper than 1, the algorithm is rewind to the previous increment and the computation is set again with the half of the arc-length. Such additional procedure is repeated until parameter ζ(t i ) reaches 1 within a numerical tolerance set.

Computational analysis of the role played by the geometrical nonlinearities in the detuned bladed-disk in presence of random mistuning and detuning optimization

For the computational analysis, we have deliberately chosen a blisk (turbomachine component comprising both rotor disk and blades) with mistuning with a relatively small number of blades so that the number of all the possible detuned configurations be also small. This will allow us

(1) to present a complete analysis of the role played by geometrical nonlinearities in the detuned blisk in presence of random mistuning,

(2) to build an exact solution (as the reference) of the combinatorial optimization problem introduced in Section 4. This reference will allow for validating the formulation of the combinatorial optimization problem (see Sections 4.1 to 4.5) and for validating the probabilistic machine learning tool for constructing a metamodel of the cost function (see Section 4.6 to 4.8). If we had chosen a larger number of blades, we would not have been able to build the reference solution of the combinatorial optimization problem, reference that is required for validation.

Description of the finite element model

The structure under consideration is a blisk with order N = 12 and with n w = N = 12 blades whose finite element computational model is constructed with hexahedral solid finite elements with 8 nodes. The finite element mesh is issued from [START_REF] Bladh | Component-mode-based reduced order modeling techniques for mistuned bladed disks. part 2: application[END_REF] for which the order of the cyclic symmetry has been modified from 24 to 12 as also used in [START_REF] Han | Optimization of intentional mistuning patterns for the mitigation of the effects of random mistuning[END_REF]. The main motivation of such reduction of the number of blades is to decrease the number of possible detuning patterns in order to have a reasonable number of detuning possibilities, as previously explained. In the present case, the aim is to constitute a full reference data basis in the nonlinear mistuned/detuned context. When the detuning concept is approached with 2 possible types of blades, there are 352 possible detuned configurations when N = 12 that reduces the number of possibilities by a factor of almost 2 N-1 2 000 with respect to a structure with cyclic symmetry of order 2N = 24. Let (0, X, Y, Z) be the Cartesian reference coordinate system for which (0, Z) coincides with the rotational axis for the blisk. Fig. 2 shows the finite element mesh of the blisk whose computational characteristics are summarized in Table 1. There are n = 27 108 degrees of freedom (dof). The blisk is made up of a homogeneous and isotropic material with mass density 7 860 Kg/m 3 ,

Modal characteristics of the tuned blisk

The eigenfrequencies (natural frequencies) of the structure are computed using the cyclic symmetry. Fig. 4 displays the eigenfrequencies of the tuned blisk with respect to its circumferential wave number. The graph exhibits several veerings. Let ν h,i be the eigenfrequency number i related to the number h of nodal diameters. In the present case, the veering corresponding to h = 2 nodal diameters corresponds to (double) eigenfrequency ν 2,4 = 4 845.18 Hz that is related to a dominant blade motion and to (double) eigenfrequency ν 2,5 = 5 091.09 Hz that is related to a blisk global motion as can be shown in Fig. 5.

Choice of the external excitation and time-frequency sampling

Since we are interested in the forced response of the structure, the presence of geometrical nonlinearities requires to solve the nonlinear equations in the time domain. The time-dependent external excitation is modeled by the R n -vector F(t) whose block decomposition according to each sector is given by

F(t) = F 0 (t), . . . , F N-1 (t) such that F j (t) = s 0 g j (t) F , j ∈ {0, 1, . . . , N -1} ,
in which s 0 is the load intensity that allows the nonlinearity rate to be calibrated, t → g j (t) is a square integrable real-valued function on R, which characterizes the time evolution of the load, and where F is a R n -vector describing the space localization of the load related to a given sector, in which integer n is equal to n/N. The Fourier transform of function t → g j (t) is the C-valued function again noted g j , but for which its argument is 2πν, such that

g j (2πν) = R e -2ιπνt g j (t) dt , ν ∈ R , in which ι is the imaginary unit. Complex-valued function ν → g j (2πν) is assumed to have a compact support B e ∪ B e ∈ R, independent of j, with B e = [-ν max , -ν min ] and B e = [ν min , ν max ].
This support is chosen as the frequency band of excitation for which ν min = 4 700 Hz and ν max = 5 200 Hz (see Fig. 7) in order that B e contains the veering region for h = 2 (see Fig. 4). Time function g j differs from one sector to another one by a constant phase shift ϕ j = (4π j)/N. This choice of the phase allows the modes with h = 2 nodal diameters to be excited for the linearized dynamical system. For all t in R, we write

g j (t) = S k=1 2 π t sin(π ∆ν t) cos(2πs k ∆ν t + ϕ j )
in which ∆ν = (ν max -ν min )/S and where

s k = ν min /∆ν + (2k -1)/2. Consequently, |g j (2πν)| = 1 if ν belongs to B e ∪ B e ∈ R
and is equal to 0 if does not. Therefore, the structure is uniformly excited in band B e . From a numerical point of view, the computation is carried out on a truncated time domain T = [t ini , t ini + T ]. The initial time is chosen as t ini = -0.06 s yielding a null initial load. The time duration T is then adjusted so that the system be returned at its equilibrium state within a given numerical tolerance for both linear and nonlinear computations. Even the fundamental eigenfrequency does not belong to excitation frequency band B e , it can indirectly be excited through the geometrical nonlinear effects. Let ξ 0 = 0.00754 be the critical damping rate at fundamental frequency ν 0 as shown in Fig. 3. Time duration is chosen as T = 0.128 s, which ensures the system to return to its equilibrium state with a negligible relative tolerance τ = e -2πξ 0 ν 0 T < 10 -12 when fundamental frequency is excited. The sample frequency ν e and the number of time steps are then chosen as ν e = 64 000 Hz and n t = 8 192 yielding a constant sampling time step δt = 15.6 × 10 -6 s and a constant sampling frequency step δν = 7.8125 Hz. Function ν → g j (2πν) is estimated by using FFT. A convergence analysis is carried out with respect to the number S + 1 of subdivisions of frequency band B e . Let Conv 1 (S ) be defined by

Conv 1 (S ) = R + |g j (2πν)| dν ν max -ν min .
Fig. 6 shows the graph of S → Conv 1 (S ). It can be seen that S = 30 yields a good representation of |g j (2πν)|. Fig. 7 displays the graphs of t → g 0 (t) and ν → |g 0 (2πν)|. The spatial repartition of the external load described by the normalized vector F is such that an external point load is applied along all directions at the excitation nodes located at the tip of each blade with a constant phase shift π/3. The load intensity s 0 is adjusted to get sufficient geometrical nonlinear effects that modify the vibrational behavior of the blisk.

Convergence analysis with respect to time sampling

The efficiency and the accuracy of the Newmark time integration scheme coupled with a Newton-Raphson procedure is investigated. The load intensity is set to s 0 = 10 N that corresponds to a situation for which nonlinear effects occur, involving displacements of magnitude 2 × 10 -3 m at the tip of blade, knowing that the blade thickness is 6.7 × 10 -3 m. The nonlinear reduced-order model is then constructed with M = 72 elastic modes corresponding to increasing eigenfrequencies up to ν 6,12 = 10 835.13 Hz. The frequency band of analysis B is chosen as B = [300 , 6 000] Hz. We are then interested in the nonlinear computations for which the nonlinear response is calculated in the truncated time domain T = [-0.06, 0.068] s. For performing the convergence analysis, the values of parameters δt, n t , and ν e are given in Table 2. T (s) 0.128×10 -3 0.128×10 -3 0.128×10 -3 0.128×10 -3 0.128×10 -3 0.128×10 First of all, it can be seen that the structure returns to its equilibrium state at the end of the simulation. The time response is strongly sensitive to the choice of δt. The response in the frequency domain displays consequent nonlinear effects, since the structure responds outside the frequency band of excitation B e . The convergence with respect to the choice of δt is relatively quick when analyzing the structure in the excitation frequency band. The fundamental frequency of the blisk is also excited through the nonlinearities and it can be seen that this response is strongly over-estimated when choosing δt = 62.5×10 -6 s and becomes stable from for time steps lower than δt = 15.6×10 -6 s, which is confirmed by the zoom of the nonlinear response in the low-frequency range displayed in Fig. 10. From now on all the computations in the time domain are done taking δt = 15.6×10 -6 s.

Convergence analysis with respect to the order of the reduced-order model

The accuracy of the NL-ROM is investigated with respect to the order M of the reduced-order model. The elastic modes are sorted by increasing order of their eigenfrequencies belonging to a family of modes that are represented in Fig. 4. The load intensity is set to s 0 = 10 N and the time step is set to δt = 15.6×10 -6 s. Denoting always u = (u 1 , . . . , u M ) as the vector of the M modal coordinates, the convergence function

M → Conv 3 (M) is defined by Conv 3 (M) = B u(2πν; M) 2 dν , in which u(2πν; M) 2 = M i=1 |u i (2πν; M)| 2 and
where the C M -vector u(2πν; M) is the Fourier transform of the R M -valued function t → u(t; M). The superscript lin is added for the computation corresponding to the response of the linearized system. Fig. 11 displays the graphs of functions M → Conv lin 3 (M) (upper graph) and M → Conv 3 (M) (lower graph). It should be noted that for each considered order M, the corresponding value of function Conv 3 (M) has been obtained by performing a whole deterministic nonlinear analysis using the NL-ROM parameterized by M. As expected, it can be seen that the converged values of both convergence functions differ, that underlines the presence of effects induced by the geometrical nonlinearities. As expected, it can also be seen that the convergence is slower in the nonlinear case. It should also be noted that the nonlinear tuned case yields response levels that are much more lower than the corresponding linearized case, which means that the geometrical nonlinearities have a stiffening effect on the vibrational behavior of the bladed-disk. Note also that the convergence analysis exhibits a maxi- mum when the nonlinear reduced-order model of order M = 48 is used for the computations. The main contribution comes from mode with ν 2,5 = 5091.09 Hz, which is the mode with the higher frequency that is located in the excitation frequency band and that is directly excited by the load vector. Such an overestimation is then corrected by including additional modes in the nonlinear reduced-order model. From now on, the order of the NL-ROM is fixed to M = 72. An analysis of the most contributing generalized coordinates is carried out over B and presented in Table 3. In the linear case, it is seen that the most contributing generalized coordinates are related to natural eigenfrequencies ν 2,5 = 5091.1 Hz and ν 2,4 = 4845.2 Hz. Moreover there are also two other contributions that are non negligible, that are located in B e , and that are not related to circumferential wave number h = 2. This is not theoretically expected since the load has been constructed to only excite circumferential wave number h = 2. Nevertheless, in the present case, the forced response is solved in the time domain. The signal is theoretically defined for t ∈ R but is numerically truncated, which explains those secondary contributions. In the nonlinear case, all the natural frequencies located in B e are contributing to the nonlinear response. It should also be noted that the ratios of the generalized coordinates related to natural frequencies ν 2,2 = 2571.5 Hz and ν 0,5 = 5266.1 Hz, that are located outside the excitation frequency band do have a non negligible participation.

linear h 2

Effects of the geometrical nonlinearities

In order to have a better understanding of the structural behavior of the blisk, it is essential to perform a sensitivity analysis of the nonlinear response with respect to load intensity s 0 . The computations are carried out with numerical parameters set to δt = 15.6 × 10 -6 s and M = 72. Parameter s 0 controls the level of nonlinear effects that can be quantified as follows. A criterion is that the nonlinear response does no longer coincide with the response of the linearized system (small deformations, small displacements) as soon as vibrations appear outside the excitation frequency band B e . As previoulsly, let u = (u 1 , . . . , u M ) be the vector of the M modal coordinates and let s 0 → Conv 4 (s 0 ) be the function defined by

Conv 4 (s 0 ) = B e u(2πν; s 0 ) 2 dν B u(2πν; s 0 ) 2 dν , in which u(2πν; s 0 ) 2 = M i=1 |u i (2πν; s 0 )| 2 and
where the C M -vector u(2πν; s 0 ) is the Fourier transform of the R M -valued function t → u(t; s 0 ). Such function is equal to unity as long as the linear assumption remains valid. Consequently, there exists a threshold s min above which geometrical nonlinearities can no longer be neglected. Fig. 12 shows the graph of function s 0 → Conv 4 (s 0 ). The threshold value is found to be s min = 1 N.

3.7.

Convergence analysis with respect to the size Q of the random germ.

As defined in Section 2.4, the global stiffness matrix [K NL ] corresponding to a reduced-order model of order M = 72 has a size (P × P) with P = 72 × 73 = 5256. Figure 13 represents the function j → [Λ Q ] j j where the eigenvalues of matrix [K NL ] are sorted by decreasing values.

One chooses Q so that eigenvalues that are greater than 10 -6 the fundamental one and their corresponding eigenvectors are kept in the reduction so that Q = 440. Note that this value defined the length of the random vector used for constructing random matrix [G(δ K )]. Let be Conv 5 (Q) be the function defined by

Conv 5 (Q) = ||∆K(Q)|| 2 ||K NL || 2
It can be shown that Conv 5 (Q = 440) = 1.98 × 10 -6 which allows a good representation of the solution to be obtained and which is such that Q = 440 P = 5256. 3.8. Convergence analysis with respect to the number of Monte Carlo numerical simulations for the mistuned blisk without intentional detuning The load intensity is set to s 0 = 10 N that corresponds to significant nonlinear effects and the numerical parameters are set to M = 72 and δt = 15.6×10 -6 s. The pure mistuning (that is to say without intentional detuning) of the blisk is considered through the nonparametric probabilistic approach with a hyperparameter δ K set to δ K = 0.1.

Since the nonparametric probabilistic approach is used for modeling the mistuning, it seems important to precise what is the physical meaning of the choice of hyperparameter δ K = 0.1 in the pure mistuning case. Figure 14 displays the probability density function of the first three eigenfrequencies with circumferential wave number h = 0 using kernel density estimation in order to provide the level of variations of the natural frequencies induced by the pure mistuning when using the nonparametric probabilistic approach. The nonlinear random responses obtained in the pure mistuning case allows for the order of magnitude of the amplification to be estimated, yielding a mean level of amplification of 1.222, a corresponding standard deviation of 0.164 and a 98th quantile of 1.629.

The convergence analysis with respect to the number n sim of Monte Carlo numerical simulations is investigated. Let U = (U 1 , . . . , U M ) be the random vector of the M modal coordinates on the probability space (Θ, T , P). Let Conv 6 (n sim ) be the function defined by 3.9. Quantification of the computational costs with regard to the detuning context of the mistuned blisk These numerical computations are essential (1) to better understand the dynamic behavior of the blisk and (2) to optimize for the best the computational costs while keeping a predictive computational model that ensures a good accuracy on the nonlinear dynamical response. From now on, the load intensity is set to s 0 = 10 N that corresponds to significant nonlinear effects (see Fig. 12). The mistuning level is characterized by a dispersion level δ K = 0.1 and the numerical parameters are set to M = 72, δt = 15.6 × 10 -6 s, n sim = 500, and Q = 440. It is essential to underline that all these convergence analyses have to be carried out meticulously. The aim of this work concerns the optimization of the detuning of the blisk in presence of mistuning. The computational model of the detuned blisk is constructed from the knowledge of two compatible meshes of two different sector types denoted as 1 and 2. The reference sector 2 is obtained from nb 1 [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF] Let us denote by nb 1 and nb 2 the number of blades of type 1 and type 2 in a detuned configuration. Let n c (nb 1 , nb 2 ) be the number of possible detuned configurations having nb 1 and nb 2 blades of type 1 and 2 respectively. Let n c be the total number of detuned configurations for all the possible values of nb 1 and nb 2 . For the tuned blisk with cyclic symmetry of order N = 12, Table 4 gives the number n c (nb 1 , nb 2 ) of detuned configurations as a function of the number nb 1 and nb 2 of blades of type 1 and 2. There are a total number of n c = 352 possible detuned configurations (patterns).

The objective is to find the optimal configuration (pattern) that will reduce for the best the response amplification induced by mistuning. Another objective is to have a complete knowledge of the nonlinear dynamical behavior of all possible detuned configurations in order to constitute a reference data basis. In the present numerical study, all the computations are distributed according to the number n sim of Monte Carlo realizations and are made on workstations with 1 536 GB RAM and 30 cores (Intel(R) Xeon(R) Platinum 8 168 CPU@2.9Ghz). The generated elapsed time is about 6 hours to perform 1 computation corresponding to the nonlinear analysis of 1 detuned pattern in presence of mistuning, which means that the construction of the full data basis requires an elapsed time of about 88 days of computations.

Choice of the observation

The pure mistuning case corresponds to the usual situation for which the blisk is described with a tuned configuration in presence of mistuning. In that particular case, the deterministic case corresponds to a structure with a perfect cyclic symmetry with order N. The usual mistuning analyses, whether it is in a linear or nonlinear context, characterize the mistuning effects by introducing a random amplification factor that is expressed with respect to the deterministic resonance of the structure with perfect cyclic symmetry.

The detuned configuration in presence of mistuning is defined by a given distribution of the blades of type 1 and 2. The detuned configuration number is then represented by a vector w c, ∈ {0, 1} n w whose component w c, j is equal to 0 or 1 whether the blade is of type 1 or of type 2. Let us recall that for a given mistuning level, the drastic mistuning effects yield a response amplification with respect to its deterministic counterpart. The objective is to find detuned configurations whose mistuning effects induce less amplification than the one obtained with the tuned configuration in presence of mistuning. The idea is then to define a scalar quantity that is able to quantify the mistuning effects of the detuned configuration with respect to the unavoidable mistuning effects of the tuned configuration.

Let X ( j, 2πν i , θ k ) be the C 3n b -vector of the realization θ k of the 3n b -displacement dofs of blade number j ∈ {1, . . . , n w } in the frequency domain, taken at frequency ν i , and corresponding to the detuned configuration number . Let [t j ] be the (3n b × n) real matrix that allows for restraining the complex vector of the n displacement-dofs of the structure to those related to blade number j, we have

X ( j, 2πν i , θ k ) = [t j ] [Φ ] U (2πν i , θ k ) ,
in which [Φ ] ∈ M n,M is the matrix of the M elastic modes for configuration and where U (2πν i , θ k ) is the realization θ k ∈ Θ of the C M -valued random variable U (2πν i ) of the modal coordinates at frequency ν i .

For blade number j of configuration , at frequency ν i , and for realization θ k , we define the scalar observation Y ( j, 2πν i , θ k ) = X ( j, 2πν i , θ k ) . Since the mistuning context requires to consider the resonance of the most responding blade, we define the R + -valued random variable A whose realization a ,k = A (θ k ) is such that

a ,k = max j {max i Y ( j, 2πν i , θ k )} . ( 36
)
We define the associated random variable J max whose realization j ,k max = J max (θ k ) is such that

j ,k max = arg max j {max i Y ( j, 2πν i , θ k )} . (37)
In order to get a robust scalar quantity for characterizing the random nonlinear dynamical behavior of the blisk, an estimate of the maximum extreme value statistics of random variable A is introduced. It should be noted that, for a fixed -configuration, the maximum is used so that all realizations are in the confidence domain with the greatest probability. This estimation is constructed as follows. The available number of Monte Carlo numerical simulations is written as n sim = ν r ν e (for n sim = 500, ν e = 10 and ν r = 50). We then define, for r ∈ {1, . . . , ν r }, a ,r M = max k∈{ν e (r-1)+1,...,rν e } a ,k , in which a ,k is defined by Eq. ( 36) and we define its mean value, a M , such that

a M = 1 ν r ν r r=1 a ,r M .
The observation of the detuned -configuration with mistuning is then defined as the amplification factor q c, with respect to its tuned counterpart with pure mistuning,

q c, = a M a t M , (38)
in which superscript t is related to the tuned configuration.

3.11. Some remarks concerning the nonlinear detuning optimization and the linearized detuning optimization In this paragraph, we investigate the deterministic response of the two tuned configurations with cyclic symmetry with order N = 12. Concerning the tuned configuration with blades of type 1, the natural eigenfrequencies related to the veering with circumferential wave number h = 2 are the double eigenfrequencies ν Hz, it can be seen that the tuned configuration with 12 blades of type 2 only presents one resonance in B e . In the linearized case, this means that there will be some detuned configurations for which there will be only one resonance in frequency band of analysis. B. Let y type α (2πν i ), with α = {1, 2} be the deterministic scalar observation defined similarly to quantity Y ( j, 2πν i , θ k ). Figure 16 displays the graphs of a zoom of ν → y type α,lin (2πν i ) and ν → y type α (2πν i ). It is clearly seen that the nonlinear and linearized response yield very different dynamical behavior due to the external load level for which geometrical nonlinearities occur. Concerning the nonlinear response, the two tuned configurations do exhibit a complex vibrational behavior, with several resonances in and outside the excitation frequency band, which is of interest for considering the nonlinear detuned optimization. Concerning the linearized response, it is seen that the second resonance for the tuned configuration with blades with type 1 has a response level about 5 times smaller than the main resonance. It is also seen that the first resonance for the tuned configuration with blades with type 2 is not captured because the corresponding natural eigenfrequency is localized below the excitation frequency band. For that reason and for the considered case, it does not seems appropriate to consider linearized detuning optimization and to compare it with the nonlinear detuned optimization. 16: Graph of a zoom of the deterministic response for the nonlinear and linearized response of the tuned configurations with blades of type 1 and of type 2. Graphs of ν → y type 1 (2πν i ) (thick line), ν → y type 1,lin (2πν i ) (thin line), ν → y type 2 (2πν i ) (thick dashed line), ν → y type 2,lin (2πν i ) (thin dashed line). Excitation frequency band B e is represented by the light pink area.

3.12. Detuned configurations that decrease the amplification of the nonlinear response induced by the mistuning A detuned configuration , which yields a nonlinear mistuned response level that is smaller than the one obtained with the tuned configuration in presence of mistuning, is characterized by q c, < 1. Such configuration belongs to the set of optimal detuned configurations. In the detuning process, the tuned configuration w c,t = [0 0 0 0 0 0 0 0 0 0 0 0] in presence of mistuning is considered and yields amplification q c,t = 1.0000. The detuning is defined around this tuned configuration. This subset corresponds to n c = 216 possible detuned configurations for which the number nb 2 of blades of type 2 is less than or equal to n w /2 = 6 (see Table 4). However, in this paragraph, we will also consider as possible configurations, the detuned configurations around the other tuned configuration w c,t 2 = [1 1 1 1 1 1 1 1 1 1 1 1] corresponding to nb 2 > 6 and for which the amplification factor is q c,t 2 = 1.1252. In such a case (see Table 4), there are 136 additional detuned configurations and consequently, the total number of configurations is n c = 352, which are randomly ordered. The tuned configuration w c,t corresponds to = 49. Fig. 17 shows the graph of j → q c, j where q c, j are sorted by increasing order for n c = 352 (upper graph) and for n c = 216 (lower graph). It can be shown that the amplification factor increases from q c, 1 = q c,104 = 0.9476 until q c, 352 = q c,9 = 1.1774 corresponding to detuned configurations defined by w c,104 = [0 0 0 1 0 0 1 0 0 1 0 1] and w c,9 = [0 1 0 0 0 1 1 1 0 0 1 1]. By comparing these two graphs, it is shown that only one improving detuned configuration has a number nb 2 of blades with type 2 greater than 6. Furthermore, it is seen that the worst detuned configuration = 9 has a number nb 2 of blades with type 2 equal to 6. It can also be seen that most of the detuned configurations having a majority of blades with type 2, have a larger amplification factor.

From now on, only detuned configurations corresponding to n c = 216 possible detuned configurations for which will be considered. With this new renumbering, it can be seen that there are 14 detuned configurations amongst all the 216 detuned configurations (representing 6.5 % of all the detuned configurations) for which the nonlinear mistuned response has smaller response levels than those obtained with the pure mistuning case. Table 5 summarizes the results obtained with these improving detuned configurations. Note that there is only one configuration that presents a subcyclicity with order s = 2. The other improving detuned configurations have no particular cyclicity.

Nonlinear detuning optimization in presence of mistuning j j s nb 1 nb 2 w c, j q c, j deterministic detuning q c, j det 2 [1 0 0 0 0 0 1 0 0 0 0 0] 0.9995 0.9788 14 [START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF] 6 6 [0 0 1 1 0 1 1 0 1 1 0 0] 0.9999 0.9903 A careful attention is also paid to the 10 detuned configurations that present a cyclic symmetry s ∈ {2, 3, 4, 6} < N. The characteristics are summarized in Table 6 and the corresponding amplification factors are represented by symbols in Fig. 17. As shown in Table 5, except for de-j j s nb 1 nb 2 w c, j q c, j 87 127 6 Figure 17: Dynamical amplification factor according to detuned configuration. Graph of function j → q c, j for the n c = 352 detuned configurations (upper graph) and for the n c = 216 detuned configurations with a number of blades with type 2 less than or equal to 6 (lower graph). Subcyclicity order s is given for s = 2 (purple square symbol), s = 3 (red diamond symbol), s = 4 (orange bullet symbol) and s = 6 (green triangle symbol).

tuned configuration = 211 with s = 2 that has an amplification factor q c,211 = 0.9995 1 that slightly inhibits the mistuning effects, these subcyclic configurations give rise to an amplification between 5% and 16% with respect to the pure mistuning situation.

3.13. Nonlinear analysis of the mistuned response of the best and the worst detuned configurations The results are analyzed for the following detuned configurations: the tuned configuration corresponding to the pure mistuning case ( = 49), the best detuned configuration ( = 104), and the worst detuned configuration ( = 9). Fig. 18 displays the graphs of k → a ,k = A (θ k ) for k = 1, . . . , n sim , which describes the maximum displacement response over the frequency and the blades for each mistuning realization θ k of the detuned configurations ∈ {49, 104, 9}. By comparing these graphs, it can be seen that the best detuned configuration is characterized by slightly lower response levels but also by less scattered realizations with respect to the pure mistuned case (left figure). The realizations of the worst detuned configurations (right graph) are clearly more scattered with higher response levels. The distribution of the blades related to these detuned configurations are displayed in Fig. 19. For a given detuned configuration and a fixed frequency ν, we introduce the real-valued random variable Y max (2πν) whose realization θ k is defined by

Y max (2πν, θ k ) = Y ( j ,k max , 2πν, θ k ) , in which j ,k
max is defined by Eq. [START_REF] Castanier | Investigation of the combined effects of intentional and random mistuning on the forced response of bladed disks[END_REF]. For detuned configurations ∈ {49, 104, 9}, Figs. 20 to 22 display the graphs of the confidence region of the real-valued random function ν → Y max (2πν) for a probability level of 0.98. As expected, it can be seen that the structural displacements are mainly located in the frequency band of excitation B e but that there are also new resonances occurring below the frequency band of excitation B e that are induced by the geometrical nonlinearities. The amplitudes of these new resonances correspond to the indirect excitation of the first blade modes of the blisk and are of a lower order of magnitude in the present case, probably because the blades are not very slender. Note that this phenomenon has already been observed in previous analyses of nonlinear mistuned blisk for which the blades were slender structures, attached to the disk [START_REF] Capiez-Lernout | Mistuning analysis and uncertainty quantification of an industrial bladed disk with geometrical nonlinearity[END_REF][START_REF] Picou | Robust dynamic analysis of detuned-mistuned rotating bladed disks with geometric nonlinearities[END_REF]. Moreover, the phenomenon for which the amplitude of the resonances located below the frequency band of excitation are dominant has been investigated in another context of sloshing in nonlinear fluid structure interaction [START_REF] Akkaoui | Revisiting the experiment of a free-surface resonance of a liquid in a vibration tank using a nonlinear fluid-structure computational model[END_REF]. Fig. 23 shows a zoom of these confidence regions around the excitation frequency band B e . It is clearly seen that the pure mistuned blisk have slightly more robust resonances that the one obtained with the detuned configurations in presence of mistuning. But the real interest is to verify that the optimal detuned configuration (middle graph) guarantees that the upper envelope of the confidence region stays below the critical level of response of the pure mistuned response. It is thus observed that the optimal detuned configuration yields a reduction of the maximum amplitude of about 12% with respect to the initial situation that corresponds to the pure mistuning case. Another observation is that the worst detuned configuration yields an amplification of about 11% with respect to the initial situation that corresponds to the pure mistuning case. Note that the resonance involved does not correspond to the main resonance of the pure mistuned system. All these observations allow for assessing the efficiency of the optimization through the output of interest q c, . 1000 2000 3000 4000 5000 6000 7000 10 3.14. Sensitivity analysis of the nonlinear mistuned response of the best and the worst detuned configurations with respect to the detuning rate and the mistuning rate It is interesting to perform a sensitivity analysis with respect to the detuning rate for the 3 configurations studied in Section 3.13. The detuning rate is described by the ratio r defined by It quantifies a ratio for the blade frequency. In this Section, observation q c, is parameterized as a function of r and δ K and is rewritten as q c, (r; δ K ). All the optimization process has been carried out for the constant mistuning defined by δ K = 0. q c, (r = 4.34 %, 0.1) = 0.9394. Note that this optimum is not necessarily a global optimum since the optimization has been carried out for r = 5.13 %. Moreover, it can also be seen that detuned configuration = 104 does not behave monotonically with respect to inhibition or amplification induced by the mistuning. There is a range of detuning rate r ∈ [3.30, 6.99] % for which the detuned configuration reduces the amplification induced by the mistuning. But we have to carefully pay attention not to have a detuning rate r smaller than 3.30 %. Note also that the detuned configuration = 9 corresponding to the worst detuned configuration in the optimization process increases the amplification induced by the mistuning whatever the detuning rate r. Then a sensitivity analysis is performed for these three detuned configurations with respect to the mistuning rate δ K , which is the hyperparameter that calibrates the uncertainty level in meta-stiffness matrix (See Section 2.4), and which contains both linear and nonlinear stiffness contributions. It is interesting to relate it to physical considerations. For instance, in the pure mistuning case and for a mistuning rate δ K = 0.1, the confidence region obtained with a probability level of 98% and related to the location of the main resonance is described by a dispersion between -1.76 % and +5.59 % around the main resonance. The performed sensitivity analysis involves a mistuning rate from δ K ∈ [0, 0.24]. The graph δ K → q c, (5.13 %; δ K ) is shown in Fig. 25. The purple line describes the pure mistuning case taken as a reference and yielding q c,49 = 1. The detuned configuration = 9 never improves the mistuning situation whatever the mistuning rate. Concerning detuned configuration = 104 that is optimal for a mistuning rate δ K = 0.1, there is a threshold value δ K = 0.08 from which the mistuning amplification effects are inhibited with respect to the pure mistuning situation. A remark concerns the deterministic detuning optimization. It can be seen that neither detuned configurations = 9 or = 104 allows for inhibiting the mistuning amplification effects as deterministic detuning is considered. Although deterministic optimization Figure 24: Analysis of the sensitivity of mistuned response for given detuned configurations with respect to the detuning rate r: graph of r → q c, (r; 0.1) for the pure mistuned situation corresponding to = 49 (square symbol, purple line), for the best detuned configuration = 104 (diamond symbol, green line), and for the worst detuned configuration = 9 (circle symbol, orange line). Parameters for optimization process (vertical red line).

r = 1 - Young 2 Young 1 . ( 39
)
has no meaning because mistuning is an unavoidable phenomenon, it has also been performed and the results show that the optimal solutions are not optimal for the detuning optimization in presence of mistuning. Let q c, det be the corresponding quantity of interest related to q c, when considering deterministic detuning (that is to say without mistuning). Table 7 displays the seven optimal deterministic detuned configurations and compares the quantity of interest when detuning optimization is carried out in presence of mistuning. It can be seen that three optimal deterministic detuned configurations belong to the set of the fourteen optimal detuned configurations (in presence of mistuning) but that such detuned configurations are in fact part of the less optimal configurations since the corresponding quantity q c, is close to 1. Table 5 also compares the quantity of interest of the detuning optimization with those obtained when no mistuning is considered. It is interesting to see that such quantity is sensitive to the presence of mistuning. This proves that mistuning cannot be neglected in the detuning optimization process. The most optimal detuned configurations are not optimal solutions when no mistuning is considered in the detuning optimization process.

Nonlinear deterministic detuning optimization (without mistuning)

j j s nb 1 nb 2 w c, j deterministic detuning q c, j det q c, j 39
Figure 25: Analysis of the sensitivity of mistuned response for given detuned configurations with respect to the mistuning δ K : graph of δ K → q c, (5.13 %; δ K ) the pure to = 49 (square symbol, purple line), for the best detuned configuration = 104 (diamond symbol, green line), and for the worst detuned configuration = 9 (circle symbol, orange line). Parameters for optimization process (vertical red line).

4. Optimization of the detuning in presence of mistuning using probabilistic learning on manifolds (PLoM) as a machine learning tool

Discrete set N c of a configuration

We consider the detuning of a blisk as the one analyzed in Section 3 having n w blades and two types of blades that are labelled by integer 0 and 1. A configuration (also called a pattern) of the blisk consists in giving a vector

w c = (w c 1 , . . . , w c n w ) , (40)
in which for k ∈ {1, . . . , n w }, w c k is equal to 0 or to 1. We then define the finite subset N c of N n w such that

N c = Π n w k=1 {0, 1} n w ⊂ N n w , (41)
and consequently, w c ∈ N c .

Discrete set C c of the n c possible configurations

The total number of all the possible configurations is denoted by n c . As already explained in Section. 2.5, n c can be very large as soon as n w is large. By excluding the rotations by an integer number of blades, the asymptotic value of n c is n c = (2 n w )/n w . For instance, if n w = 48, then n c = 5.86 × 10 12 . This asymptotic formula includes the configuration with n w blades of type 0 and the configuration with n w blades of type 1. If the configurations are restricted to the case for which the number of blades of type 1 is less than or equal to n w /2, then n c < (2 n w )/n w . It will be the case for the configurations of the detuned blisk considered in this section, for which n w = 12 and therefore, n c = 216 (see also Section 3.9).

The finite set C c of the n c possible configurations is defined by

C c = {w c,1 , . . . , w c,n c } , w c, ∈ N c . (42)
4.3. Quantity of interest (QoI) computed for a given configuration For a given configuration w c ∈ N c , the HFCM is used (see Section 3) for computing the QoI that is the dynamic amplification factor q c ∈ R (defined by Eq. ( 38)) of the detuned disk in presence of random mistuning. Consequently, there exists a mapping w c → f HFCM (w c ) defined on N c with values in R such that q c = f HFCM (w c ) .

Mapping f HFCM is completely defined by the set {(q c, , w c, ), = 1, . . . , n c } of the n c points in R × N c in which q c, = f HFCM (w c, ) with w c, ∈ C c .

w c ∈C c J c (w c ) , (44)
in which w c → J c (w c ) is chosen for minimizing the dynamic amplification factor defined by Eq. ( 43), that is to say,

J c (w c ) = f HFCM (w c ) , w c ∈ N c (45)
(i) The exhaustive search of the combinatorial optimization problem on finite set defined by ( 44)- [START_REF] Sanchez-Alvarez | Key action mechanisms of intentional mistuning[END_REF], is not tractable and unachievable for large value of n c . There are many methods in the literature for solving such a problem: the theory of linear programming [START_REF] Karloff | Linear Programming[END_REF] including the polynomial-time algorithms [START_REF] Renegar | A polynomial-time algorithm, based on Newton's method, for linear programming[END_REF], integer programming [START_REF] Williams | Integer programming[END_REF] that deals with graph structures [START_REF] Diestel | Graph Theory[END_REF], the operations research [START_REF] Taha | Operations Research: an Introduction[END_REF], and the computational complexity theory [START_REF] Arora | Computational Complexity: a Modern Approach[END_REF]. Certainly, the most popular methods are the branch-and-bound [START_REF] Brusco | Branch-and-Bound Applications in Combinatorial Data Analysis[END_REF], the branch-and-cut [START_REF] Mitchell | Branch-and-cut algorithms for combinatorial optimization problems[END_REF], the tabu search and greedy [START_REF] Gendreau | Tabu search[END_REF], and the dynamic programming [START_REF] Aboueisha | Extensions of Dynamic Programming for Combinatorial Optimization and Data Mining[END_REF].

(ii) When the numerical cost for the evaluation of J c (w c ) at any point w v in C c is negligible, then the previous listed methods and algorithms can be used for estimating an optimal solution.

(iii) The specificity of the combinatorial optimization problem under consideration is the high numerical cost of one evaluation q c, = f HFCM (w c, ) using HFCM. In practice, only N d n c evaluations can be carried out (for instance, N d = 50 or 100 while n c can be 10 12 ). We therefore need to construct a metamodel of J c for computing J c (w c ) for any w c proposed by the optimization algorithm in order that each evaluation has a negligible numerical cost. We propose to use a machine learning tool for evaluating J c (w c ) for any w c in C c .

Transforming the combinatorial optimization problem in another one adapted to probabilistic machine learning tools

The combinatorial optimization problem defined by Eqs. ( 44)-( 45) is reformulated in a probabilistic framework, which will be adapted to a machine learning solver in order to limit the number of evaluations using HFCM.

Let W c = (W c 1 , . . . , W c n w ) be the N c -valued random variable, defined on a probability space (Θ, T , P), for which the n c elements w c,1 , . . . w c,n c of C c are n c independent realizations. Let Q c be the real-valued random variable, defined on (Θ, T , P) such that

Q c = f HFCM (W c ) . (46)
Consequently, {q c, = f HFCM (w c, ), = 1, . . . , n c } is the set of the n c realizations of Q c . The optimization problem defined by Eqs. ( 44)-( 45) is then replaced by the following one: find w c,opt ∈ N c such that w c,opt = arg min

w c ∈C c J c (w c ) , (47)
in which the cost function w c → J c (w c ) is chosen as the conditional mathematical expectation of

Q c given W c = w c , J c (w c ) = E{Q c | W c = w c } , w c ∈ N c . (48)
For w c, ∈ C c , the conditional random variable {Q c W c = w c, } is equal to the deterministic value q c, = f HFCM (w c, ). We thus have J c (w c, ) = J c (w c, ) for ∈ {1, . . . , n c }, which proves that Eqs. ( 47)-( 48) is a rewriting of Eqs. ( 44)-( 45) for w c ∈ C c . Therefore, this optimization problem is still a combinatorial optimization problem, but J c (w c ) will be computed using a probabilistic machine learning tool. Training set D d is made up of N d points (q c, , w c, ) such that q c, = f HFCM (w c, ) for ∈ {1, . . . , N d } with N d n c . This means that the construction of D d requires N d calls to HFCM, and consequently, N d will be small. The training set is therefore, by hypothesis, a small data set (in opposite to a big data set that is generally requires for the machine learning algorithms).

(ii) Learned set D ar . The second step will thus consist in using a machine learning method for generating, from the training set D d , the learned set D ar made up of N ar learned realizations {(q ar , w ar ), = 1, . . . , N ar } of an R × R n w -valued random variable (Q, W) (defined in Section 4.7) without calling HFCM. The number N ar of points is chosen sufficiently large (N ar N d ) to obtain a converged estimate of the conditional expectation introduced in Eq. ( 48) for computing J c (w c ) for any w c ∈ C c (see Eq. ( 42)). So D ar can be viewed as a big data set. 4.7. Probabilistic learning on manifolds (PLoM) used as a machine learning tool and approximate combinatorial optimization problem There are many methods available in the literature for building metamodels. Each method has its own validity assumptions with respect to the input-output dimensions, the complexity of the function to be represented, the number of points in the training set (small or big), the presence or not of uncertainties, a formulation performed in a deterministic or stochastic framework. It is not possible here to propose a review of all these methods, knowing that the use of each requires a precise analysis of its domain of validity, which depends on the hypotheses with which it was constructed. Nevertheless, we refer the reader to [START_REF] Ghanem | Handbook of Uncertainty Quantification[END_REF] for methods based on probability theory and mathematical statistics, including polynomial chaos expansion methodology, to [START_REF] Zimmermann | Gradient-enhanced surrogate modeling based on proper orthogonal decomposition[END_REF][START_REF] Constantine | Active subspace methods in theory and practice: applications to kriging surfaces[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF][START_REF] Bhosekar | Advances in surrogate based modeling, feasibility analysis, and optimization: A review[END_REF][START_REF] Kalogeris | Diffusion maps-based surrogate modeling: An alternative machine learning approach[END_REF][START_REF] Mignolet | Compressed principal component analysis of non-gaussian vectors[END_REF] for surrogate based modeling, to [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Carlberg | A low-cost, goal-oriented compact proper orthogonal decomposition basis for model reduction of static systems[END_REF][START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF][START_REF] Amsallem | Design optimization using hyper-reduced-order models[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Soize | Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics[END_REF][START_REF] Swischuk | Projection-based model reduction: Formulations for physics-based machine learning[END_REF][START_REF] Wang | Discussion on a non-intrusive model-order reduction of geometrically nonlinear structural dynamics using modal derivatives[END_REF] for projectionbased model reduction, and to [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design[END_REF][START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][START_REF] Byrd | On the use of stochastic Hessian information in optimization methods for machine learning[END_REF][START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF] for optimization of expensive functions.

The specific difficulty of the problem considered here is the small number of points available in training set D d (a few tens or even hundreds for industrial turbomachinery problems) and the complexity of the function f HFCM . For instance, the very efficient ANN approach [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF] cannot be used, a priori, for constructing a metamodel of J c because N d is too small. There are methods for sampling underlying distributions on manifolds [START_REF] Talwalkar | Large-scale manifold learning[END_REF][START_REF] Gorissen | A surrogate modeling and adaptive sampling toolbox for computer based design[END_REF][START_REF] Öztireli | Spectral sampling of manifolds[END_REF][START_REF] Marzouk | Sampling via measure transport: An introduction[END_REF][START_REF] Parno | Transport map accelerated markov chain Monte Carlo[END_REF][START_REF] Perrin | Data-driven kernel representations for sampling with an unknown block dependence structure under correlation constraints[END_REF][START_REF] Tsilifis | Bayesian adaptation of chaos representations using variational inference and sampling on geodesics[END_REF][START_REF] Kevrekidis | Manifold learning for parameter reduction[END_REF]. Among all these existing methods in machine learning, there is the probabilistic learning method (PLoM), which has specifically been developed for small non-Gaussian data (small value of N d ) in arbitrary dimension [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF], with possibility to take into account additional constraints coming from experiments [START_REF] Soize | Physics-constrained non-gaussian probabilistic learning on manifolds[END_REF] or from nonlinear partial differential equations [START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF], to construct a polynomial chaos representation of databases on manifolds [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF], to construct Bayesian posteriors in high dimension [START_REF] Soize | Sampling of bayesian posteriors with a non-gaussian probabilistic learning on manifolds from a small dataset[END_REF], and which has been used for complex optimization problems under uncertainties [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF][START_REF] Ghanem | Optimal well-placement using probabilistic learning[END_REF] and challenging applications [START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF][START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF][START_REF] Arnst | Computation of sobol indices in global sensitivity analysis from small data sets by probabilistic learning on manifolds[END_REF].

(i) Generating the learned set D ar with the PLoM method. We thus propose to use the PLoM method for generating the learned set D ar from the training set D d . The PLoM method is adapted to a continuous random variable X = (Q, W), defined on (Θ, T , P), with values in R n = R × R n w with n = 1 + n w , while the combinatorial optimization problem defined by Eqs. ( 47)-( 48) is related to the random variable X c = (Q c , W c ), in which Q c = f HFCM (W c ) (see Eq. ( 46)) and where W c is a random variable with values in a finite set N c . From the training set D d = {(q c, , w c, ), = 1, . . . , N d }, the PLoM method is used for generating the learned set D ar = {(q ar , w ar ), = 1, . . . , N ar }.

In order to facilitate the reading of this paper, the reader will find in Section A.1 of Appendix A a summary of the PLoM algorithm. We give this summary, because the proposed algorithm is the assembly of ingredients, which are distributed in three different papers with slightly different notations: basic algorithm of PLoM [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF], novel algorithm to estimate the optimal value of the parameter of the kernel for the calculation of the diffusion-maps basis [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF], and taking into account of the normalization constraints [START_REF] Soize | Physics-constrained non-gaussian probabilistic learning on manifolds[END_REF].

(ii) Approximate combinatorial optimization problem. In order to use the learned set D ar , we introduce the following approximation of the combinatorial optimization problem defined by Eqs. ( 47)-( 48): find w opt in N c such that

w opt = arg min w c ∈C c J ar (w c ) , (49)
in which the learned cost function w c → J ar (w c ) is defined by

J ar (w c ) = E{Q | W = w c } , w c ∈ N c , (50)
in which (Q, W) is the random variable with values in R × R n w , which approximates (for the learning step) random variable (Q c , W c ) with values in R × N c and for which the learned realizations are {(q ar , w ar ), = 1, . . . , N ar }. Since random variable (Q, W) has been defined by its learned realizations, the joint probability density function p Q,W of (Q, W) on R × R n w will be constructed by using the multidimensional Gaussian Kernel Estimation (KDE) method with the learned realizations of D ar (thus this pdf exists by construction, see Section 4.8). With such a hypothesis, the right-hand side of Eq. ( 50) can be written as

E{Q|W = w c } = 1 p W (w c ) R q p Q,W (q, w c ) dq, (51)
in which p W (w c ) = R p Q,W (q, w c ) dq.

(H3) Definition of the best optimal configuration. Let n s be fixed and let w c, opt be the configuration identified by Eq. ( 55) . If we have

J c (w c, opt ) < min =1,...,N d f HFCM (w c, ) , (56)
then w c, opt is identified as the optimal configuration in W opt n s (N ar ), which is better than all the configurations {w c, , = 1, . . . , N d } used for constructing the training set D d .

(H4) If Eq. ( 56) does not hold, that is to say, a configuration has not be found for improving the knowledge corresponding to the training set, then we can:

(a) Increase the value of n s (N d remaining fixed at its value). Note that the value of n s should be kept small because the number of evaluations with HFCM is n s . (c) It should be noted that to build the learned set from the training set, the PLoM algorithm is initialized with the points of the training set, which means that the points of the training set belong also to the learned set. Consequently, if the best optimal configuration belongs to the training set, then this optimal configuration will be identified with Eq. ( 56).

Estimating the cost function using the learned set

As explained in Section 4.7-(i), the PLoM method allows for generating N ar N d realizations {x ar , = 1, . . . , N ar } with x ar = (q ar , w ar ) of random variable X = (Q, W) with values in R × R n w , n = 1 + n w , without calling HFCM. The classical Gaussian KDE method of the pdf x → p X (x) on R n of random variable X, with respect to the Lebesgue measure dx, is written [START_REF] Duong | Cross-validation bandwidth matrices for multivariate kernel density estimation[END_REF][START_REF] Duong | Feature significance for multivariate kernel density estimation[END_REF][START_REF] Filippone | Approximate inference of the bandwidth in multivariate kernel density estimation[END_REF][START_REF] Zougab | Bayesian estimation of adaptive bandwidth matrices in multivariate kernel density estimation[END_REF] as

p X (x) = 1 N ar N ar =1 1 ( √ 2πs SBX ) n σ 1 ...σ n × exp{- 1 2s 2 SBX n k=1 ( x k -x ar σ k ) 2 } , (57)
in which σ k is the empirical estimate of the standard deviation of X k calculated with {x ar,k , = 1, . . . , N ar } and where s SBX is the Silverman bandwidth that is written [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis: The Kernel Approach With S-Plus Illustrations[END_REF] as

s SBX = 4 N ar (2 + n) 1/(4+n) . (58)
Since W with values in R n w is a representation of W c with values in the finite set {0, 1} n w , a modification of Eq. ( 57) is proposed as follows.

(i) Random variable Q is normalized in a random variable Q. Let q and σ Q be the empirical mean value and standard deviation of Q estimated with {q ar , = 1, . . . , N ar }. We then have,

Q = q + σ Q Q , Q = (Q -q)/σ Q , (59)
whose realizations of Q are

q ar = (q ar -q)/σ Q , ∈ {1, . . . , N d } , (60)
Consequently, Eq. ( 50) can be rewritten as

J ar (w c ) = q + σ Q E{ Q | W = w c } , w c ∈ N c . (61)
In Eq. ( 61), the conditional mathematical expectation is written as

E{ Q|W = w c } = 1 p W (w c ) R q p Q,W ( q, w c ) d q, (62)
in which p W (w c ) = R p Q,W ( q, w c ) d q and where p Q,W is the pdf on R × R n w of random variable ( Q, W) with respect to d q dw, whose following construction is directly deduced from Eq. ( 57).

(ii) The PLoM algorithm that is summarized in Section A.1 of Appendix A is appropriated to the case of a random vector X that belongs to an uncountable set R n . In the present case, the quantity of interest Q has values in an uncountable set whereas the control variable W is a representation of W c that takes its values in a finite discrete set. This mixed uncountable -countable case can be approached elegantly with PLoM by taking an adapted value of the Silverman bandwidth of the nonparametric representation of the probability density function of W in order to properly separate the bimodality linked to the two discrete values {0} and {1}. This is done and explained below.

Since W c is a random variable with values in the finite set {0, 1} n w and since the R n w -valued random variable W is a representation of W c , we do not normalize W and we choose its bandwidth as s w = c 0 s SBW for constructing the Gaussian KDE of pdf p W in order to obtain a possible bimodality centered in 0 and in 1, in which s SBW is the following Silverman bandwidth of W,

s SBW = 4 N ar (2 + n w ) 1/(4+n w ) . (63)
The numerical experiments have shown that c 0 = 0.2 is a good value for obtaining a clear separation (see Fig. 26). Therefore, we have

p W (w) = 1 N ar N ar =1 1 (s w √ 2π) n w exp{- 1 2s 2 w w -w ar 2 }, (64)
in which . is the usual Euclidean norm on R n w .

(iii) Using (i) and (ii), the Gaussian KDE of p Q,W is defined as

p Q,W ( q, w) = 1 N ar N ar =1 1 s q √ 2π exp{- 1 2s 2 q ( q -q ar ) 2 × 1 (s w √ 2π) n w exp{- 1 2s 2 w w -w ar 2 } , (65)
in which s q is chosen as s q = (4/(3N ar )) 1/5 . From Eqs. ( 62) and ( 65), it can be deduced that 64) with s w = c 0 s SBW , c 0 = 0.2, and s SBW given by Eq. ( 63).

E{ Q|W = w c } = N ar =1 q ar exp{-1 2s 2 w w c -w ar 2 } N ar =1 exp{-1 2s 2 w w c -w ar 2 } . (66)
4.9. Data for the validation of the optimization method based on a probabilistic machine learning tool In subsections (i)-(iv), there are a few repetitions of results and discussions that have already be given in the previous sections. These small repetitions have voluntary be introduced in order to facilitate the reading of the paper in order to avoid going back and forth in the paper.

(i) Hypothesis and objective for the validation. The number n w of blades in the blisk (See Section 3) has deliberately be chosen so that the number n c of elements in C c (all the possible configurations) is also small in order to build the exact solution of the combinatorial optimization problem defined by Eqs. ( 44)- [START_REF] Sanchez-Alvarez | Key action mechanisms of intentional mistuning[END_REF], in order to be able to validate the proposed approach.

For n w = 12, the number of all possible configurations for which the blade of type 1 is less than or equal to n w /2 = 6 is n c = 216. Consequently, 216 calls to HFCM has been carried out for constructing the solution w c,opt un of Eqs. ( 44)- [START_REF] Sanchez-Alvarez | Key action mechanisms of intentional mistuning[END_REF]. We wish to insist here on the fact that the proposed optimization method has been developed to analyze bladed-disks whose number of blades is such that n c leads to a combinatorial optimization problem for which the cost function could only be evaluated with HFCM for N d configurations with N d n c . In such a case, the construction of the exact solution of Eqs. ( 44)-( 45) is unachievable and therefore, a reference solution would not be available for the validation. For instance, for n w = 24, we would have n c 350, 000, which would require 350, 000 call to HFCM, that is unachievable. (iii) Computation of the dynamic amplification factor for all the configurations in C c . In order to construct the exact optimal solution w c,opt un of the combinatorial optimization problem on finite set C c , defined by Eqs. ( 44)-( 45), the n c = 216 configurations {w c, , = 1, . . . , n c } of C c are considered and the n c dynamic amplification factors {q c, = f HFCM (w c, ), = 1, . . . , n c } (see Sections 3.12 and 3.13) have been computed with HFCM. Fig. 27 shows the graph → q c, for ∈ {1, . . . , n c }.

(1) The configuration w c,49 is the configuration without detuning in presence of random mistuning for which w c,49 k = 0 for k = 1, . . . , n w and q c,49 = 1.

(2) There are 14 configurations {w c, j , j = 1, . . . , 14} for which q c, j < 1. The three smallest values of q c, are reached for = 104, 123, and 166 and are q c, = 0.9476, 0.9529, and 0.9537.

(3) It can be deduced that there are 202 configurations for which q c, is greater than or equal to 1 (iv) Reference Optimal solution of the combinatorial optimization problem defined by Eqs. ( 44)-( 45) on finite set C c . From the previous paragraph (iii), it can be deduced that the reference optimal solution w c,opt un of Eqs. ( 44)-( 45) is w c,104 = [0 0 0 1 0 0 1 0 0 1 0 1], which is such that q c,104 = f HFCM (w c,104 ) = 0.9476.

(v) Defining the subset C c,q≥1 of C c for constructing the training set. In order to present a validation of the proposed method, which is based on the use of PLoM as a probabilistic learning tool, we consider the most difficult case that consists only of configurations w c, for which the dynamic amplification factor q c, = f HFCM (w c, ) is greater than or equal to 1. These configurations will constitute a subset C c,q≥1 defined by

C c,q≥1 = {{w c, j , j = 1, . . . , n c,q≥1 } | q c, j = f HFCM (w c, j ) ≥ 1} . (67)
Fig. 28 shows the graph → q c, for the n c,q≥1 configurations of C c,q≥1 . 52)) is carried out and the best solution w c, opt is deduced (see Section 4.7-(vi)).

D d = {(q c, j , w c, j ), j = 1, . . . , N d } ⊂ C c,q≥1 , (68) 43
(2) The third value N d = 100 is then considered in order to compare the corresponding best solution with the one obtained for N d = 75.

(3) Finally, the last value, N d = n c,q≥1 = 202, has just been considered to check the consistency of the approach proposed for solving the combinatorial optimization problem.

(vii) Values of the parameters of the PLoM algorithm. These values are given in Section A.2 of Appendix A.

Results and discussion

(i) Convergence analysis of the subset of optimal solutions with respect to the number N ar of points in the learned set D ar . This convergence analysis with respect to N ar (see Section 4.6-(ii)) has been defined in Section 4.7-(v) and is carried out for each training set D d defined by N d points. Let {(q ar , w ar ), = 1, . . . , N ar } be the N ar points of D ar generated by PLoM. The criteria for the convergence analysis are the empirical estimates conv Q (N ar ) and conv W (N ar ) of the mean-square norm of Q and W, 1.089133 (ii) Illustration of the detailed results for N d = 75 using PLoM. In order to limit the number of figures, we illustrate the results obtained by using PLoM for N d = 75. Qualitatively, the results for the three other values of N d are similar. (a) Generation of the learned set with PLoM. Figs. 31 and 32 show realizations q ar and w ar,1 (first component of w ar ) for ∈ {5 240 000, . . . , 5 245 000} of D ar generated by using PLoM (we have presented a zoom in order to avoid a view having an absence of space resolution when the 10 7 realizations are plotted). The pdf w 1 → p W 1 (w 1 ) of component W 1 of W, estimated using {w ar,1 , = 1, . . . , 10 7 } is shown in Fig. 33.

conv Q (N ar ) = { 1 N ar N ar =1 (q ar ) 2 } 1/2 , ( 69
) conv W (N ar ) = { 1 N ar N ar =1 ar 2 } 1/2 . ( 70
(b) Cost function computed with the learned set. Fig. 34 shows the graph of function → J ar (w c, ) for ∈ {1, . . . , n c } (with n c = 216) in which J ar (w c, ) is computed by using Eqs. ( 61) and [START_REF] Williams | Integer programming[END_REF] with N ar = 10 7 . This graph is superimposed to the graph of function → q c, (the reference shown in Fig. 27). It can be seen that the two graphs have similar fluctuations.

(c) Computation of the best optimal solution. For n s = 6, the indices { j , j = 1, . . . , 6} of the 6 configurations w c, j , which define subset W opt n s of optimal solutions (see Eqs. ( 52) and ( 53)) are 7, 34, 49, 55, 1, and 166, for which the corresponding values of q c, j (see Section 4.9-(iii) and Fig. 27) are respectively, 0.9974, 0.9882, 1.0000 (the tuned configuration), 1.0015, 1.0067, and 0.9529. The use of Eq. ( 56) yields the best optimal solution w c, opt = [1 0 1 1 0 0 1 0 0 1 0 0] for which opt = 166 and q c, opt = 0.9529.

(d) Numerical cost. The numerical cost in terms of Elapsed Time (in seconds) of the analysis performed with a workstation 3072 GB RAM and 110 cores (Intel(R) Xeon(R) Platinium 8280 CPU@2.706 Hz) is the following: 21 607 s for the generation of the learned dataset by the PLoM with constraints and 1 144 s for computing the best optimal solution. (iii) Checking the stability and consistency of the proposed methodology and algorithm. The proposed methodology and algorithm allow for estimating the best optimal solution w c, opt given by Eq. ( 56), in which W opt n s is constructed using Eqs. ( 52)-( 53), and where J ar (w c, j ) is computed by Eqs. ( 61) and [START_REF] Williams | Integer programming[END_REF], in which the learned set D ar = {(q ar , w ar ), = 1, . . . , N ar } is generated using PLoM algorithm from the training set D d defined by Eq. ( 68) for a given value of N d . In order to check the stability/consistency of the algorithm, we have analyzed the case for which and9 for n s = 6.

(1) The best optimal solution obtained for N d = 50 improves the dynamic amplification factors {q c, j , j = 1, . . . , 50} used in the training set D d defined by Eq. ( 68) for which q c, j ≥ 1, which have been computed with HFCM. However, although q c, opt = 0.9822 < 1, it is not a value close to the optimum that is 0.9476. As explained in Section 4.7-(vi)-(H4)-(a) and (b), this solution can possibly be improved by increasing the value of N d as it can be seen in Table 8.

(2) For N d = 75 and, a fortiori, for N d = 100 and 202, the best optimal solutions q c, opt are 0.9529 and 0.9538 that are close to the optimum 0.9476, and which constitute very good predictions. 133 104 0.9476 216 q c, j 0.9476 0.9529 0.9538 0.9822 0.9834 0.9824 Table 9: Best optimal configuration w c, opt listed in Table 8. opt q c, opt w c, opt = (w c, opt 1

, . . . , w c, opt 12 ) 104 0.9476 0 0 0 1 0 0 1 0 0 1 0 1 123 0.9537 1 0 0 1 0 0 0 0 0 0 0 0 141 0.9822 1 0 0 0 0 1 0 0 0 0 0 0 166 0.9529 1 0 1 1 0 0 1 0 0 1 0 0

Conclusion and discussion

We have presented an approach for the optimization of the detuning in presence of random mistuning and geometrical nonlinearities for bladed-disks, based on the use of high-fidelity computational models. This very challenging problem has given rise to very little published work and remains an open subject. The difficulties we have addressed are related to the developments of an efficient computational methodology for reducing the computational cost, to the physics understanding of such stochastic nonlinear dynamical systems, and to the probabilistic formulation of the detuning optimization problem. We have presented a deep computational analysis on a bladed-disk that is representative of industrial turbomachines, in order to understand the role played by the geometrical nonlinearities on the dynamical behavior and to exhibit the consequences on the detuning effects. Several results can be put forward.

A careful convergence analysis with respect to all numerical parameters must be performed to obtain a predictive solution with an optimal computational cost. This is essential for exploring a relatively large number of detuned configurations in the context of the stochastic nonlinear mistuning analysis. It should be noted that the sensitivity analysis with respect to the time step used in the numerical resolution of the nonlinear coupled differential equations is particularly delicate, yielding a fast convergence in the excitation frequency band but a slow convergence in the low-frequency range that is not directly excited by the external loading.

Another point of attention concerns the choice of the quantities of interest used in the formulation of the detuning optimization problem from a scalar-valued highly nonlinear cost function. It is found that there is a few number of optimal solutions with respect to the number of possible detuned configurations. It is also proved that the optimization problem is well-posed, yielding to robust optimal detuned configurations that are identified after post-processing analysis as detuned configurations for which the mistuning amplification effects are inhibited with respect to the pure mistuning situation.

The detuning patterns yielding the best optimal detuned configurations and yielding the worst amplification response levels have no particularly structure in terms of number of blades of different types and of blade distribution, which deserve further investigations to understand the complex mechanisms induced by the detuning. As expected, the sub-cyclic detuned configurations that present a cyclic symmetry with a lower cyclic order are sensitive to the mistuning and are not part of the optimal solutions.

For the considered bladed-disk, the sensitivity analysis with respect to the detuning rate or the mistuning rate shows that the detuned configuration yielding the worst case in terms of response amplification remains a wrong detuned configuration whatever the detuning or mistuning rate. There also exists a threshold in terms of detuning rate and mistuning rate from which a detuned configuration belongs to the set of optimal detuned configurations. There is no guarantee that an optimal solution remains optimal when the detuning rate increases too much. Note that the mistuning rate corresponds to the uncertainty level related to the global stiffness matrix with a huge dimension and that the involved magnitude orders are much larger than the uncertainty levels describing the usual frequency mistuning modeling but corresponds in fact to small mistuning levels.

The chosen external load corresponds to a narrow-frequency band of excitation so that the linearized detuning optimization problem has no real interest. It cannot be properly achieved since the natural frequencies of most detuned configurations do not belong to this narrow excitation frequency range and since the structure only can respond in the excitation frequency band. This is to underline that the nonlinear detuned optimization problem is of different nature because of the complex vibratory behavior induced by the geometric nonlinearities that appear themselves as an excitation that spreads the nonlinear response on the whole frequency band of analysis.

Concerning the detuning optimization of bladed-disks in presence of random mistuning, for nonlinear responses regimes involved by geometrical nonlinearities and using high-fidelity computational models, there are no published methods. When the number of detuned configurations is large, we have proposed a reformulation in a probabilistic framework of the combinatorial optimization problem, which is adapted to a probabilistic machine learning tool in order to limit the number of evaluations of the cost function with the high-fidelity computational model. The optimization must be carried out in a probabilistic framework in order to give robustness to the identification of an optimal detuned configuration. It should be noted that the proposed approach is unambiguous as it allows for building a subset of optimal detuned configurations and the search for the optimal solution in this subset is done using the high-fidelity computational model. The methodology proposed has been validated for a 12-bladed-disk structure for which the exact optimal detuned configuration in presence of random mistuning has been identified. A good prediction has been obtained for this very difficult case, which constitutes one validation of the proposed method.

Figure 1 :

 1 Figure 1: Reference configuration for i = 1, 2, 3 and using the classical convention for summations over repeated indices, as

Figure 2 :

 2 Figure 2: Finite element mesh of the blisk: blade sector (upper figure), full blisk (lower figure)

Figure 3 :Figure 4 :

 34 Figure 3: Critical damping rate: graph of ν → ξ(2πν). Square symbol is related to fundamental frequency ν 0 = 977.32 Hz and ξ 0 = ξ(2πν 0 ) = 0.00754.

Figure 5 :

 5 Figure 5: Graph of the modal shapes related to double eigenfrequencies ν 2,4 = 4 845.18 Hz (upper figures corresponding to blade dominant motion) and to ν 2,5 = 5 091.09 Hz (lower figures corresponding to blisk global motion)

Figure 6 :Figure 7 :

 67 Figure 6: Graph of convergence function S → Conv 1 (S ).

Figure 8 :

 8 Figure 8: Convergence analysis with respect to time sampling δt. Graph of function δt → Conv 2 (δt). Let x tip (t) = (x tip 1 (t), x tip 2 (t), x tip 3 (t)) be the displacement vector of the node located at the tip of blade number 0 and let x tip i (2πν) be its Fourier transform. Fig. 9 displays the graphs of t → x tip 3 (t) and ν → |x tip 3 (2πν)|.First of all, it can be seen that the structure returns to its equilibrium state at the end of the simulation. The time response is strongly sensitive to the choice of δt. The response in the frequency domain displays consequent nonlinear effects, since the structure responds outside the frequency band of excitation B e . The convergence with respect to the choice of δt is relatively quick when analyzing the structure in the excitation frequency band. The fundamental frequency of the blisk is also excited through the nonlinearities and it can be seen that this response is strongly over-estimated when choosing δt = 62.5×10 -6 s and becomes stable from for time steps lower than δt = 15.6×10 -6 s, which is confirmed by the zoom of the nonlinear response in the low-frequency range displayed in Fig.10. From now on all the computations in the time domain are done taking δt = 15.6×10 -6 s.

Figure 9 :Figure 10 :

 910 Figure 9: Graph of function t → x tip 3 (t) representing the nonlinear response at the tip of blade 0 in the transverse direction for different values of δt (upper figure). Graph of ν → |x tip 3 (2πν)| in the frequency analysis B = [300, 6 000] Hz.

Figure 11 :

 11 Figure 11: Convergence analysis with respect to the order M of the NL-ROM. Graph of function M → Conv lin 3 (M) corresponding to the linearized case (upper graph) and graph of M → Conv 3 (M) (lower graph).

Figure 12 :

 12 Figure 12: Sensitivity analysis with respect to load intensity s 0 . Graph of function s 0 → Conv 4 (s 0 ).

Figure 13 :

 13 Figure 13: Convergence analysis with respect to the number Q of eigenvectors to be kept in the simulation. Graph of function j → [Λ Q ] j j .

Figure 14 :Figure 15 :

 1415 Figure14: Estimation of the probability density function p ν 0,i (ν) for i = {1, 2, 3} of the first three natural frequencies with circumferential wave number h = 0 using kernel density estimation. Corresponding natural frequencies for the tuned case (red line).

type1 2 , 4 =

 24 4845.18 Hz and ν type1 2,5 = 5091.09 Hz. Concerning the tuned configuration with blades of type 2, the natural eigenfrequencies related to the veering with circumferential wave number h = 2 are the double eigenfrequencies ν type2 2,4 = 4596.54 Hz and ν type1 2,5 = 4829.83 Hz. Since the excitation frequency band B e is a narrow excitation frequency band chosen as B e = [4700 , 5200]

Figure

  Figure16: Graph of a zoom of the deterministic response for the nonlinear and linearized response of the tuned configurations with blades of type 1 and of type 2. Graphs of ν → y type 1 (2πν i ) (thick line), ν → y type 1,lin (2πν i ) (thin line), ν → y type 2 (2πν i ) (thick dashed line), ν → y type 2,lin (2πν i ) (thin dashed line). Excitation frequency band B e is represented by the light pink area.

5 Figure 18 :Figure 19 :

 51819 Figure 18: Analysis of the mistuned response for given detuned configurations: graph of realizations k → a ,k = A (θ k ) for the tuned configuration (pure mistuning) = 49 (left figure), the best detuned configuration = 104 (middle figure), and the worst detuned configuration = 9 (right figure).

Figure 20 :

 20 Figure 20: Analysis of the mistuned response for the tuned configuration: graph of the confidence region of ν → Y max (2πν) for a probability level of 0.98, corresponding to the pure mistuned configuration = 49 (purple area). Excitation frequency band B e is represented by the light pink area.

Figure 21 :

 21 Figure 21: Analysis of the mistuned response for given detuned configurations: graph of the confidence region of ν → Y max (2πν) for a probability level of 0.98, corresponding to the best detuned configuration = 104 (green area). Excitation frequency band B e is represented by the light pink area.

Figure 22 :

 22 Figure 22: Analysis of the mistuned response for given detuned configurations: graph of the confidence region of ν → Y max (2πν) for a probability level of 0.98, corresponding to the worst detuned configuration = 9 (orange area). Excitation frequency band B e is represented by the light pink area.

Figure 23 :

 23 Figure 23: Analysis of the mistuned response for given detuned configurations: for a probability level of 0.98, zoom of the confidence region of ν → Y max (2πν) for the pure mistuned configuration = 49 (upper graph), for the best detuned configuration = 104 (middle graph), and for the worst detuned configuration = 9 (lower graph). Critical level defined by the maximum of the upper envelope corresponding to the pure mistuned situation (red line).

4. 4 .

 4 Definition of a combinatorial optimization problem on finite set C c Such a combinatorial optimization problem [61, 62, 63] on finite set C c would consist in finding w c,opt un ∈ N c such that w c,opt un = arg min

4. 6 .

 6 Definition of the training set D d and learned set D ar for a probabilistic machine learning tool Taking into account Section 4.5, we will use a probabilistic machine learning tool for constructing a probabilistic metamodel of cost function J c in order to avoid a large number of evaluations of J c (w c ) with HFCM for w c ∈ C c . (i) Training set D d . The first step consists in generating the training set D d by using HFCM.

  (b) If (a) is not a success, this means that the value of N d is not large enough for obtaining a good learned set D ar from the training set D d . It is then necessary to increase the value of N d .

Figure 26 :

 26 Figure 26: Scheme of the graph of the pdf w k → p W k (w k ) of a component W k of W using Eq. (64) with s w = c 0 s SBW , c 0 = 0.2, and s SBW given by Eq. (63).

(

  ii) Defining the n c configurations of C c . The n c = 216 configurations {w c, , = 1, . . . , n c } of C c are extracted from the 352 configurations introduced in Section 3.9.

Figure 27 :

 27 Figure 27: For the n c = 216 configurations of C c , dynamic amplification factor q c, as a function of index ∈ {1, . . . , n c } of configuration w c, , computed with HFCM.

(

  vi) Defining the training set D d for the validation of the machine learning tool. Let N d be the number of points in the training set D d (see Section 4.6-(i). For the validation that we propose, 4 values of N d are considered. We recall that any convergence analysis with respect to N d has no object (see Section 4.7-(vi)-(H1)). The reason for this choice of these four values are the following. The construction of the training set D d is then performed as follows. For each considered value of N d , the first N d points of C c,q≥1 are selected,

Figure 28 :

 28 Figure28: For the n c,q≥1 = 202 configurations of C c,q≥1 , graph of the dynamic amplification factor q c, j = f HFCM (w c, j ) ≥ 1 as a function of index j ∈ {1, . . . , n c,q≥1 } computed with HFCM.

)For N d = 50 ,

 50 the convergence is obtained for N ar = 2 × 10 7 , while for N d = 75, 100, and 202, for N ar = 10 7 . For illustration and for N d = 75, Fig. 29 displays the graph of function N ar → conv Q (N ar ) defined by Eq. (69) while Fig. 30 displays the graph of function N ar → conv W (N ar ) defined by Eq. (70). The graphs of the three other values of N d are similar.

Figure 30 :

 30 Figure 30: For N d = 75, graph of function N ar → conv W (N ar ).

Figure 31 :

 31 Figure 31: For N d = 75, realizations q ar for ∈ {5.240, . . . , 5.245} × 10 6 of Q generated by PLoM.

Figure 32 :

 32 Figure 32: For N d = 75, realizations w ar,1 of the first component of w ar for ∈ {5.240, . . . , 5.245} × 10 6 of W generated by PLoM.

Figure 33 :

 33 Figure 33: For N d = 75, pdf w 1 → p W 1 (w 1 ) of component W 1 of W, estimated using {w ar,1 , = 1, . . . , 10 7 } generated by PLoM.

Figure 34 :

 34 Figure 34: For N d = 75, graph of function → J ar (w c, ) for ∈ {1, . . . , n c }, with n c = 216 (red thick line) and graph of function → q c, , the reference shown in Fig. 27 (black thin line).

Figure 36 :

 36 Figure 36: For N d = 75, graph of function α → λ α (ε opt ).

Figure 37 :

 37 Figure 37: For N d = 75, error function err(i) as a function of iteration number i of the iterative algorithm for the computation of the Lagrange multipliers.(v) Quantifying the concentration of the probability measure of random matrix [H m opt ]. The use of Eq. (75) yields d 2 N d (m opt ) = 0.073 2, which shows the concentration of the measure is kept.

Table 1 :

 1 3. Computational characteristics of the finite element model

		Elements Nodes	dof
	Sector	476	846	2151
	Full model	5712	9 036 27 108

Table 2 :

 2 Numerical parameters for Newmark algorithmLet u = (u 1 , . . . , u M ) be the vector of the M modal coordinates. Let δt → Conv 2 (δt) be the convergence function defined by Fig.8shows the graph of δt → Conv 2 (δt). It can be seen that a reasonable convergence is obtained for δt = 15.6×10 -6 s yielding n t = 8 192 time steps, ν e = 64 000 Hz, and a frequency resolution δν = 7.8 Hz.

	-3

M

-valued function t → u(t; δt).

Table 3 :

 3 Analysis of the most contributing generalized coordinates

			2	1	3
	ν(Hz)	4845.2 5091.1 4840.0 4865.6
	ratio	1.0000 0.3660 0.2108 0.1927
	nonlinear
	h	2	2	1	3	0	5	2	4	6	0
	ν(Hz)	5091.1 4845.2 4840.0 4865.6 5266.1 4872.1 2571.5 4871.3 4872.1 4768.6
	ratio	1.0000 0.8001 0.5421 0.3519 0.3129 0.2694 0.2653 0.2431 0.2375 0.2092

Table 4 :

 4 Number n c (nb 1 , nb 2 ) of detuned configurations as a function of the number nb 1 and nb 2 of blades of type 1 and 2. sector 1 by decreasing the Young modulus of the blade from Young 1 = 2.00 × 10 11 N/m 2 to Young 2 = 1.80×10 11 N/m 2 .

			11 10 9	8	7	6	5	4	3	2	1	0
	nb 2	0	1	2	3	4	5	6	7	8	9 10 11 12
	n c (nb 1 , nb 2 )	1	1	6 19 43 66 80 66 43 19 6	1	1

Table 5 :

 5 Characteristics of the the improving detuned configurations.

Table 6 :

 6 Characteristics of the detuned configurations with a subcyclicity s.

	1.2
	1.1
	1
	0.9
	0.8
	0	50	100	150	200	250	300	350
	1.2
	1.1
	1
								s=2
	0.9							s=3 s=4
								s=6
	0.8
	0	50	100	150	200	250	300	350

Table 7 :

 7 Characteristics of the the improving deterministic detuned configurations without mistuning

	0	0.05	0.1	0.15	0.2

  Figure 29: For N d = 75, graph of function N ar → conv Q (N ar ).

	1.0891355
	1.089135
	1.0891345
	1.089134
	1.0891335
	2.197
	2.1968
	2.1966
	2.1964
	2.1962
	2.196
	0	2	4	6	8	10
						10 6

Acknowledgments

The authors thank Pr. Christophe Pierre from the Illinois University concerning the use of the bladed-disk finite element model.

(iii) Difficulties induced by the fluctuations (variations) of the cost function. There is still a difficulty with the combinatorial optimization problem defined by Eqs. ( 49)- [START_REF] Desceliers | Non-linear viscoelastodynamic equations of three-dimensional rotating structures in finite displacement and finite element discretization[END_REF], induced by the physics of the nonlinear stochastic dynamics of detuned bladed-disks in presence of random mistuning. Function w c → f HFCM (w c ) defined for w c ∈ C c has many local minima with very close values (very similar dynamic amplification factor q c = f HFCM (w c ) for many different configurations w c ). This implies that function w → J ar (w) defined by Eq. ( 50) and extended on R n w , is not convex on R n w , and has also many local minima with very close values.

(iv) Reformulation of the combinatorial optimization problem. Taking into account the need to introduce the approximation defined by Eq. ( 50) to build a metamodel of cost function J ar and taking into account the specificities of its fluctuations (variations) involved by physics, the formulation defined by Eqs. ( 49)-( 50) is not suitable and must be adapted. Estimating a single optimal configuration w opt , which yields the smallest dynamic amplification factor, is not sufficiently robust. We therefore propose to reformulate Eq. ( 49) as follows.

Subset W opt n s of optimal solutions. Let n s n c be an integer that is fixed a priori (chosen as a few units). Let W opt n s be the subset of C c made up of the n s configurations w c, 1 , . . . , w c, ns in C c with { 1 , . . . , n s } ⊂ {1, . . . , n c }, which correspond to the n s first smallest values J ar (w c, 1 ), . . . , J ar (w c, ns ) of the set {J ar (w c ), w c ∈ C c }. We then have,

J ar (w c, 1 ) ≤ J ar (w c, 2 ) ≤ . . . ≤ J ar (w c, ns ) .

The subset of optimal solutions of the combinatory optimization problem is therefore defined as the subset W opt n s of C c . (v) Convergence analysis of the subset of optimal solution. Since the conditional mathematical expectation in Eq. ( 50) will be estimated with the N ar points of the learned set D ar , the subset of optimal solution W opt n s depends on N ar and will be rewritten as W opt n s (N ar ) when necessary. A convergence analysis will be carried out in Section 4.10 with respect to N ar N d .

(vi) Searching a best optimal configuration in the subset W opt n s of optimal solutions. The hypotheses for constructing the best optimal configuration are the following.

(H1) The number N d of points in the training set is fixed. This number corresponds to the maximum numerical effort that can be made for performing N d calls to HFCM. As N d is assumed to be small, convergence with respect to N d has no object. Once set N d , we seek to improve the values of the dynamic amplification factors that have been computed with HFCM for constructing the training set.

(H2) Let us assume that N d and N ar are fixed. As previously explained, taking into account the complexity of the variations of function f HFCM , which is not explicitly known in C c , but only point-by-point using HFCM, it is not possible to demonstrate that subset W opt n s (N ar ) contains an optimal configuration that improves the dynamic amplification factors computed for constructing the training set D d . The only way to do this is to compute again, with HFCM,

for the n s configurations of subset W opt n s (N ar ), which have been identified by PLoM (see Eqs. ( 52)-( 53)). Therefore, we define the configuration w c, opt with opt ∈ { 1 , . . . , n s } as the one that satisfies, w c, opt = arg min A. Probabilistic learning on manifolds (PLoM), its parameterization and numerical results

A.1. Summary of the PLoM algorithm

The PLoM approach [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF] starts from a training set D d made up of a relatively small number N d of points. For the supervised case, it is assumed that the training set is related to an underlying stochastic manifold related to a R n -valued random variable X = (Q, W) in which X, Q (quantity of interest), and W (control parameter) are R n -, R n q -, and R n w -valued random variables defined on a probability space (Θ, T , P) with n = n q + n w . Let U (uncontrolled parameter) be another R n u -valued random variable defined on (Θ, T , P). Random variable Q is written as Q = f(U, W) = F(W) in which the measurable mapping f is not explicitly known (unknown) and F is such that F = f(U, •) is a random mapping. The probability distributions of the vector-valued random variables W and U are assumed to be given. The stochastic manifold is defined by the unknown random graph {w, F(w)} for w belonging to an admissible set C w that is the support of the probability distribution of W. With the PLoM construction, its is assumed that this stochastic manifold cannot directly be described. Under these conditions, the non-Gaussian probability measure of X is concentrated in a region of R n for which the only available information is the cloud of the points of the training set. The PLoM method makes it possible to generate the learned set D ar whose n MC N d points (learned realizations) are generated by the non-Gaussian probability measure that is estimated from the training set. The concentration of the probability measure is preserved thanks to the use of a diffusion-maps basis that allows to enrich the available information from the training set. It should also be noted that the estimate of the unknown probability measure of X cannot be performed from the training set by using an arbitrary estimator. It must be parameterized in a manner that permits convergence to any probability measure as its number of points in the training set goes towards infinity. The PLoM method therefore does not only consist in generating points that belong to the region in which the measure is concentrated, but also allows these learned points to be realizations of the estimate probability measure with the convergence properties evoked above. The choice of the kernel density estimation method for estimating the probability measure of X from the training dataset guarantees that this required fundamental property is satisfied [START_REF] Soize | Probabilistic learning on manifolds[END_REF]. Considering all of these assumptions in the construction and in the methodology of PLoM makes PLoM as a general approach for small datasets for arbitrary non-Gaussian measures. Using the learned set D ar , PLoM allows for carrying out any conditional statistics such as w → E{Q|W = w} from C w in R n q , and consequently, to directly construct metamodels in a probabilistic framework.

The training set

The PLoM method allows for generating the learned set D ar made up of N ar N d learned realizations {x ar , = 1, . . . , N ar } of random vector X. As soon as the learned set has been constructed, the learned realizations for Q and W can be extracted as (q ar , w ar ) = x ar for = 1, . . . , N ar .

[X 1 , . . . , X N d ] be the random matrix with values in M n,N d , whose columns are N d independent copies of random vector X. Therefore, [x d ] is one realization of random matrix [X]. The normalization of random matrix [X] is performed by using a principal component analysis (PCA) of X allowing for defining the random matrix

in which [µ] is the (ν × ν) diagonal matrix of the ν positive eigenvalues of the empirical estimate of the covariance matrix of X (computed using

be smaller that ε PCA . If ν < n, then there is a statistical reduction.

(ii) Construction of a reduced-order diffusion-maps basis (ROB). To identify the subset around which the points of the training set are concentrated, the PLoM relies on the diffusion-maps method [START_REF] Coifman | Diffusion maps[END_REF][START_REF] Lafon | Diffusion maps and coarse-graining: A unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF]. This is an algebraic basis of vector space R N d , which is constructed using the diffusion maps. Let [K] and [b] be the matrices such that, for all i and j in {1, . . . ,

The eigenvector ψ 1 associated with λ 1 = 1 is a constant vector. For a given integer κ ≥ 0, the diffusion-maps basis {g 1 , . . . , g α , . . . , g N d } is a vector basis of R N d defined by g α = λ κ α ψ α . For a given integer m with 3 ≤ m ≤ N d , the reduced-order diffusion-maps basis of order m is defined as the family {g 1 , . . . , g m } that is represented by the matrix

) and [g m ] α = g α . This ROB depends on two parameters, ε DM and m, which have to be identified. It is proven in [START_REF] Soize | Probabilistic learning on manifolds[END_REF], that the PLoM method does not depend of κ that can therefore be chosen to 0.

For estimating the optimal values ε opt of ε DM and m opt of m, we use the algorithm proposed in [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF], which is summarized hereinafter. We have to find the value m opt ≤ N d of m and the smallest value

with an amplitude jump equal to an order of magnitude (a factor 10 as demonstrated in [START_REF] Soize | Probabilistic learning on manifolds[END_REF]) between λ m opt (ε opt ) and λ m opt +1 (ε opt ). This property means that we have to find m opt ≤ N d and the smallest positive value ε opt in order (1) to have λ 2 (ε opt ) < 1 (one must not have several eigenvalues in the neighborhood of 1) and ( 2) to obtain a plateau for λ 2 (ε opt ) to λ m opt (ε opt ) with a jump of amplitude 10 between λ m opt (ε opt ) and λ m opt +1 (ε opt ). A further in-depth analysis makes it possible to state the following criterion and algorithm to easily estimate ε opt and m opt . Let ε DM → Jump(ε DM ) be the function on ]0, +∞[ defined by

The algorithm is the following:

-set the value of m to m opt = ν + 1; -identify the smallest possible value ε opt of ε DM in order that Jump(ε opt ) ≤ 0.1 and such that Equation (73) be verified.

Remark concerning the choice of factor 10 for the jump of the eigenvalues. This choice is justified in [START_REF] Soize | Probabilistic learning on manifolds[END_REF], which gives the mathematical results in support of PLoM (we refer the reader to Theorem 7.8 of this paper and its Lemmas 7.5, 7.6, and 7.7 on which the proof is based). In this paper, it is proven that Eq. ( 73) must hold for preserving the concentration of the probability measure, which is quantified by the square of the L 2 -distance (see Eq. ( 75)), and which is estimated with the N ar N d learned realizations. Note that there are two important hypotheses in Eq. ( 73) that must be verified: one is the existence of the plateau for m ∈ {2, . . . , m opt } and the other one is the existence of a jump between m = m opt and m = m opt +1. In [START_REF] Soize | Probabilistic learning on manifolds[END_REF], it is proven that the plateau and the jump on the eigenvalues are directly related to the variations of the function m → ε DM (m) (denoted as m → ε d (m) in [START_REF] Soize | Probabilistic learning on manifolds[END_REF]), which has to be rapidly decreasing in m in the neighborhood of m opt by the lower integer values and which has to remain much lower than 1 for m ≥ m opt (see Fig. 1 of [START_REF] Soize | Probabilistic learning on manifolds[END_REF]). When Eq. ( 73) holds, then ε DM (m) 1 for m ≥ m opt , that is a fundamental property used for proving Theorem 7.8. This property ε DM (m) 1 for m ≥ m opt has been quantified by choosing one order of magnitude for the jump of the eigenvalues.

Let M opt = {m opt , m opt + 1, . . . , N d } in which m opt is the optimal value of m previously defined. Theorem 7.8 of [START_REF] Soize | Probabilistic learning on manifolds[END_REF] shows that min m∈M opt d 2

, which means that the PLoM method, for m = m opt and [g m opt ] is a better method than the usual one corresponding to

(v) Generation of learned realizations {η ar , = 1, . . . , N ar } of random vector H. The generation of learned realizations [z 1 ar ], . . . , [z n MC ar ] of random matrix [Z m opt ] is carried out by using the MCMC generator detailed in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF], which is based on a reduced-order Itô stochastic differential equation (ISDE) that is constructed as the projection on the ROB of the ISDE related to a dissipative Hamiltonian dynamical system for which the invariant measure is the pdf of random matrix [H] constructed with the Gaussian kernel-density estimation method and [η d ]. Let {([Z(t)], [Y(t)]), t ∈ R + } be the unique asymptotic (for t → +∞) stationary diffusion stochastic process with values in M ν,m opt ×M ν,m opt , of the following reduced-order ISDE (stochastic nonlinear second-order dissipative Hamiltonian dynamical system), for t > 0,

(

)] in M ν,N d is defined, for all k = 1, . . . , ν and for all j = 1, . . . , N d , by

p(u j ) = 1

(2) [W wien (t)] = [W wien (t)] [a] where {[W wien (t)], t ∈ R + } is the M ν,N d -valued normalized Wiener process.

(3) [N ] is the M ν,N d -valued normalized Gaussian random matrix that is independent of process [W wien ].

(4) The free parameter f 0 , such that 0 < f 0 < 4/ s ν , allows the dissipation term of the nonlinear second-order dynamical system (dissipative Hamiltonian system) to be controlled in order to kill the transient part induced by the initial conditions. A common value is f 0 = 4 (note that s ν < 1).

( 6) The learned realizations {x ar , = 1, . . . , N ar } of random vector X are then calculated (see Eq. ( 71)) by

(vi) Constraints on the second-order moments of the components of H if loss of normalization occurs. In general, the mean value of H estimated using the N ar learned realizations {η ar , = 1, . . . , N ar }, is sufficiently close to zero. Likewise, the estimate of the covariance matrix of H, which must be the identity matrix, is sufficiently close to a diagonal matrix. However, sometimes the diagonal entries of the estimated covariance matrix can be lower than 1 (for instance 0.8). Such a case is application-dependent. Normalization can be recovered by imposing constraints

in the algorithm presented in paragraph (v). For that, we use the method and the iterative algorithm presented in [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF] (that is based on Sections 5.5 and 5.6 of [START_REF] Soize | Physics-constrained non-gaussian probabilistic learning on manifolds[END_REF]). The constraints are imposed by using the Kullback-Leibler minimum cross-entropy principle. The resulting optimization problem is formulated using a Lagrange multiplier v = (v 1 , . . . , v ν ) associated with the constraints. The optimal solution of the Lagrange multiplier is computed using an efficient iterative algorithm. At each iteration, the MCMC generator detailed in paragraph (v) is used. The constraints are rewritten as E{h(H)} = b , in which the function h = (h 1 , . . . , h ν ) and the vector b = (b 1 , . . . , b ν ) are such that h k (H) = (H k ) 2 and b k = 1 for k in {1, . . . , ν}. To take into account the constraints in the algorithm of paragraph (v), Eq. ( 76) is replaced by the following one,

The iteration algorithm for computing v i+1 as a function of v i is the following,

)}] (the covariance matrix), and where α i is a relaxation function (less than 1) that is introduced for controlling the convergence as a function of iteration number i. For given i 2 ≥ 2, for given β 1 and β 2 such that 0 < β 1 < β 2 ≤ 1, α i can be defined by: 54

The convergence of the iteration algorithm is controlled by the error function i → err(i) defined by err

At each iteration i, E{h(H v i )} and [cov{h(H v i )}] are estimated by using the N ar learned realizations of H m opt (v i ) obtained by reshaping the learned realizations. (i) PCA of random vector X = (Q, W). We have n q = 1 and n w = 12. Consequently, X is a R nvalued random variable with n = 1 + 12 = 13. The PCA of X is constructed using the training set D d defined by Eq. ( 68), with a relative tolerance ε PCA = 10 -6 . Using Eq. ( 72) yields ν = n = 13 (so there is no statistical reduction).

(ii) Computation of the diffusion-maps basis. The graph of function ε DM → Jump(ε DM ) defined by Eq. ( 74) is shown in Fig. 35 and allows for identifying the optimal value ε opt = 19.6 of ε DM . The function α → λ α (ε opt ) of the eigenvalues of the transition matrix for ε DM = ε opt is shown in Fig. 36 and shows that the optimal value m opt of m is 14. The diffusion-maps basis is made up of the columns of matrix [g m opt ] ∈ M ν,m opt . (iii) Values of the parameters that control the integration scheme of the reduced-order Itô stochastic differential equation. The dissipation coefficient is f 0 = 4 and the integration step is ∆r = 0.018399. The number of realizations of random matrix [H] with values in M ν,N d is n MC = 133 400 yielding N ar 10 7 realizations for the R ν -valued random variable H.

(iv) Constraints. The constraints {E{H 2 k } = 1, k = 1, . . . , ν} are applied for generating the N ar realizations {η ar , = 1, . . . , N ar } of random vector H. The parameters of the relaxation function are β 1 = 0.01, β 2 = 0.2, and i 2 = 20. The graph of the error function i → err(i) defined by Eq. ( 78) is shown in Fig. 37. The convergence is obtained at iteration i = 45 for which