domain: stat math

text:
Remy Priem

Nathalie Bartoli

Youssef Diouane

Sylvain Dubreuil

Paul Saves

Rémy Priem

High-dimensional efficient global optimization using both random and supervised embeddings

I. Nomenclature

= Gaussian distribution σ(𝑙) = variance function of a Gaussian process conditioned by 𝑙 sampled points μ(𝑙) = mean function of a Gaussian process conditioned by 𝑙 sampled points 𝜽 (𝑙) = hyper-parameters of a Gaussian process conditioned by 𝑙 sampled points 𝛼 (𝑙) = acquisition function 𝑦 (𝑙) 𝑚𝑖𝑛 = minimum output of the objective function in a set of 𝑙 sampled points 𝑓 A = objective function reduced in the linear subspace A

A Artificial Intelligence Engineer for Aeronautical Systems, Technical Unit, remy.priem@intradef.gouv.fr B Professor, Mathematics and Industrial Engineering Department, youssef.diouane@polymtl.ca, AIAA member C Senior researcher, Information Processing and Systems Department, nathalie.bartoli@onera.fr, AIAA MDO TC Member. D Researcher, Information Processing and Systems Department, sylvain.dubreuil@onera.fr, AIAA Member. E PhD Student, Information Processing and Systems Department, paul.saves@onera.fr, Student AIAA Member.

1 𝑨 = transfer matrix A = linear subspace B = hyper-cube in the linear subspace 𝑇 = number of dimension reduction methods R = dimension reduction method 𝛾 = reverse application 𝒖 = vector of design variables in the linear subspace 𝑔 = constraint in the linear subspace

II. Introduction

M ultidisciplinary Design Analysis and Optimization (MDAO) [1] methodology is an incoming pillar in many industries especially in aircraft design. It considers a set of interacting disciplines required to design an aircraft at the beginning of the design process. For instance, one can tackle the aerodynamics, propulsion and structural coupling considering the design variables of these three disciplines. Because of the environmental pressure on aviation, new cutting edge aircraft architectures are studied in order to decrease drastically the aircraft consumption [START_REF] Schmollgruber | Multidisciplinary exploration of DRAGON: an ONERA hybrid electric distributed propulsion concept[END_REF]. In this way, it is not possible to perform MDAO procedures with low fidelity models developed for classical tube and wing aircraft configurations and designers have to use high fidelity models, even at the beginning of the process, to assess the performance of the regarded aircraft. The resulting MDAO procedures are eventually performed on computationally time consuming problems to find the best architecture and furthermore, disciplines are often black-boxes because of their complexity. This means that no information excepted the output values is available and therefore, the MDAO cannot be completed with classical gradient based or evolutionary optimization algorithms [START_REF] Schmollgruber | Multidisciplinary exploration of DRAGON: an ONERA hybrid electric distributed propulsion concept[END_REF].

In the last two decades, this issue is tackled with Bayesian Optimization (BO) framework [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF][START_REF] Močkus | On Bayesian Methods for Seeking the Extremum[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF] relying on iterative enrichments of surrogate models (i.e. a Gaussian Process (GP)) to seek for the optimum in a given design space Ω ⊂ R 𝑑 where 𝑑 ∈ N + is the number of design variables. Because of the characteristics of the aircraft design, like safety constraints and the important number of design variables, the resulting MDAO procedures need specific BO algorithms to be solved. Constrained Bayesian Optimization (CBO) [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Gelbart | Bayesian Optimization with Unknown Constraints[END_REF][START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF][START_REF] Diouane | TREGO: a trust-region framework for efficient global optimization[END_REF] deals with the constraints whereas High Dimensional Bayesian Optimization (HDBO) undertakes the high number of design variables [START_REF] Eriksson | Scalable Global Optimization via Local Bayesian Optimization[END_REF][START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF][START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF][START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF][START_REF] Kandasamy | High Dimensional Bayesian Optimisation and Bandits via Additive Models[END_REF] in these problems. HDBO problems can be mathematically described as follows:

min 𝒙∈Ω 𝑓 (𝒙), (1)
where 𝑓 : R 𝑑 ↦ → R is the objective function, the design space Ω = [-1, 1] 𝑑 is a bounded domain and 𝑑 ≫ 20 corresponds to a high number of design variables.

In the context of BO, two approaches are mainly investigated to handle the high number of design variables. The first one is typically based on a specific adaptation of the GP structure [START_REF] Kandasamy | High Dimensional Bayesian Optimisation and Bandits via Additive Models[END_REF][START_REF] Bouhlel | Efficient Global Optimization for High-Dimensional Constrained Problems by Using the Kriging Models Combined with the Partial Least Squares Method[END_REF][START_REF] Eriksson | Scalable Constrained Bayesian Optimization[END_REF] to deal with the large number of design variables, e.g., EGO-KPLS [START_REF] Bouhlel | Efficient Global Optimization for High-Dimensional Constrained Problems by Using the Kriging Models Combined with the Partial Least Squares Method[END_REF] where a partial least squares (PLS) method is used to reduce the number of the GP hyper-parameters. The second approach focuses on the use of dimension reduction methods to scale down the design space [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF][START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF][START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF], e.g., REMBO [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF] where a random linear embedding of the initial space is used. The dimension of the obtained subspace can be hence much lower than the original one. Most of the methods based on the REMBO paradigm have difficulties particularly to derive the new bounds for the reduced-dimension optimization problem [START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF]. An important computational effort in [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF] is in general needed to complete the optimization.

Two drawbacks of most existing HDBO methods are the computational effort needed to perform the optimization process and the construction of accurate bounds (over the subspace) in which the optimization is performed. In this paper, a new HDBO method, named Efficient Global Optimization coupled with Random and Supervised Embedding (EGORSE) is introduced to overcome the challenges previously detailed. First, the standard BO framework is described in Section III. Then, in Section IV, the proposed methodology is detailed. A sensitivity analysis study with respect to EGORSE hyper-parameters is detailed in Section V. Using academic benchmark problems, a comparison with state-of-art HDBO methods shows the high potential of EGORSE both in terms the efficiency and the global computational effort (see Sections VI and VII). Conclusions and perspectives are drawn in Section VIII.

III. Bayesian optimization

A. The general framework

To solve the unconstrained optimization problem (1), the Bayesian optimization (BO) framework [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF][START_REF] Močkus | On Bayesian Methods for Seeking the Extremum[END_REF] builds a surrogate model (using a Gaussian process [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Krige | A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand[END_REF]) of the objective function 𝑓 using a set of 𝑙 sampled points in the design domain Ω, known as the Design of Experiments (DoE). The optimal solution is estimated by iteratively enriching a Gaussian process (GP) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Krige | A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand[END_REF] via a search strategy that balances the exploration of the design space Ω and the minimization of the surrogate model of 𝑓 . Namely, at each iteration of a given BO method, the search strategy requires solving a maximization sub-problem where the objective is referred as the acquisition function [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF][START_REF] Wang | Max-Value Entropy Search for Efficient Bayesian Optimization[END_REF][START_REF] Tran | aphBO-2GP-3B: A Budgeted Asynchronously-Parallel Multi-Acquisition for Known/Unknown Constrained Bayesian Optimization on High-Performing Computing Architecture[END_REF]. The acquisition function being fully defined using the GP, the search strategy is computationally inexpensive and straightforward compared to the original optimization problem [START_REF] Cramer | Problem Formulation for Multidisciplinary Optimization[END_REF] in which the function 𝑓 is assumed to be expensiveto-evaluate. The DoE is updated iteratively using the optimal solutions of the sub-problems. The same process is repeated until a maximum number of evaluations is reached. The main steps of the BO framework, when applied to solve the optimization problem [START_REF] Cramer | Problem Formulation for Multidisciplinary Optimization[END_REF], are summarized in Algorithm 1.

Algorithm 1

The Bayesian optimization framework. input : Objective function, initial DoE, a maximum number of iterations max_nb_it 1: for 𝑙 = 0 to max_nb_it -1 do

2:

Build the surrogate model using a GP 3:

Find 𝒙 (𝑙+1) a solution of the enrichment maximization sub-problem 4:

Evaluate the objective function at 𝒙 (𝑙+1)

5:

Update the DoE 6: end for output : The best point found in the DoE In the next two sections, we give more details on the GP and the enrichment process.

B. Gaussian process

In the context of a BO process, a scalar GP [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Krige | A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand[END_REF] is a surrogate model whose distribution is fully described by a prior mean function, a covariance kernel and a set of sampled points. The global behavior of the GP is depicted by the prior mean function whereas the covariance kernel characterizes the similarities between two sampled points of the design space Ω. Let a non-conditioned scalar GP defined by a prior mean function 𝜇 : R 𝑑 ↦ → R and a covariance kernel 𝑘 : R 𝑑×2 ↦ → R. A conditioned scalar GP of 𝑓 by the DoE of 𝑙 sampled points D (𝑙) = 𝒙 (𝑘 ) , 𝑦 (𝑘 ) 𝑘=0,...,𝑙-1 , where 𝒙 (𝑘 ) ∈ Ω and 𝑦 (𝑘 ) = 𝑓 𝒙 (𝑘 ) , defines a Gaussian distribution N μ(𝑙) , σ(𝑙) 2 for each 𝒙 ∈ Ω. The mean μ(𝑙) : R 𝑑 ↦ → R and standard deviation σ(𝑙) : R 𝑑 ↦ → R are expressed as follows:

μ(𝑙) (𝒙) = 𝜇(𝒙) + 𝒌 (𝑙) (𝒙) ⊤ 𝑲 (𝑙) -1 𝒀 (𝑙) -𝝁 (𝑙) , (2)
σ(𝑙) (𝒙) = 𝑘 (𝒙, 𝒙) -𝒌 (𝑙) (𝒙) ⊤ 𝑲 (𝑙) -1 𝒌 (𝑙) (𝒙) 1 2 , (3)
where 𝝁 (𝑙) = 𝜇 𝒙 (0) , . . . , 𝜇 𝒙 (𝑙-1) ⊤ is the prior mean vector computed on the sampled points of D (𝑙) , 𝒌 (𝑙) = 𝑘 𝒙, 𝒙 (0) , . . . , 𝑘 𝒙, 𝒙 (𝑙-1) ⊤ is the covariance vector between 𝒙 and the sampled points of D (𝑙) , 𝑲 (𝑙) = 𝑘 𝒙 (𝑖) , 𝒙 ( 𝑗 ) 𝑖, 𝑗=0,...,𝑙-1 is the covariance matrix computed on the D (𝑙) , and 𝒀 (𝑙) = 𝑦 (0) , . . . , 𝑦 (𝑙-1) ⊤ is a vector of outputs of 𝑓 . Note that there is a wide range of prior mean functions and covariance kernels [START_REF] Duvenaud | Structure Discovery in Nonparametric Regression through Compositional Kernel Search[END_REF] and their selection is very case dependent. Most of these functions depend on hyper-parameters, denoted by 𝜽 (𝑙) ∈ [R + ] 𝑛 , that need to be estimated to explain the best the DoE of the objective function 𝑓 . To estimate the hyper-parameters 𝜽 (𝑙) of the GP at each iteration, a maximum likelihood estimator is typically used [START_REF] Gelman | Bayesian Data Analysis[END_REF].

C. The enrichment sub-problem

The BO framework combines the information provided by the GP (namely, μ(𝑙) and σ(𝑙) ) to build the enrichment strategy. The latter is guided by the following maximization sub-problem:

max 𝒙∈Ω 𝛼 (𝑙) (𝒙), (4)
where 𝛼 (𝑙) : R 𝑑 ↦ → R is the acquisition function. There are numerous acquisition functions in the literature [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF][START_REF] Wang | Max-Value Entropy Search for Efficient Bayesian Optimization[END_REF], their choice is essential for the enrichment process. The Expected Improvement (EI) [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] acquisition function is the most used in BO. Considering the 𝑙 th iteration of the BO framework, the expression 𝛼 (𝑙) 𝐸 𝐼 depends on the μ(𝑙) and σ(𝑙) . For a given point 𝒙 ∈ Ω, if σ(𝑙) (𝒙) = 0, then 𝛼 (𝑙) 𝐸 𝐼 (𝒙) = 0. Otherwise,

𝛼 (𝑙) 𝐸 𝐼 (𝒙) = 𝑦 (𝑙) 𝑚𝑖𝑛 -μ(𝑙) (𝒙) Φ 𝑦 (𝑙) 𝑚𝑖𝑛 -μ(𝑙) (𝒙) σ (𝑙) (𝒙) + σ (𝑙) (𝒙)𝜙 𝑦 (𝑙) 𝑚𝑖𝑛 -μ(𝑙) (𝒙) σ (𝑙) (𝒙) , (5)
where the functions Φ and 𝜙 are, respectively, the cumulative distribution function and the probability density function of the standard normal distribution. The current minimum is given by 𝑦 (𝑙) 𝑚𝑖𝑛 = min 𝒀 (𝑙) . In this framework, it is possible to tackle problems with non linear constraints [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Priem | Optimisation bayésienne sous contraintes et en grande dimension appliquée à la conception avion avant projet[END_REF] using different mechanisms with different computational costs. However, the classical BO process can not handle high dimensional problem because of the GP model. Indeed, building a conventional GP in high-dimension can be problematic due to the likelihood maximization step used to estimate the hyper-parameters. Moreover, classical GPs tend to miss some information in high dimension as the distance between points increase. In the next section, we propose a solution to help overcome these challenges.

IV. Supervised linear embeddings for BO A. Method description

In this paper, the objective function is supposed to depend only on the effective dimensions 𝑑 𝑒 ≪ 𝑑 where we typically assume that it exists a function 𝑓 A : R 𝑑 𝑒 ↦ → R such as 𝑓 A ( 𝑨𝒙) = 𝑓 (𝒙) with 𝑨 ∈ R 𝑑 𝑒 ×𝑑 , A = {𝒖 = 𝑨𝒙 ; ∀𝒙 ∈ Ω} and Ω = [-1, 1] 𝑑 [START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF][START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF]. The idea is then to perform the optimization procedure in the reduced linear subspace A so that the number of hyper-parameters to estimate and dimension of the design space are reduced to 𝑑 𝑒 instead of 𝑑. This allows to build the inexpensive GPs and will ease the acquisition function optimization. Using a subspace (based on 𝑨) for the optimization requires finding the effective dimension of the reduced design space B ⊂ R 𝑑 𝑒 as well as the backward application 𝛾 : B ↦ → Ω.

The proposed method focuses on the definition of the optimization problem when a linear subspace is used as well as on efficient construction procedure of such embeding subspaces. Most existing HDBO methods rely on random linear subspaces meaning that no information is used to incorporate a priori information from the optimization problem within the embeding space which may slow down the optimization process. In this work, a recursive search, with 𝑇 ∈ N supervised reduction dimension methods, is performed to find supervised linear subspace so that the most important search directions for the exploration of the objective function are included. Using an initial DoE, one can use a Partial Least Squares (PLS) regression [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF] to build such linear embedding prior to the optimization process. Furthermore, the new search design of the optimization problem (within the linear embeding subspace) is a necessary step. Most methods rely on a classic optimization problem formulation that may limit the process performance due to very restricted new design space. Here, once an appropriate linear subspace is found, the optimization problem is turned into a constrained optimization problem to limit the computational cost of the algorithm. The constrained optimization problem can be solved using a classical CBO method [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF][START_REF] Priem | Optimisation bayésienne sous contraintes et en grande dimension appliquée à la conception avion avant projet[END_REF][START_REF] Priem | Upper Trust Bound Feasibility Criterion for Mixed Constrained Bayesian Optimization with Application to Aircraft Design[END_REF].

B. Definition of the reduced search spaces

To define the optimization problem in the low dimensional spaces, transfer matrices from the initial to the low dimensional spaces and optimization domains must be defined. The definition of these transfer matrices relies on a set of 𝑇 ∈ N dimension reduction methods. For each of the 𝑇 dimension reduction methods R (𝑡 ) , the transfer matrix 𝑨 (𝑡 ) ∈ R 𝑑 𝑒 ×𝑑 is build using R (𝑡 ) where 𝑡 ∈ {1, . . . , 𝑇 }. In this way, we propose to use supervised dimension reduction algorithms, like the PLS [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF], to guide the optimization process through highly varying linear subspaces of the objective function. This procedure allows to tackle the issue of Wang et al. [START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF], Binois et al. [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF] by relying on random linear subspaces defined with random Gaussian transfer matrices. In their works, the optimization can be performed in a subspace in which the objective function is not varying, meaning the optimum of the objective function could not be discovered. In the BO framework, the optimization process is usually performed in an hyper-cube. However A (𝑡 ) is not an hypercube. As Ω = [-1, 1] 𝑑 , it is possible to compute B (𝑡 ) ⊂ R 𝑑 𝑒 [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF] the smallest hyper-cube containing all points of A (𝑡 ) such as

B (𝑡 ) = - 𝑑 ∑ 𝑖=1 𝐴 (𝑡 ) 1,𝑖 , 𝑑 ∑ 𝑖=1 𝐴 (𝑡 ) 1,𝑖 × • • • × - 𝑑 ∑ 𝑖=1 𝐴 (𝑡 ) 𝑑 𝑒 ,𝑖 , 𝑑 ∑ 𝑖=1 𝐴 (𝑡 ) 𝑑 𝑒 ,𝑖 . (6)
Performing the optimization process in B (𝑡 ) leads to define a backward application from B (𝑡 ) to Ω to compute the objective function on the desired point.

C. The backward application

We use the backward application introduced by [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF]. Namely, a bijective application 𝛾 (𝑡 ) 𝐵 :

A (𝑡 ) ⊂ R 𝑑 𝑒 ↦ → Ω ⊂ R 𝑑 such that 𝛾 (𝑡 ) 𝐵 (𝒖) = arg min 𝒙∈Ω 𝒙 -𝑨 (𝑡 ) + 𝒖 2 , s.c. 𝑨 (𝑡 ) 𝒙 = 𝒖 , (7)
where 𝑨 (𝑡 ) + = 𝑨 (𝑡 ) ⊤ 𝑨 (𝑡 ) 𝑨 (𝑡 ) ⊤ -1 the pseudo-inverse of 𝑨 (𝑡 ) . This problem requires to solve a quadratic optimization problem. As A (𝑡 ) ⊂ B (𝑡 ) , 𝛾 (𝑡 ) 𝐵 defines an injection from B (𝑡 ) to Ω, meaning some points of B (𝑡 ) do not have any image in Ω. For instance, Figure 1 displays Ω ⊂ R 10 projected in a domain B (𝑡 ) ⊂ R 2 . The A (𝑡 ) domain is in white while the points in the black domain do not have any image in Ω by 𝛾 (𝑡 ) 𝐵 .

2 1 0 1 2 x 1 2 1 0 1 2 x 2
Fig. 1 Illustration of B (𝑡 ) and A (𝑡 ) . In black, the points of B (𝑡 ) without image in Ω by 𝛾 (𝑡 ) 𝐵 ; in white, the points of A (𝑡 ) corresponding to the points of B (𝑡 ) with an image in Ω by 𝛾 (𝑡 ) 𝐵 .

The method developed by [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF] is using the 𝛾 (𝑡 ) 𝐵 even if such backward application is not defined all over B (𝑡 ) by itself. In fact, the function 𝑓 (𝑡 ) (𝒖) = 𝑓 𝛾 (𝑡 ) 𝐵 (𝒖) is only defined on A (𝑡 ) , meaning the optimization problem in the linear subspace is given by: min

𝒖 ∈ A (𝑡 ) 𝑓 (𝑡 ) (𝒖) . (8)
However, standard BO algorithms are only solving optimization problems with an objective function defined on an hyper-cube. Binois et al. [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF] used B (𝑡 ) , see [START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF], the smallest hyper-cube containing A (𝑡 ) . This way, some points reachable by the BO algorithm might not be evaluated on the objective function as they do not have any image in Ω using 𝛾 (𝑡 ) 𝐵 . To fix this issue, an additional modification of the EI acquisition function is thus introduced in [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF]. In fact, over A (𝑡 ) , the acquisition function 𝛼 (𝑙) 𝐸 𝐼 𝑒𝑥𝑡 (𝒖) = 𝛼 (𝑙) 𝐸 𝐼 (𝒖) (5) while on B (𝑡 ) \A (𝑡 ) , the acquisition function 𝛼 (𝑙) 𝐸 𝐼 𝑒𝑥𝑡 (𝒖) = -∥𝒖∥ 2 . As the EI acquisition function is positive (i.e., lower bounded), the optimization of the acquisition function provides a point of A (𝑡 ) . This quadratic problem [START_REF] Gelbart | Bayesian Optimization with Unknown Constraints[END_REF] must be solved at each acquisition function evaluation to know if 𝒖 belongs to A (𝑡 ) . Depending on the number of evaluations of the 𝛼 (𝑙) 𝐸 𝐼 𝑒𝑥𝑡 acquisition function, an optimization process performed with the RREMBO algorithm can be expensive in CPU time. Actually the cost in CPU time grows with the number of design variables 𝑑. Figure 2a shows the objective function of 10 design variables projected in a linear subspace of 2 dimension. The grey area shows that an important part of B (𝑡 ) does not have any image in Ω. In Figure 2b, one can see the grey domain is replaced by negative values increasing towards the center of B (𝑡 ) allowing convergence to points having an image by 𝛾 (𝑡 ) 𝐵 in Ω.

D. The optimization problem

Section IV.C shows that the backward application 𝛾 (𝑡 ) 𝐵 is bijective from A (𝑡 ) to Ω and injective from B (𝑡 ) to Ω. Thus some points of B (𝑡 ) do not have image by 𝛾 (𝑡 ) 𝐵 in Ω. To avoid points of B (𝑡 ) \A (𝑡 ) , we define an optimization problem with a constraint which is feasible on A (𝑡 ) and unfeasible on B (𝑡 ) \A (𝑡 ) . For the value of the constraint when 𝒖 ∉ A (𝑡 ) , we chose the opposite of the 2-norm of 𝒖 . The constraint value is negative and tends to zero when 𝒖 is getting closer to A (𝑡 ) in norm. To define the constraint value in the feasible zone, ones can rely on [START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF] and on Ω = [-1, 1] 𝑑 . The equation shows that points on the edge of B (𝑡 ) , having an image in Ω by 𝛾 (𝑡 ) 𝐵 , are in the corners of Ω. The corners of a domain are considered to be the points whose components are -1 or 1. However, the corners of the domain are the farthest points from the center 𝑥 𝑐 = 0 ∈ R 𝑑 and have the same norm. All the images from points of B (𝑡 ) by 𝛾 (𝑡 ) 𝐵 have a lower norm than the corners of Ω. Hence the constraint of the optimization problem is defined as the difference between the norm of the corners Ω and the norm of 𝛾 (𝑡 ) 𝐵 (𝒖) when 𝒖 ∈ A (𝑡 ) . Thus, the constraint value tends to zero when 𝒖 is getting closer to B (𝑡 ) \A (𝑡 ) . Eventually, the constraint function is normalized to provide values in [-1, 1]. This normalization on the bounds of the constraint function balances the importance of the feasible and unfeasible domain in the optimization process.

To summarize, the constraint is given by 𝑔 (𝑡 ) (𝒖) ≥ 0 where:

𝑔 (𝑡 ) (𝒖) =          1 - 𝛾 (𝑡 ) 𝐵 (𝒖) 2 2 𝑑 if 𝒖 ∈ A (𝑡 ) -𝒖 𝑨 (𝑡 ) 2 2 otherwise (9) with 𝑢 𝑨 (𝑡 ) 𝑖 = 𝑢 𝑖 / 𝑑 𝑗=1 𝐴 (𝑡 )
𝑖, 𝑗 . The terms 𝑢 𝑨 (𝑡 )

𝑖

and 𝑢 𝑖 are respectively the 𝑖 th components of the 𝒖 𝑨 (𝑡 ) and 𝒖 vectors. One remarks that 𝑓 (𝑡 ) (𝒖) = 𝑓 𝛾 (𝑡 ) 𝐵 (𝒖) function is not defined all over B (𝑡 ) . That is why one seeks to give a value to 𝑓 (𝑡 ) on B (𝑡 ) \A (𝑡 ) to obtain a function not only defined on B (𝑡 ) . The extension of 𝑓 (𝑡 ) on B (𝑡 ) \A (𝑡 ) is given by 𝑓 (𝑡 ) (𝒖) = 𝑓 𝛾 (𝑡 ) 𝑊 (𝒖) with

𝛾 (𝑡 ) 𝑊 (𝒖) ∈ arg min 𝒙∈Ω ∥𝒙 -𝑨 (𝑡 ) + 𝒖∥ 2 . ( 10
)
Indeed, the 𝛾 (𝑡 ) 𝑊 application exists for all points 𝒖 ∈ B (𝑡 ) although it can provide the same 𝒙 ∈ Ω for different 𝒖 ∈ B (𝑡 ) . These points of Ω are moreover reachable with points of A (𝑡 ) using the 𝛾 (𝑡 ) 𝐵 backward application. Thus, one of the objective function minima is necessarily in A (𝑡 ) . Eventually, the objective function 𝑓 (𝑡 ) is given on B (𝑡 ) by

𝑓 (𝑡 ) (𝒖) =        𝑓 𝛾 (𝑡 ) 𝐵 (𝒖) if 𝒖 ∈ A (𝑡 )
𝑓 𝛾 (𝑡 ) 𝑊 (𝒖) otherwise [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF] In fact, one solves the following constrained optimization problem in the B (𝑡 ) hypercube:

min 𝒖 ∈ B (𝑡 ) 𝑓 (𝑡 ) (𝒖) s.c. 𝑔 (𝑡 ) (𝒖) ≥ 0 . ( 12
)
where 𝑓 (𝑡 ) (𝒖) and 𝑔 (𝑡 ) (𝒖) are respectively defined by [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF] and [START_REF] Diouane | TREGO: a trust-region framework for efficient global optimization[END_REF]. To solve this problem, a standard CBO algorithm [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Priem | Optimisation bayésienne sous contraintes et en grande dimension appliquée à la conception avion avant projet[END_REF], like SEGOMOE [START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF][START_REF] Priem | Upper Trust Bound Feasibility Criterion for Mixed Constrained Bayesian Optimization with Application to Aircraft Design[END_REF] can be used. Figure 3a shows the optimization problem [START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF] for the 10 design variable function given in Figure 2a and projected in a 2 dimensional linear subspace. One sees that the grey unfeasible area corresponds to the grey area of Figure 2a. Figure 3b shows the optimization sub-problem of the SEGOMOE algorithm using the EI acquisition function where the green and red squares are points of the initial DoE. The predicted unfeasible zone, in grey in Figure 3b, contains the unfeasible points of the initial DoE. Recall, with SEGOMOE, only the mean of the GP modeling the constraint is used to defined the feasible zones. On the contrary of RREMBO no quadratic problem solving is needed to solve the optimization sub-problem. In fact, the introduced process avoids the quadratic problem solving in the optimization sub-problem by defining functions existing all over B (𝑡 ) , which is not the case of RREMBO (see Section IV.C). The quadratic problem is only solved when the objective function is called at each iteration of the SEGOMOE algorithm. Thus, the introduced EGORSE algorithm should be faster in CPU time than RREMBO.

E. Adaptive learning of the linear subspace

To ease the convergence process, we use a linear subspace discovered by supervised dimension reduction methods, i.e. taking points of the considered domain and the associated function values as inputs. In our case, linear dimension reduction methods are used to provide the transfer matrix 𝑨 (𝑡 ) ∈ R 𝑑×𝑑 𝑒 from Ω to A (𝑡 ) . For instance, the PLS [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF] or the Marginal Gaussian Process (MGP) [START_REF] Garnett | Active Learning of Linear Embeddings for Gaussian Processes[END_REF] methods can be considered to find the transfer matrix. In order to find the minimum of the objective function, we solve the optimization problem [START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF] with a CBO algorithm [START_REF] Frazier | A Tutorial on Bayesian Optimization[END_REF][START_REF] Shahriari | Taking the Human Out of the Loop: A Review of Bayesian Optimization[END_REF][START_REF] Priem | Optimisation bayésienne sous contraintes et en grande dimension appliquée à la conception avion avant projet[END_REF] in B (𝑡 ) generated by the 𝑨 (𝑡 ) matrix with a maximum of max_nb_it_sub iterations. However, in larger space, the linear subspace approximation can lack of accuracy if the number of points used by the dimension reduction methods is not sufficient. To improve the generated subspace, the previous process is iterated using all the evaluated points during the previous iterations. Note that points outside of the subspace A (𝑡 ) are added to provide information which is non biased by the subspace selection. In that way, points coming from an optimization performed in a subspace defined by other dimension reduction methods can be used. For instance, a transfer matrix defined by a random Gaussian distribution is chosen. The convergence properties, given by Binois et al. [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF], are thus preserved. Finally, we generalize this approach by using several dimension reduction methods, which can be unsupervised (like random Gaussian transfer matrices or hash tables [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF][START_REF] Nayebi | A Framework for Bayesian Optimization in Embedded Subspaces[END_REF]), or supervised (like PLS or MGP [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF][START_REF] Garnett | Active Learning of Linear Embeddings for Gaussian Processes[END_REF]), to consider their advantages. The overall process of EGORSE is finally given by Algorithm 2 whereas the flow chart of the method is described by the eXtended Design Structure Matrix (XDSM) [START_REF] Lambe | Extensions to the Design Structure Matrix for the Description of Multidisciplinary Design, Analysis, and Optimization Processes[END_REF] 𝑓 , a maximum number of iterations max_nb_it, a maximum number of iterations by subspace max_nb_it_sub, a number of active directions 𝑑 𝑒 , a list R = R (1) , . . . , R (𝑇 ) of 𝑇 ∈ N + supervised or unsupervised dimension reduction methods 1: for 𝑖 = 0 to max_nb_it -1 do 2:

for 𝑡 = 1 to 𝑇 do Build 𝑓 (𝑡 ) (see [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF]

) 6: Build 𝑔 (𝑡 ) (see (9)) 7:
Solve the optimization problem min

𝒖 ∈ B (𝑡 )
𝑓 (𝑡 ) (𝒖) s.c. 𝑔 (𝑡 ) (𝒖) ≥ 0 with a CBO algorithm using max_nb_it_sub maximum iterations.

D (𝑖+1) 𝑓 = D (𝑖)
𝑓 ∪ {Points already evaluated on 𝑓 at iteration 𝑖} 10: end for output : The best point regarding the value of 𝑓 in D (max_nb_it) 𝑓 contributions:

• The restriction of the number of quadratic problems solved with a new formulation of the subspace optimization problem. • The use of supervised dimension reduction methods promoting the search in highly varying directions of the objective function in Ω. • The adaptive learning of the linear subspace using all the evaluated points.

n it , D (0) f R (t) f ,d e Ω f ,Ω n sub D (nit) f 0,7→1: Loop 1 : D (i) f 1: Build the matrix A (t) 2 : A (t) 3 : A (t)
2: Compute the bounds of the research domain B (t) 3 :

B (t) 4 : B (t)
3: Build the objective function f (t) and the constraint g (t) 5 : f (t) , g (t) 7 :

D (i+1) f 4,6→5: CBO 5 : u (l) 6 : u (l+1) 5: Evaluate f (t)
and g (t)

Fig. 4 An XDSM diagram of the EGORSE framework.

V. A sensitivity analysis over EGORSE hyper-parameters

A. Implementation details EGORSE method is implemented with Python 3.8. A CBO process (see Section IV.D) is performed at each iteration of the EGORSE algorithm using SEGO [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF] (from the SEGOMOE toolbox [START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF]). The SEGOMOE toolbox uses the SMT [START_REF] Bouhlel | A Python Surrogate Modeling Framework with Derivatives[END_REF], a Python package used to build the GP model. In the CBO process, the EI acquisition function is optimized in two steps. First, a good starting point is found by solving the sub-optimization problem (4) with the ISRES [START_REF] Runarsson | Search biases in constrained evolutionary optimization[END_REF] from the NLOPT [START_REF] Johnson | The NLopt nonlinear-optimization package[END_REF] Python toolbox. ISRES is an evolutionary optimization algorithm able to solve multi-modal optimization problems with equality and inequality constraints. The algorithm explores the domain to find an optimal area maximizing the acquisition function and respecting the feasibility criteria. However, such algorithm needs many function evaluations to converge. To limit this number of evaluations, the solution provided by ISRES is refined with a gradient based optimization algorithm as the analytical derivatives of the acquisition function and the feasibility criteria are easily available from the GP approximations and provided as outputs from SMT [START_REF] Bouhlel | A Python Surrogate Modeling Framework with Derivatives[END_REF]. The gradient based algorithm used is SNOPT [START_REF] Gill | SNOPT: An SQP algorithm for large-scale constrained optimization[END_REF] from the PyOptSparse [START_REF] Perez | pyOpt: a Python-based object-oriented framework for nonlinear constrained optimization[END_REF] Python toolbox whose initial guess is given by ISRES. To solve the quadratic problem [START_REF] Gelbart | Bayesian Optimization with Unknown Constraints[END_REF], the CVXOPT [START_REF] Andersen | CVXOPT: Python software for convex optimization[END_REF] Python toolbox is chosen. A point 𝒖 ∈ B (𝑡 ) 𝐵 is considered to belong to A (𝑡 ) if the CVXOPT optimization status 'optimal' is reached meaning that the quadratic problem has a solution.

B. On the setting of EGORSE hyper-parameters

EGORSE is controlled by a set of hyper-parameters. To select the value of these hyper-parameters, a parametric study is performed on two optimization problems. The hyper-parameters considered in this study are the following.

• The number of points in the initial DoE. It impacts the supervised dimension reduction methods. On the overall EGORSE versions, three sizes of initial DoE are tested: 5 points, 𝑑 points and 2𝑑 points where 𝑑 is the number of design variables.

• The supervised dimension reduction method. It changes the behavior of the algorithm favoring different directions of the domain. The PLS [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF] and MGP [START_REF] Gardner | Discovering and Exploiting Additive Structure for Bayesian Optimization[END_REF] methods are considered in this study. The PLS method is coming from the scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] Python toolbox. The MGP method is implemented within SMT [START_REF] Bouhlel | A Python Surrogate Modeling Framework with Derivatives[END_REF]. • The unsupervised dimension reduction method relies on a random Gaussian transfer matrix like in Binois et al.

[13] work or on hash table like in Nayebi et al. [START_REF] Nayebi | A Framework for Bayesian Optimization in Embedded Subspaces[END_REF] work. In what comes next, the following 6 possible variants of EGORSE will be tested and compared:

1) EGORSE Gaussian: it uses Gaussian random transfer matrix.

2) EGORSE Hash: it uses random matrix defined by Hash table.

3) EGORSE PLS: it uses transfer matrix defined by the PLS method. 4) EGORSE PLS+Gaussian: it uses Gaussian random transfer matrix and transfer matrix defined by the PLS. 5) EGORSE MGP: it uses transfer matrix defined by the MGP method. 6) EGORSE MGP+Gaussian: it uses Gaussian random transfer matrix and transfer matrix defined by the MGP.

C. Study on two analytical problems

Definition of the two problems

The considered class of problems is an adjustment of the Modified Branin (MB) problem [START_REF] Parr | Infill Sampling Criteria for Surrogate-Based Optimization with Constraint Handling[END_REF] whose number of design variables is artificially increased. This problem is commonly used in the literature [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF][START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF][START_REF] Nayebi | A Framework for Bayesian Optimization in Embedded Subspaces[END_REF] and it is defined as follows: min

𝒖 ∈Ω 1 𝑓 1 (𝒖), (13)
where [START_REF] Bouhlel | Efficient Global Optimization for High-Dimensional Constrained Problems by Using the Kriging Models Combined with the Partial Least Squares Method[END_REF] and

Ω 1 = [-5, 10] × [0,
𝑓 1 (𝒖) =       𝑢 2 - 5.1𝑢 2 1 4𝜋 2 + 5𝑢 1 𝜋 -6 2 + 10 - 10 8𝜋 cos (𝑢 1 ) + 1       + 5𝑢 1 + 25 15 . ( 14
)
The modified version of the Branin problem is selected because it count three local minima including a global one. The value of the global optimum is about MB_𝑑 𝑚𝑖𝑛 = 1.1. Furthermore, the problem is normalized to have 𝒖 ∈ [-1, 1] 2 .

To artificially increase the number of design variables, a random matrix 𝑨 𝑑 ∈ R 2×𝑑 is generated such that for all 𝒙 ∈ [-1, 1] 𝑑 , 𝑨 𝑑 𝒙 = 𝒖 belongs to [-1, 1] 2 . An objective function MB_𝑑, where 𝑑 is the number of design variables, is defined such that MB_𝑑 (𝒙) = 𝑓 1 ( 𝑨 𝑑 𝒙). Eventually, we solve the following optimization problem:

min 𝑥 ∈ [ -1,1] 𝑑 MB_𝑑 (𝒙) = min 𝑥 ∈ [-1,1] 𝑑 𝑓 1 ( 𝑨 𝑑 𝒙). (15)
In the following numerical experiments are conducted on two functions of respective dimension 10 and 100. These two test functions are denoted MB_10 and MB_100.

Convergence plots

To study the EGORSE hyper-parameter impact, convergence plots are obtained for both problems. Ten independent optimizations are thus performed on each of the problems, for each of the EGORSE versions, using 10 initial DoE. The number of iterations is imposed to 10 and the number of evaluations by sub-space optimization is set to 20𝑑 𝑒 and 𝑑 𝑒 = 2. For EGORSE PLS, EGORSE MGP, EGORSE Hash and EGORSE Gaussian, the number of iterations is doubled to keep a fixed number of evaluations, i.e. 800 evaluations per run. Indeed, at each iteration, the problem is evaluated only 20𝑑 𝑒 for EGORSE composed of a unique dimension reduction method against 2 × 20𝑑 𝑒 for EGORSE composed of two dimension reduction methods. All of the EGORSE versions are then compared displaying the evolution of the mean and standard deviation on the 10 optimizations of the lower value of the objective function evaluated with respect to the number of evaluations. For readability, the standard deviation is displayed with a reduction factor of four.

Results analysis

In this section, the convergence robustness and speed of the six EGORSE versions are analyzed with convergence plots (see Section V.C.2). The convergence plots of the considered EGORSE versions for the MB_10 problems are displayed in Fig. 5. In Fig. 5a, one can see that the EGORSE Gaussian algorithm offers the best performance in term of convergence speed and robustness for an initial DoE of 5 points. However, Figures 5b and5c show that EGORSE PLS+Gaussian outperforms EGORSE Gaussian in term of convergence speed and robustness for larger initial DoE. Figure 5d finally compares these different versions and they are hardly distinguishable for the MB_10 problem. Thus, all the six EGORSE algorithms seem to provide equivalent performance for low dimensional problems. Figure 6 much more distinguishable on Fig. 6. In fact, Figures 6a, 6b and6c show that EGORSE PLS+Gaussian provides the best convergence speed and robustness trade-off for all the initial DoE sizes tested. To select the best number of initial DoE points, Fig. 6d displays the convergence plots of EGORSE PLS+Gaussian for the three tested number of points in the initial DoE. One can easily see that EGORSE PLS+Gaussian with an initial DoE of 𝑑 points provides the best performance in term of convergence speed and robustness. in terms of evaluations number, the use of an initial DoE of 𝑑 points allows the algorithm to explore the domain in an interesting direction more quickly than a 2𝑑 points initial DoE. On the contrary, using an initial DoE of 5 points forces the algorithm to seek for the best direction for a long time.

To conclude, choosing the PLS and Gaussian dimension reduction methods with an initial DoE of 𝑑 points seems the most suitable to obtain the best performance of EGORSE. Nevertheless, EGORSE capabilities have to be compared with HDBO algorithms to validate its usefulness as it is proposed in the following.

VI. Comparison with state-of-the-art HDBO methods

A. BO algorithms and setup details

EGORSE is now compared to the following state-of-the-art algorithms:

• TuRBO [10]: a HDBO algorithm using confidence regions to favor the exploitation of the DoE data. Tests are performed with the TuRBO A [START_REF] Eriksson | Scalable Global Optimization via Local Bayesian Optimization[END_REF] Python toolbox. • EGO-KPLS [START_REF] Bouhlel | Efficient Global Optimization for High-Dimensional Constrained Problems by Using the Kriging Models Combined with the Partial Least Squares Method[END_REF]: an HDBO method relying on the reduction of the number of GP hyper-parameters. This allows to speed up the GP building. The SEGOMOE [START_REF] Bartoli | Adaptive Modeling Strategy for Constrained Global Optimization with Application to Aerodynamic Wing Design[END_REF] Python toolbox is used. All the hyper-parameters of this algorithm are the default ones. The number of principal components for the KPLS model is set to two. • RREMBO [START_REF] Binois | On the Choice of the Low-Dimensional Domain for Global Optimization via Random Embeddings[END_REF]: a HDBO method using the random Gaussian transfer matrix to reduce the number of dimensions of the optimization problem. RREMBO B implementation of this algorithm is used. The parameters are also set by default.

• HESBO [START_REF] Nayebi | A Framework for Bayesian Optimization in Embedded Subspaces[END_REF]: a HDBO algorithm using Hash tables to generate the transfer matrix. We use the HESBO C Python toolbox with the default parameters. For EGORSE, the version showing the best performance in term of convergence speed and robustness in Section V.C.3 is selected, i.e. EGORSE PLS + Gaussian with an initial DoE of 𝑑 points. To achieve this comparison, 10 optimizations for each problem and for each studied method are completed to analyze the statistical behavior of these BO algorithms. Because of the different strategies implemented in the previously introduced algorithms, a specific test plan must be adopted for each of them.

• EGORSE: The test plan of Section V.C.2 is implemented.

• TuRBO: Five trust regions are used with a maximum of 800 evaluations of the objective function. Note that it is not possible to provide the same initial DoE used for EGORSE in TuRBO. The number of points generated at the beginning of the algorithm is thus imposed to 𝑑. The EI acquisition function is chosen. • EGO-KPLS: The optimization is performed with a maximum evaluation number of 800, with the initial DoE used in EGORSE and with the EI acquisition function.

• RREMBO & HESBO: 20 optimizations of 20𝑑 𝑒 evaluations are performed for each of the initial DoE for RREMBO and HESBO. These 20 optimizations are then concatenated and considered as a unique optimization process. The number of effective directions are imposed to 𝑑 𝑒 = 2 and the acquisition function is EI. RREMBO and HESBO are equivalent to EGORSE without supervised dimension reduction method.

B. Results analysis

In this section, a comparison between EGORSE and the four studied algorithms is performed in term of robustness, convergence speed both in CPU time and in number of iterations. Fig. 7 provides the iteration convergence plots, as introduced in Section V.C.2, and the time convergence plots drawing the evolution of the means of the best discovered function values against the CPU time.

Fig. 7a shows that TuRBO and EGO-KPLS are converging the fastest and with a low standard deviation. Moreover, the convergence plots of EGORSE, RREMBO and HESBO are hardly distinguishable. Fig. 7b displays that EGO-KPLS converges the fastest to the lowest values with a low standard deviation. TuRBO is also providing good performance even if it converges slower than EGO-KPLS. Regarding the three methods using dimension reduction procedure, EGORSE is converging to the lowest value with a relatively low standard deviation. The good performance of EGO-KPLS and TuRBO is certainly due to the ability of these algorithms to search all over Ω, which is not the case for other methods. However, when the dimension of Ω increases, the ability to search all over Ω becomes a drawback. In fact, a complete search in Ω is intractable in time is this case.

Figures 7c and7d depict the convergence CPU time necessary to obtain the regarded value. First, the RREMBO, TuRBO and EGO-KPLS complete the optimization procedure in more than 8 hours on the MB_100 problem against an hour on the MB_10 problem. This suggests that RREMBO, TuRBO and EGO-KPLS are intractable in time for larger problems. Then, one can easily see that EGORSE is converging the fastest in CPU time than the other algorithms on the MB_100 problem. In fact, the computation time needed to find the enrichment point is much lower than the one for TuRBO, EGO-KPLS and RREMBO. This was sought in the definition of the enrichment sub-problem introduced in Section IV.D. Finally, EGORSE is converging to a lower value than HESBO in a similar amount of time. Thus, EGORSE seems more interesting to solve HDBO problems than the studied algorithms. Note that only HESBO and EGORSE are able to perform an optimization procedure on HDBO problems.

To conclude this section, several sets of hyper-parameters of EGORSE have been tested on two problems of dimension 10 and 100. The analysis of the obtained results has shown that EGORSE PLS+Gaussian with an initial DoE of 𝑑 points is performing the best. A comparison of EGORSE with HDBO algorithms has also been carried out. It has pointed out that TuRBO, EGO-KPLS and RREMBO are intractable in time for HDBO problems. Furthermore, EGORSE has appeared to be the most suitable to solve HDBO problem efficiently.

VII. Evaluation of EGORSE on a high dimensional planning optimization

In this section, the EGORSE algorithm is evaluated to find an optimal path planning problem using 600 design optimization variables.

A. Problem definition and implementation details

The Rover_600 path planning problems relies on the same idea than MB_d problems except that the objective function is a adjustment of the Rover_60 [START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF] problem. It consists on a robot routing from a starting point 𝑥 𝑠𝑡 𝑎𝑟𝑡 to a goal point 𝑥 𝑔𝑜𝑎𝑙 in a forest. The robot trajectory is a spline defined by 30 control points. These points, including the starting and the goal ones, are the design variables of the optimization problem that belong to Ω = [0, 1] 60 . The objective function is minimal when the robot follows the shortest trajectory without meeting a tree. The minimum of the function is 𝑓 𝑚𝑖𝑛 = -5. Fig. 8 gives an example of trajectory. To increase the number of design variables, the problem is normalized in Ω = [-1, 1] 60 , a random matrix 𝑨 𝑑 ∈ R 60×600 is generated such that all 𝒙 ∈ [-1, 1] 𝑑 , 𝑨 𝑑 𝒙 = 𝒖 ∈ [-1, 1] 60 . An objective function Rover_600, where 𝑑 = 600 is the number of design variables, is defined such that Rover_600(𝒙) = Rover_60( 𝑨 𝑑 𝒙). Eventually, we solve the following optimization problem:

min 𝑥 ∈ [ -1,1] 6 00 Rover_600(𝒙) = min 𝑥 ∈ [-1,1] 60 Rover_60( 𝑨 𝑑 𝒙). (16)
We note that, in our tests, TuRBO, EGO-KPLS and RREMBO algorithms are not included anymore in this comparison as they are intractable in time for optimization problems with more than a hundred design variables (see Section VI.B). So, only EGORSE PLS+Gaussian is compared to HESBO. The same test plan as in Section VI.A is used with 200 optimizations (instead of 20) of 20𝑑 𝑒 evaluations.

B. Result analysis

Time and iteration convergence plots of EGORSE and HESBO are depicted in Fig. 9. Fig. 9b clearly shows that EGORSE is converging fast to the lowest objective value with a very low standard deviation. However, the obtained objective value is larger than the known optimal one (i.e. 𝑓 𝑚𝑖𝑛 = -5). This is due to two main reasons that may be addressed in further research:

• The number of effective directions used in EGORSE (i.e. 𝑑 𝑒 = 2) is much lower than the actual number of effective directions (i.e. 𝑑 𝑒 = 60). The effective search space is not covering the space in which Rover_600 is varying. • The dimension reduction method PLS is global. The local variations of the function, in which the global optimum can be located, are thus deleted. Even if this problem is tackled by searching in randomly generated subspace, EGORSE cannot provide better results. Figure 9a shows that HESBO is performing the optimization procedure faster than EGORSE. In fact, HESBO does not solve any quadratic problem at each iteration. However, one can see that the time difference is not significant.

VIII. Conclusion

This paper introduces EGORSE, a high-dimensional efficient global optimization using both random and supervised embeddings to tackle expensive to compute black-box optimization problems. EGORSE shows a high potential to tackle the two main drawbacks of most existing HDBO methods with a very competitive computational effort and with an accurate definition of the new bounds related to the reduced optimization problem. A parametric study on the hyper-parameters of EGORSE has shown that combining both PLS and the random Gaussian reduction methods provides the best results. On very large-scale optimization problem, i.e. d=600, EGORSE gives a lower minimal value than HESBO in an equivalent computational time. Thus EGORSE outperforms by far all the state-of-the-art HDBO solvers. Nonetheless, EGORSE is still unable to find the global optimum in many cases because the effective dimension of the problem is much larger than the effective dimension used in the algorithm. The effective directions employed in the optimization process seem not sufficient to cover the complete variability of the function. Many perspectives could be envisaged to follow this work:

• an automatic choice of the number of effective dimensions, inspired by [START_REF] Saves | Bayesian optimization for mixed variables using an adaptive dimension reduction process: applications to aircraft design[END_REF],

• an extension of EGORSE for constrained expensive-to-compute problems • an extension of EGORSE to improve its convergence, for instance with trust region strategies [START_REF] Diouane | TREGO: a trust-region framework for efficient global optimization[END_REF].

IX. Appendix: review on PLS and MGP methods

In this section, a recall on PLS and MGP is provided. More details can be found in [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF][START_REF] Garnett | Active Learning of Linear Embeddings for Gaussian Processes[END_REF].

A. Partial least squares (PLS) based method:

The PLS [START_REF] Helland | On the Structure of Partial Least Squares Regression[END_REF] method searches for the most important 𝑑 𝑒 orthonormal directions 𝒂 (𝑖) ∈ R 𝑑 , 𝑖 ∈ {1, . . . , 𝑑 𝑒 } of Ω in regards of the influence of the inputs 𝑿 (𝑙) (0) = 𝒙 (0) ⊤ , . . . , 𝒙 (𝑙) ⊤ ⊤ ∈ R 𝑑×𝑙 on the outputs 𝒀 (𝑙) 𝑠, (0) = 𝒀 (𝑙) 𝑠 ∈ R 𝑙 . These directions are recursively defined as follows:

𝑿 (𝑙) (𝑖+1) = 𝑿 (𝑙) (𝑖) -𝒕 (𝑖) 𝒑 (𝑖) , 𝒀 (𝑙) 𝑠, (𝑖+1) = 𝒀 (𝑙) 𝑠, (𝑖) -𝑐 (𝑖) 𝒕 (𝑖) ,

where 𝑿 (𝑙) (𝑖+1) ∈ R 𝑑×𝑙 and 𝒀 (𝑙) 𝑠, (𝑖+1) ∈ R 𝑑 are the residuals of the projection of 𝑿 (𝑙) (𝑖) and 𝒀 (𝑙) 𝑠, (𝑖) on the i th principal component 𝒕 (𝑖) ∈ R 𝑑 . Finally, 𝒑 (𝑖) ∈ R 𝑑 and 𝑐 (𝑖) ∈ R are respectively the regression coefficients of 𝑿 (𝑙) (𝑖) and 𝒀 (𝑙) 𝑠, (𝑖) on the i th principal component 𝒕 (𝑖) for 𝑖 ∈ {1, . . . , 𝑑 𝑒 }. In fact, the square of the covariance between 𝒂 ′ (𝑖) and 𝑿 (𝑙) (𝑖) ⊤ 𝒀 (𝑙) 𝑠, (𝑖) is recursively maximized.

B. Marginal Gaussian process (MGP) based method:

With the MGP [START_REF] Garnett | Active Learning of Linear Embeddings for Gaussian Processes[END_REF], the vectored matrix vect( 𝑨) ∈ R 𝑑•𝑑 𝑒 is considered to be a realization of a Gaussian distribution P( 𝐴) = N vect 𝑨 𝑝 , 𝚺 𝑝 where vect 𝑨 𝑝 and 𝚺 𝑝 are respectively the prior mean and the covariance matrix. The density function of P( 𝐴) is noted 𝑝( 𝑨). Then, the posterior probability law of vect( 𝑨) with respect to D (𝑙) 𝑠 , P 𝑨|D (𝑙) 𝑠 , is estimated with the following Laplace approximation:

P 𝑨|D (𝑙) 𝑠 ≈ N vect( Â), Σ , (21)
𝑝 𝑨|D (𝑙)

𝑠 ≈ 𝑝( 𝑨) • L 𝒀 (𝑙) 𝑠 , 𝑨 , (22)
where vect  is a local maximum of 𝑝 𝑨|D where ∇ 2 is the Hessian operator with respect to vect( 𝑨). Although the derivatives and Hessian of the log prior are trivial, the log likelihood ones are not and are provided in Garnett et al. [START_REF] Garnett | Active Learning of Linear Embeddings for Gaussian Processes[END_REF] work. The remaining details on the MGP construction are not necessary to understand the definition of the transfer matrix. In this way, they are not introduced but can be found in Garnett et al. [START_REF] Garnett | Active Learning of Linear Embeddings for Gaussian Processes[END_REF]. In fact, the transfer matrix used when referring to the MGP dimension reduction method is given by Â.

Ω

  = design space 𝑑 = dimension of the design space 𝑓 = objective function 𝒙 = vector of design variables 𝜇 = prior mean function 𝑘 = covariance kernel D (𝑙) = design of experiments of 𝑙 sampled points 𝑦 = output of the objective function N

  (a) Projected objective function. (b) 𝛼 (𝑙) 𝐸 𝐼𝑒𝑥𝑡 .

Fig. 2

 2 Fig. 2 Projection of an objective function of 10 design variables into a linear subspace of 2 dimensions and the corresponding 𝛼 (𝑙) 𝐸 𝐼 𝑒𝑥𝑡 acquisition function. The grey area is the unfeasible domain, the green squares are DoE points of A (𝑡 ) , the red squares are DoE points of B (𝑡 ) \A (𝑡 ) and the green star is a solution of the optimization sub-problem.

  (a) Optimization problem[START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF].(b) SEGOMOE optimization sub-problem associated to problem[START_REF] Wang | Batched High-Dimensional Bayesian Optimization via Structural Kernel Learning[END_REF].

Fig. 3

 3 Fig. 3 An optimization problem (12) in a 2 dimensional linear subspace and the associated SEGOMOE optimization sub-problem. The grey area is the unfeasible domain, the green squares are DoE points of A (𝑡 ) , the red squares are DoE points of B (𝑡 ) \A (𝑡 ) and the green star is a solution of the optimization sub-problem.

diagram of Figure 4 .

 4 To conclude, the EGORSE algorithm includes three main Algorithm 2 The EGORSE process. input : an objective function 𝑓 , an initial DoE D (0)

3 : 4 :

 34 Build 𝑨 (𝑡 ) ∈ R 𝑑 𝑒 ×𝑑 with R (𝑡 ) using or not D (𝑖) 𝑓 Define B(𝑡 ) (see (6)) 5:

  best valid values of the objective EGORSE Gaussian (DOE:5 pts) EGORSE PLS + Gaussian (DOE:D pts) EGORSE PLS + Gaussian (DOE:2D pts) (d) Best versions for each DoE size.

Fig. 5

 5 Fig. 5 Convergence plots of 6 versions of EGORSE applied on the MB_10 problem. The grey vertical line shows the size of the initial DoE.

  DoE of 𝑑 points.

  best valid values of the objective EGORSE PLS + Gaussian (DOE:5 pts) EGORSE PLS + Gaussian (DOE:D pts) EGORSE PLS + Gaussian (DOE:2D pts) (d) Best versions for each DoE size.

Fig. 6

 6 Fig. 6 Convergence plots of 6 versions of EGORSE applied on the MB_100 problem. The grey vertical line shows the size of the initial DoE.

A

  https://github.com/uber-research/TuRBO B https://github.com/mbinois/RRembo C https://github.com/aminnayebi/HesBO MB_100 iteration convergence plot.

Fig. 7

 7 Fig. 7 Iteration and time convergence plots for 5 HDBO algorithms on the MB_10 and MB_100 problems. The grey vertical line shows the size of the initial DoE.

Fig. 8

 8 Fig. 8 Example of a robot trajectory in a forest.

Fig. 9

 9 Fig. 9 CPU time convergence plots of HESBO, RREMBO, TuRBO, EGO-KPSL on the Rover_600 problem with an initial DoE of 𝑑 points. The grey vertical line shows the size of the initial DoE.

  𝒂 ′(𝑖) ∈ arg max𝒂 ′ ∈R 𝑑 𝒂 ′⊤ 𝑿 (𝑙) (𝑖) ⊤ 𝒀 (𝑙) 𝑠, (𝑖) 𝒀 (𝑙) 𝑠, (𝑖) ⊤ 𝑿 (𝑙) (𝑖) 𝒂 ′ , 𝒂 ′⊤ 𝒂 ′ = 1 ,(17)𝒕 (𝑖) = 𝑿 (𝑙) (𝑖) 𝒂 ′ (𝑖) , 𝒑 (𝑖) = 𝑿 (𝑙) (𝑖) ⊤ 𝒕 (𝑖) , 𝑐 (𝑖) = 𝒀 (𝑙)𝑠, (𝑖) ⊤ 𝒕 (𝑖)

𝑨 ′ = [𝒂 ′ ( 1 )

 1 , . . . , 𝒂 ′ (𝑑 𝑒 ) ], 𝑷 = [ 𝒑 ′ (1) , . . . , 𝒑 ′ (𝑑 𝑒 ) ], 𝑨 = 𝑨 ′ (𝑷 ⊤ 𝑨 ′ ) -1 ,

  function of P 𝑨|D (𝑙) 𝑠 and L the likelihood. The covariance matrix Σ is estimated by the inverse of the logarithm of 𝑝 𝑨|D (𝑙) 𝑠 Hessian matrix evaluated at vect  Σ-1 = -∇ 2 log 𝑝 𝑨|D (𝑙) 𝑠 𝑨=  ≈ -∇ 2 log L 𝒀 (𝑙) 𝑠 , 𝑨 𝑨=  -∇ 2 log 𝑝 ( 𝑨) 𝑨=  (23)

  displays the convergence plots of the six EGORSE versions applied to the MB_100 problem. One can see that the different EGORSE versions are

		EGORSE Hach EGORSE MGP EGORSE MGP + Gaussian	EGORSE Gaussian EGORSE PLS EGORSE PLS + Gaussian
	30
	15 20 25 Mean of the best valid values of the objective
	10
	0	100	200	300 Number of Evaluations 400 500	600	700	800

Acknowledgments

This work is part of the activities of ONERA -ISAE -ENAC joint research group. The authors would like to express their special thanks to Nina Moello for her test on the implementation of the log-likelihood gradient and Hessian in the python toolbox SMT and for her implementation of the HESBO algorithm in our computer environment.


R Ravaud

G Lemarquand

V Lemarquand

Force and Stiffness of Passive Magnetic Bearings Using Permanent Magnets. Part 2 : Radial Magnetization

Keywords: Magnetic forces, analytical calculation, ring permanent magnet, magnetic bearing

This paper deals with the calculation of the force and the stiffness between two ring permanent magnets whose polarization is radial. Such a configuration corresponds to a passive magnetic bearing. The magnetic force exerted between ring permanent magnets is determined by using the coulombian model. The expressions obtained are semi-analytical and we show that it is not possible to find an exact analytical expression of the force between two ring permanent magnets. Then, thanks to these semi-analytical calculations, the ring dimensions are optimized in order to have a great force or a great stiffness. Moreover, we show that the relative position of the rings for which the force is the strongest depends on the air gap dimension. This result is new because the curvature effect is taken into account in this paper.

We can say that such semi-analytical expressions are more precise than the numerical evaluation of the magnetic forces obtained with the finite element method. Moreover, semi-analytical expressions have a low computational cost whereas the finite element method has a high one. Thereby, as shown in this paper, such calculations allow an easy optimization of quadripolar lenses or devices using permanent magnets.

I. INTRODUCTION

The first studies concerning passive magnetic bearings with permanent magnets have been done by Yonnet [START_REF] Yonnet | Passive magnetic bearings with permanent magnets[END_REF] [START_REF] Yonnet | Permanent magnet bearings and couplings[END_REF]. Such passive magnetic bearings used only permanent magnets radially or axially magnetized. Permanent magnets are commonly used in many electrical devices and engineering applications. Most engineering applications need several ring permanent magnets and the determination of the magnetic force between them is thus required. It is to be noted that several way of obtaining the magnetic forces between rings are possible. Authors generally use either numerical means or 2D analytical calculations for determining magnetic fields or magnetic forces created between ring permanent magnets [3]- [START_REF] Akoun | 3d analytical calculation of the forces exerted between two cuboidal magnets[END_REF]. One of the major problem of a numerical approach lies in the fact that it has a high computational cost. Moreover, numerical algorithms are not so precise as analytical calculations. Alternative solutions are thus required. As it is difficult to obtain fully analytical expressions of the magnetic field created by ring permanent magnets [START_REF] Kim | Mutual inductance of noncoaxial circular coils with constant current density[END_REF]- [START_REF] Kim | Restoring force between two noncoaxial circular coils[END_REF], the determination of the forces between them is still more difficult. Consequently, authors use semi-analytical expressions which are very important steps in order to evaluate either the magnetic field created by ring permanent magnets or the magnetic forces exerted between them. The determination of the magnetic field created by ring permanent magnets have been studied a lot with such approaches [START_REF] Babic | New expression for calculating magnetic fields due to current-carrying solid conductors[END_REF]- [START_REF] Rakotoarison | Using coulombian approach for modeling scalar potential and magnetic field of a permanent magnet with radial polarization[END_REF].

Another way of calculating the forces between ring permanent magnets can be done by using the 2D analytical expressions of the magnetic field created by infinite parallelepipedic magnets. Assuming that the magnetic field created by ring permanent magnets can be approximated by the magnetic field created by infinite parallelepipedic magnets either to the faces of the magnets (2D analytical 1. Representation of the two ring permanent magnets which are radially centered. The inner radius of the outer ring is denoted r in , the outer one is rout, its height is h. The inner radius of the inner ring is denoted r in2 , the outer one is r out2 , its height is z b -z a .

approach) [START_REF] Lemarquand | Nonsinusoidal torque of permanent-magnet couplings[END_REF]- [START_REF] Zhu | Comparison of halbach magnetized brushless machines based on discrete magnet segments or a single ring magnet[END_REF], the magnetic forces can be calculated more easily with the 2D analytical approach.

However, these formula are not valid when the ring radius is small [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]- [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. In any case, the magnetic field created by ring permanent magnets can be determined in terms of elliptic integrals [START_REF] Masson | Electrostatique[END_REF] for ring permanent magnets axially magnetized. Such an approach is appropriate because the algorithms used to calculate elliptic integrals are both very robust and fast. Consequently, the expressions obtained have a very low computational cost (less than 0.2 s to determine the magnetic components of the field created by ring permanent magnets whose polarization is axial).

First, this paper presents useful semi-analytical expressions of the force and the stiffness exerted between two ring permanent magnets whose polarization is radial. Then, this paper explains why we cannot reduce some numerical integrations of the semianalytical expressions of the force between two ring permanent magnets. Eventually, this paper deals with the optimization of the ring dimensions in order to have either a great force or a great stiffness between the rings. We show that the relative position between the rings is of great importance in the design of devices using ring permanent magnets. All the expressions determined in this paper are available online [27].

II. CALCULATION OF THE AXIAL FORCE BETWEEN TWO RING PERMANENT MAGNETS WHOSE POLARIZATION IS RADIAL

This section presents a semi-analytical calculation of the force exerted between two ring permanent magnets whose polarization is radial and which are radially centered. Such configuration corresponds to an axial passive magnetic bearing using ring magnets radially magnetized for mutual attraction. We can say that the devices realized with ring permanent magnets radially magnetized were the first to be built. Moreover, they are usually made of several sections which are stacked together.

A. Notation and geometry

The geometry considered is shown in Fig 1 . A two dimensional representation of the passive magnetic bearing is shown in Fig 2.

The outer radius of the outer ring is r out and the inner one is r in . The outer ring height is h. The outer radius of the inner ring is r out2 and the inner one is r in2 . The inner ring height is z b -z a . It is to be noted that the coulombian model of a permanent magnet is used. Consequently, each ring permanent magnet is represented by both two curved planes which correspond to the inner and outer Fig. 2. Radial magnetic bearing using magnets radially magnetized for mutual attraction. The inner radius of the outer ring is denoted r in , the outer one is r out , its height is h. The inner radius of the inner ring is denoted r in2 , the outer one is r out2 , its height is z b -z a .

faces of the rings which are charged with a surface magnetic pole density σ * and a magnetic pole volume density σ * v . For each case, the inner face is charged with the surface magnetic pole density +σ * and the outer one is charged with the surface magnetic pole density -σ * . It is noted that all the illustrative calculations are done with σ * = J. n = 1T where J is the magnetic polarization vector and n is the unit normal vector which is directed towards 0. Moreover, it is noted that the magnetic pole volume density exits for ring permanent magnets whose polarization is radial in order to to have a charge equilibrium of the ring magnet.

B. Semi-analytical expression of the magnetic force

The two ring permanent magnets which form an axial passive bearing are supposed to be radially centered. Consequently, there is only the axial component of the magnetic force which is exerted between the two rings. We can call it the axial force F z . This axial force can be determined by integrating the magnetic field created by the outer ring on the contributions of the inner one. We must take into account both the magnetic pole surface density and the magnetic pole volume density of each ring. Consequently, as there are two magnetic pole surface densities and one magnetic pole volume density for each ring permanent magnet, we have nine terms for determining the axial force between the two ring permanent magnets. By denoting H z the axial component of the magnetic field produced by the outer ring permanent magnet, Thus, this axial force F z can be written as follows:

F z = (S in ) H z σ * 2 d S - (S out ) H z σ * 2 d S + (V ) H z σ * 2 r 2 d Ṽ (1)
where σ * 2 is the magnetic pole surface density owing to the inner ring permanent magnet, (S in ) is the inner face of the inner ring permanent magnet and (S out ) is the outer face of the inner ring permanent magnet and (V ) is the volume of the inner ring permanent magnet. Therefore, the axial force F z can be expressed as follows:

F z = - 2π θ 1 =0 h z 1 =0 2π θ 2 =0 z b z 2 =z a a(r in , r out2 )dz 1 dθ 1 dz 2 dθ 2 + 2π θ1=0 h z1=0 2π θ2=0 z b z2=za a(r out , r out2 )dz 1 dθ 1 dz 2 dθ 2 + 2π θ 1 =0 h z 1 =0 2π θ 2 =0 z b z 2 =z a a(r in , r in2 )dz 1 dθ 1 dz 2 dθ 2 - 2π θ1=0 h z1=0 2π θ2=0 z b z2=za a(r out , r in2 )dz 1 dθ 1 dz 2 dθ 2 - r out r 1 =r in 2π θ 1 =0 h z 1 =0 2π θ 2 =0 z b z 2 =z a b(r out2 , r 1 )dr 1 dz 1 dθ 1 dz 2 dθ 2 + rout r1=rin 2π θ1=0 h z1=0 2π θ2=0 z b z2=za b(r in2 , r 1 )dr 1 dz 1 dθ 1 dz 2 dθ 2 + r out2 r 2 =r in2 2π θ 1 =0 h z 1 =0 2π θ 2 =0 z b z 2 =z a b(r in , r 2 )dz 1 dθ 1 dr 2 dz 2 dθ 2 - rout2 r2=rin2 2π θ1=0 h z1=0 2π θ2=0 z b z2=za b(r out , r 2 )dz 1 dθ 1 dr 2 dz 2 dθ 2 + r out r 1 =r in r out2 r 2 =r in2 2π θ 1 =0 h z 1 =0 2π θ 2 =0 z b z 2 =z a c(r 1 , r 2 )dr 1 dz 1 dθ 1 dr 2 dz 2 dθ 2 (2) with a(α, β) = σ * 1 σ * 2 4πµ 0 (z 2 -z 1 )αβ (α 2 + β 2 -2αβ + (z 2 -z 1 ) 2 ) 3 2 (3)
where σ * 1 is the magnetic pole surface density owing to the inner ring permanent magnet. and

b(α, β) = a(α, β) β (4) c(α, β) = a(α, β) αβ (5)
The next step is to evaluate (2). If we determined directly F z for given values with all the integrals presented in (2), the computational would be too high because several numerical integrals should be determined at a time. Consequently, we must reduce the number of integrals by integrating analytically a(α, β), b(α, β) and c(α, β) according to the integral variables. Unfortunately, it is not possible to find a fully analytical expression of the magnetic force between two ring permanent magnets but we can use a useful semi-analytical expression which is expressed as follows:

F z = σ * 1 σ * 2 2µ 0 2π θ 1 =0 Sdθ 1 + σ * 1 σ * 2 2µ 0 2π θ 1 =0 M dθ 1 + σ * 1 σ * 2 2µ 0 2π θ1=0 rout2 r2=rin2 V dθ 1 dr 2 (6)
where S corresponds to the magnetic interaction between the magnetic pole surface densities of each ring permanent magnet.

Then, M corresponds to the magnetic interaction between the magnetic pole surface densities of one ring permanent magnet and the magnetic pole volume density of the other one. At least, V corresponds to the magnetic interaction between the magnetic pole volume densities of each ring permanent magnet.

The first contribution S is given by [START_REF] Selvaggi | Calculating the external magnetic field from permanent magnets in permanent-magnet motors -an alternative method[END_REF].

S = -f (z a , z b , h, θ 1 , r in , r out2 ) +f (z a , z b , h, θ 1 , r out , r out2 ) +f (z a , z b , h, θ 1 , r in , r in2 ) -f (z a , z b , h, θ 1 , r out , r in2 ) (7)
with

f (α 1 , α 2 , α 3 , θ 1 , α 5 , α 6 ) = α 5 α 6 log α 3 -α 1 + α 2 5 + α 2 6 + (α 3 -α 1 ) 2 -2α 5 α 6 cos(θ 1 ) +α 5 α 6 log α 1 + α 2 5 + α 2 6 + α 2 1 -2α 5 α 6 cos(θ 1 ) -α 5 α 6 log α 3 -α 2 + α 2 5 + α 2 6 + (α 3 -α 2 ) 2 -2α 5 α 6 cos(θ 1 ) -α 5 α 6 log α 2 + α 2 5 + α 2 6 + α 2 2 -2α 5 α 6 cos(θ 1 ) (8)
The second contribution M is given by [START_REF] Conway | Noncoaxial inductance calculations without the vector potential for axisymmetric coils and planar coils[END_REF].

M = -t(r in , r out , r out2 , h, z a , z b , θ 1 ) +t(r in , r out , r in2 , h, z a , z b , θ 1 ) +t(r in , r out , r in , h, z a , z b , θ 1 ) -t(r in , r out , r out , h, z a , z b , θ 1 ) (9)
with

t(β 1 , β 2 , β 3 , β 4 , β 5 , β 6 , θ 1 ) = t (2) β 1 , β 2 , β 4 -β 5 , β 2 3 + (β 4 -β 5 ) 2 , 2β 3 cos(θ 1 )
+t (2)

β 1 , β 2 , β 5 , β 2 3 + β 2 5 , 2β 3 cos(θ 1 ) -t (2) β 1 , β 2 , β 4 -β 6 , β 2 3 + (β 4 -β 6 ) 2 , 2β 3 cos(θ 1 ) -t (2) β 1 , β 2 , β 6 , β 2 3 + β 2 6 , 2β 3 cos(θ 1 ) (10)
and

t (2) (β 1 , β 2 , q, d, f ) = t (3) (β 2 , q, d, f ) -t (3) (β 1 , q, d, f ) (11)
and

t (3) (s, q, d, f ) = -s + 4d -f 2 -4q 2 2 arctan -f + 2s 4d -f 2 -4q 2 - f 4 log d -q 2 -f s + s 2 +s log q + d -f s + s 2 + q log -f + 2(s + d -f s + s 2 ) - 4d -f 2 -4q 2 + f η log [u 1 ] 4η - -4d + f 2 + 4q 2 + f η log [u 2 ] 4η (12
)

u 1 = - 2 f 2 + 4f q 2 -f 2 (η + 2s) + 4d(-f + η + 2s) q 2 (-4d + f 2 + 4q 2 -f η)(-f + η + 2s) - 8q(-2qs + η d -f s + s 2 ) q 2 (-4d + f 2 + 4q 2 -f η)(-f + η + 2s) ( 13
)
u 2 = - 2 f 2 + 4f q 2 -f 2 (η -2s) -4d(f + η -2s) q 2 (-4d + f 2 + 4q 2 + f η)(f + η -2s) - -8q(2qs + η d -f s + s 2 ) q 2 (-4d + f 2 + 4q 2 + f η)(f + η -2s) (14)
with

η = -4d + f 2 + 4q 2 (15)
The third contribution V is given by [START_REF] Azzerboni | Magnetic field evaluation for disk conductors[END_REF].

V = th (1) (r out , r 2 , z a , z b , h, θ 1 ) -th (1) (r in , r 2 , z a , z b , h, θ 1 ) (16)
with

th (1) = t (3) (r 1 , h -z a , r 2 2 + (h -z a ) 2 , 2r 2 cos(θ 1 )) +t (3) (r 1 , z a , r 2 2 + z 2 a , 2r 2 cos(θ 1 )) -t (3) (r 1 , h -z b , r 2 2 + (h -z b ) 2 , 2r 2 cos(θ 1 )) -t (3) (r 1 , z b , r 2 2 + z 2 b , 2r 2 cos(θ 1 )) (17)

C. Expression of the axial stiffness between two rings whose polarization is radial

The stiffness K exerted between two ring permanent magnets can be determined by calculating the derivative of the axial force with respect to z a . We set z b = z a + b where b is the height of the inner ring permanent magnet. Thus, the axial stiffness K can be calculated with [START_REF] Rakotoarison | Using coulombian approach for modeling scalar potential and magnetic field of a permanent magnet with radial polarization[END_REF].

K = - ∂ ∂z a F z ( 18
)
where F z is given by (2). We obtain :

K = K S + K M + K V ( 19
)
where K S corresponds to the stiffness determined with only the magnetic pole surface densities of each ring permanent magnet. Then, K M corresponds to the stiffness determined with the magnetic interaction between the magnetic pole surface densities of one ring permanent magnet and the magnetic pole volume density of the other one. At least, K V corresponds to the stiffness determined with the magnetic interaction between the magnetic pole volume densities of each ring permanent magnet. Thus, the first contribution K S is expressed as follows:

K S = η 31 1 √ α 31 K * - 4r 3 r 1 α 31 - 1 √ β 31 K * - 4r 3 r 1 β 31 + 1 √ δ 31 K * - 4r 3 r 1 δ 31 - 1 √ γ 31 K * - 4r 3 r 1 γ 31 +η 41 1 √ α 41 K * - 4r 4 r 1 α 41 - 1 √ β 41 K * - 4r 4 r 1 β 41 + 1 √ δ 41 K * - 4r 4 r 1 δ 41 - 1 √ γ 41 K * - 4r 4 r 1 γ 41 +η 32 1 √ α 32 K * - 4r 3 r 2 α 32 - 1 √ β 32 K * - 4r 3 r 2 β 32 + 1 √ δ 32 K * - 4r 3 r 2 δ 32 - 1 √ γ 32 K * - 4r 3 r 2 γ 32 +η 42 1 √ α 42 K * - 4r 4 r 2 α 42 - 1 √ β 42 K * - 4r 4 r 2 β 42 + 1 √ δ 42 K * - 4r 4 r 2 δ 42 - 1 √ γ 42 K * - 4r 4 r 2 γ 42 (20)
with

η ij = 2r i r j σ * µ 0 ( 21
)
α ij = (r i -r j ) 2 + z 2 a ( 22
)
β ij = (r i -r j ) 2 + (z a + h) 2 (23)
γ ij = (r i -r j ) 2 + (z a -h) ( 24
)
δ ij = (r i -r j ) 2 + (b -h) 2 + z a (2b -2h + z a ) (25)
K * [m] = π 2 0 1 1 -m sin(θ) 2 dθ ( 26
)
The second contribution K M is expressed as follows:

K M = σ * 1 σ * 2 2µ 0 2π θ=0 udθ ( 27
)
with

u = f (r in , r out , r in2 , h, z a , b, θ) -f (r in , r out , r out2 , h, z a , b, θ) +f (r in2 , r out2 , r in , h, z a , b, θ) -f (r in2 , r out2 , r out , h, z a , b, θ) (28)
and

f (α, β, γ, h, z a , b, θ) = -γ log α -γ cos(θ) + α 2 + γ 2 + z 2 a -2αγ cos(θ) +γ log α -γ cos(θ) + α 2 + γ 2 + (z a + b) 2 -2αγ cos(θ) +γ log α -γ cos(θ) + α 2 + γ 2 + (z a -h) 2 -2αγ cos(θ) -γ log α -γ cos(θ) + α 2 + γ 2 + (b -h) 2 + 2z a (b -h) + z 2 a -2αγ cos(θ) +γ log β -γ cos(θ) + β 2 + γ 2 + z 2 a -2αγ cos(θ) -γ log β -γ cos(θ) + β 2 + γ 2 + (z a + b) 2 -2αγ cos(θ) +γ log β -γ cos(θ) + β 2 + γ 2 + (z a -h) 2 -2αγ cos(θ) -γ log β -γ cos(θ) + β 2 + γ 2 + (b -h) 2 + 2z a (b -h) + z 2 a -2αγ cos(θ) (29)
The third contribution K V is expressed as follows:

K V = σ * 1 σ * 2 2µ 0 2π θ=0 r out r 1 =r in δdθ (30) with δ = -log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + z 2 a -2r 1 r in2 cos(θ) + log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + (z a + b) 2 -2r 1 r in2 cos(θ) -log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + (b -h) 2 + 2bz a -2hz a + z 2 a -2r 1 r in2 cos(θ) + log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + (z a -h) 2 -2r 1 r in2 cos(θ) + log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + z 2 a -2r 1 r out2 cos(θ) -log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + (z a + b) 2 -2r 1 r out2 cos(θ) -log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + (z a -h) 2 -2r 1 r out2 cos(θ) + log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + (b -h) 2 + 2bz a -2hz a + z 2 a -2r 1 r out2 cos(θ) (31)
We can say that the expression of the axial stiffness can be determined analytically if we only take into account the magnetic pole surface densities of each ring.

III. DISCUSSION ABOUT THE POSSIBILITY OF REDUCING THE NUMBER OF NUMERICAL INTEGRATIONS FOR THE AXIAL FORCE

EXPRESSION

The aim of this section is to explain why we cannot reduce the number of numerical integrations of the semi-analytical expression of the force exerted between two ring permanent magnets radially magnetized. They correspond physically to three kinds of interactions between two ring permanent magnets whose polarization is radial.

A. Interaction between the contributions of the surface densities

The first kind of physical interaction is due to the surface contributions of each ring permanent magnet. This physical interaction corresponds to the case when the surface densities of the outer ring are integrated with the surface densities of the inner one. In Eq.

(2), these surface contributions correspond to the integrand denoted a(α, β) where α and β can be r in , r out , r in2 and r out2 . Let us consider the integrand a(α, β). As we can see in Eq. ( 2), the integration variables depend on θ 1 , θ 2 , z 1 and z 2 . The integration according θ 2 does not change the form of the integrand a(α, β) because a(α, β) does not depend on θ 2 . Consequently, we can say that the form of the force F s between the surface density contributions of each ring permanent magnet is given as follows :

F s = θ1 z1 z2 a 1 (z 2 -z 1 ) (a 2 -a 3 cos(θ 1 ) + (z 2 -z 1 ) 2 ) 3 2 dθ 1 dz 1 dz 2 (32)
where a 1 , a 2 and a 3 are constant. For example, we can have

a 1 = 2π(h-z a ), a 2 = 2π(r 2 in +r 2 out +(h-z a ) 2
) and a 3 = 2π(-2r in r out2 ) for our illustration here. In short, these parameters are given by (2) and correspond to the case when four integrals must be determined.

After having integrated according z 1 and z 2 , we obtain a semi-analytical expression with only one numerical integration whose variable depends on θ 1 :

F s = θ 1 log a 4 + a 5 + a 6 cos(θ 1 ) dθ 1 (33)
These parameters are given by ( 8) and they are not a function of the angle θ. It is to be noted that (33) cannot be integrated analytically.

B. Interaction between the contributions of the surface and volume densities

The second kind of physical interaction is due to both the surface and the volume contributions of each ring permanent magnet.

This physical interaction corresponds to the case when the surface densities of each ring permanent magnet are integrated with the volume densities of the other ones. In Eq. ( 2), these surface contributions correspond to the integrand denoted b(α, β) where α and β can be r in , r out , r in2 , r out2 , r 1 and r 2 . Let us consider the integrand b(α, β). As we can see in Eq. ( 2), the integration variables depend on θ 1 , θ 2 , z 1 , z 2 , r 1 and r 2 . The integration according θ 2 does not change the form of the integrand b(α, β) because b(α, β)

does not depend on θ 2 . Consequently, we can say that the form of the force F vs between the surface density contributions and the volume density contributions of each ring permanent magnet is given as follows :

F vs = θ 1 r i z 1 z 2 b 1 (z 2 -z 1 ) (b 2 -b 3 cos(θ 1 ) + (z 2 -z 1 ) 2 ) 3 2 dθ 1 dz 1 dz 2 dr i ( 34
)
where b 1 , b 2 and b 3 are constant. These parameters are given by ( 2) and correspond to the case when five integrals must be determined.

After having integrated according z 1 , z 2 and r 1 , we obtain a semi-analytical expression with only one numerical integration whose variable depends on θ 1 : this semi-analytical expression is in fact given by ( 33) where s depends on θ 1 . We see that (34) owns a term which has the same form as the one presented in Eq.( 33). Consequently, it cannot be integrated analytically as well.

C. Interaction between the contributions of the volume densities

The third kind of physical interaction is due to the volume contributions of each ring permanent magnet. This physical interaction corresponds to the case when the volume densities of each ring permanent magnet are integrated together. In Eq. ( 2), this volume contribution correspond to the integrand denoted c(α, β) where α and β can be r 1 or r 2 . Let us consider the integrand c(α, β). As we can see in Eq. ( 2), the integration variables depend on θ 1 , θ 2 , z 1 , z 2 , r 1 and r 2 . The integration according θ 2 does not change the form of the integrand c(α, β) because c(α, β) does not depend on θ 2 . Consequenlty, we can say that the form of the force F v between the volume density contributions is given as follows :

F v = θ 1 r 1 r 2 z 1 z 2 (z 2 -z 1 ) (c 1 -c 2 cos(θ 1 ) + (z 2 -z 1 ) 2 ) 3 2 dθ 1 dz 1 dz 2 dr 1 dr 2 (35)
where c 1 and c 2 are constant. These parameters are given by ( 2) and correspond to the case when six integrals must be determined.

After having integrated according z 1 , z 2 and r 1 , we obtain the same integrand as the one obtained previously in the case of the study of the force exerted between the surface and volume charge contributions. Consequently, we deduct that we cannot integrate analytically F v according to θ 1 . Moreover, the analytical integration according to r 2 does not seem possible. As a conclusion, we can say that the obtaining of a fully analytical expression of the force between two ring permanent magnets whose polarization is radial does not seem possible but a semi-analytical expression can be used to determine this axial force.

IV. OPTIMIZATION OF THE INNER RING PERMANENT MAGNET DIMENSIONS

This section discusses the optimal dimensions of the rings in order to have either a great force or a great stiffness.

A. Influence of the air gap dimension on the force and the stiffness

First, we study the influence of the air gap dimension on the force and the stiffness between the two ring permanent magnets. For this purpose, we represent the axial force versus the axial displacement of the inner ring for different air gaps in Fig. 3. It is noted that, in our configuration, the air gap corresponds to the difference between r in and r out2 . Furthermore, the width and the height of each ring permanent magnet are constant. We take r in = 0.025m, r out = 0.028m, r out2 -r int2 = 0.003m, J = 1T, h = 0.003m, zb -z a = 0.003m. Fig. 3 shows three important points.

First, we see that the smaller the air gap between the ring permanent magnets is, the greater the axial force is. Consequently, it is necessary to have the smaller air gap between ring permanent magnets if a great force is searched. This result is well-known. It was shown with the two dimensional approach.

Second, we see that the exact position of the maximal force exerted between two ring permanent magnets depends slightly on the air gap dimension. This result is new because our study uses a three-dimensional approach of the magnetic force whereas the previous ones used a two-dimensional approach. The magnet curvature must be taken into account in order to obtain precisely the position of the maximal force exerted between two ring magnets. The exact position of the maximal force is represented in Sti ness [N/m] Fig. 5. Representation of the axial stiffness of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different air gaps; r in = 0.025m, r out = 0.028m,r out2 -r int2 = 0.003m, J = 1T, h = 0.003m, zb -z a = 0.003m.

B. Determination of the optimal height of the inner ring permanent magnet

Another parameter which can be optimized is the height of the inner ring permanent magnet. To do so, the axial component of the magnetic force is plotted versus the axial displacement of the inner ring permanent magnet for several inner ring heights in Fig. 6.

The values taken for the parameters are still the same as the previous ones. Fig. 6 shows that the axial component of the magnetic force is the greatest if its height equals the outer ring height. Indeed, if the ring inner height is smaller than the ring outer one, the smaller the inner ring width is, the smaller the axial component of the magnetic force is. If the ring inner height is greater than the ring outer one, the greater the inner ring width is, the smaller the axial component of the magnetic force is. Consequently, if a great force is searched, the two ring permanent magnets must have the same height. Furthermore, we see that the gradient of the curves is the more important when the inner ring height is the same as the outer one. This result is consistent with Fig 7 where the axial stiffness is represented versus the axial displacement of the inner ring permanent magnet.

C. Determination of the optimal width of the inner ring permanent magnet

The third parameter which can be optimized in our configuration is the width of the inner ring permanent magnet. To do so, the axial component of the magnetic force is plotted versus the axial displacement of the inner ring permanent magnet for several inner ring widths in Fig. 8. The values taken for the parameters are the same as the previous ones. Fig. 8 shows that the greater the inner ring width is, the greater the axial component of the magnetic force is. However, it is noted that a compromise in the ring dimensions must be found because the cost of the magnet must be taken into account. A good compromise can be found as follows : if the inner ring width equals two times its height, the axial component of the magnetic force is 72N when z = 0.0015m. If the inner ring width equals three times its height, the axial component of the magnetic force is 73N. Consequenlty, we deduct that it is not necessary to have an inner ring width which is greater than two times its height.

The optimal stiffness depends also on the inner ring width. To see that, we have represented in

V. OBTAINING THE BEST CONFIGURATION

The previous section shows that the ring dimensions must be optimized in order to create a very good passive magnetic bearing.

The air gap must be the smallest, the ring heights must be the same and we have shown that the inner ring width must equal two times its height. All these parameters have been determined with a outer ring whose cross-section is a square. However, we can also optimize the outer ring in order to improve the passive bearing. By taking into account the optimal dimensions found in the previous section, we can compare three configurations. The first one, shown in Fig 10-A, consists of two rings whose cross-section is a square.

VI. CONCLUSION

This paper has presented new three-dimensional semi-analytical expressions allowing us to determine both the axial force and the axial stiffness between two ring permanent magnets whose polarization is radial. This paper also discusses the reason why we cannot find a fully analytical expression of the axial component of the magnetic force. Eventually, we have discussed the optimal dimensions of the ring permanent magnets which allow us to have either a great axial force or a great axial stiffness. We have shown that a good compromise can be found when the cross-section of a ring is a rectangle. However, the relative position of the two rings is not constant but depends on the air gap dimension. Such results can be very useful for people involved in the design of magnetic bearings.

  Fig.1. Representation of the two ring permanent magnets which are radially centered. The inner radius of the outer ring is denoted r in , the outer one is rout, its

Fig. 3 .Fig. 4 .

 34 Fig.3. Representation of the axial component of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different air gaps; r in = 0.025m, r out = 0.028m,r out2 -r int2 = 0.003m, J = 1T, h = 0.003m, zb -z a = 0.003m.

Fig 4 .

 4 Such result is very useful because it clearly shows that if a great axial force is searched, the relative height between the two rings inner ring depends on the air gap dimension. Eventually, Fig 3 shows that the stiffness depends greatly on the air gap dimension. Fig 3 shows that the smaller the air gap dimensionis, the greater the stiffness is because the gradient of the curve is the more important for small air gaps. This result is consistent with the representation of the axial stiffness versus the axial displacement of the inner ring permanent magnet (Fig5). Indeed, we see that the smaller the air gap dimension is, the greater the axial stiffness is. Moreover, we see that when the axial force is maximal in Fig 3, the axial stiffness equals zero in Fig 5, which is still consistent.

Fig. 6 .

 6 Fig.6. Representation of the axial component of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different heights; r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, za = 0m.

Fig. 7 .Fig. 8 .

 78 Fig.7. Representation of the axial component of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different heights; r in = 0.025m, r out = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, z a = 0m.

  Fig 9 the axial stiffness versus the inner ring width when z a = 0. Fig 9 shows that the larger the inner ring width is, the greater the axial stiffness is.

Fig. 9 .Fig. 10 .

 910 Fig.9. Representation of the axial stiffness exerted between two ring permanent magnets versus inner ring width r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, zb -z a = 0.003m.

Figs 11 and 12

 12 Figs 11 and 12 show that the best configration is the configuration B presented in Fig 10. However, the relative height between the two ring permanent magnets depends on the air gap dimension (Fig 4). Consequently, this last parameter must be taken into account in the design of passive bearings using ring permanent magnets.

 12. Representation of the axial stiffness exerted between two ring permanent magnets versus inner ring width r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, zb -z a = 0.003m.

[3] B. Azzerboni, E. Cardelli, and A. Tellini, "Computation of the magnetic field in massive conductor systems," IEEE Trans. Magn., vol. 25, pp. 4462-4473, november 1989.