Svenja Hermann
email: svenja.hermann@femto-st.fr

Pauline Butaud

Rafael De O Teloli

Jean-François Manceau

Maxime Savary

Gaël Chevallier
email: gael.chevallier@femto-st.fr

Magnetically induced friction damping based on magnetoactive elastomers -A proof of concept

Keywords: friction damping, magnetoactive elastomer (MAE), sandwich structure, backbone identification

In this work we present a new concept for friction damping which is based on the permanent magnetization and the flexibility of magnetoactive elastomers (MAE). The concept is called "magnetofriction" and is studied on a sandwich structure in the present work. The assembled structure consists of a MAE core and ferromagnetic skins. Due to the magnetization of the MAE, the skins are attracted to the core and the sandwich is assembled without the need of glue or other means of fixation. The magnetic forces generate the normal loading in the contact interfaces between skins and core while a tangential loading is induced by forced vibrations applied to the sandwich beam. Random and harmonic vibration tests reveal the nonlinear dynamic behavior of the sandwich structure and highlight the damping capacity of magnetofriction. The harmonic balance method is used on a nonlinear dynamic model of the system to identify the parameters related to the resonance frequency and damping ratio. The results obtained from the identification state damping ratios up to 9 % under harmonic excitation, which is remarkable for friction damping. In addition, a Finite Element-based model of the structure is implemented to study the local contact states in the sandwich. Multi-physical simulations are performed to analyze the magnetic and mechanical loading in the contact interface. The local results are used to evaluate the overall contact states in the interfaces which determine the absence or presence of friction damping. The simulations show that the vibration amplitude changes the contact state.

Introduction

Excessive vibrations in mechanical structures may generate undesired noise or even mechanical failure of components in the long term. Since quiet operation and durability are priority requirements during project design stages, different damping technologies have been developed during the last years. These technologies include passive devices, e.g. viscoelastic dampers, or active devices, e.g. piezoelectric actuators coupled with a control circuit. From an energetic point of view, passive devices have advantages over active devices, since they do not require an external power supply. In this work we present an original passive concept for vibration damping based on friction in an assembled structure which contains magnetoactive elastomers.

Magnetoactive elastomers are viscoelastic materials and thus capable to dissipate vibrations into heat efficiently. The viscoelastic damping is independent from the vibration magnitude and is therefore classified as a linear mechanism. However, the main drawback of viscoelastic materials is their sensitivity to temperature and frequency variations [START_REF] Butaud | In-core heat distribution control for adaptive damping and stiffness tuning of composite structures[END_REF]. Until today, magnetoactive elastomers have been studied exclusively in the context of the adaptive vibration control by external magnetic fields [START_REF] Zhou | Use of magnetorheological elastomer in an adaptive sandwich beam with conductive skins. part ii: Dynamic properties[END_REF][START_REF] Nayak | Dynamic analysis of magnetorheological elastomer-based sandwich beam with conductive skins under various boundary conditions[END_REF][START_REF] Nayak | Vibration analysis of a three-layer magnetorheological elastomer embedded sandwich beam with conductive skins using finite element method[END_REF][START_REF] Hu | Experimental investigation of the vibration characteristics of a magnetorheological elastomer sandwich beam under non-homogeneous small magnetic fields[END_REF][START_REF] Szmidt | Adaptive damping of a double-beam structure based on magnetorheological elastomer[END_REF][START_REF] Dyniewicz | Semi-active control of a sandwich beam partially filled with magnetorheological elastomer[END_REF].

Friction is one of the main sources of damping in assembled structures. It is induced by relative motion between mechanically linked components. In the contact interface of the joints, the normal stress distribution is governed by the tightening load [START_REF] Popp | Vibration damping by friction forces: theory and applications[END_REF]. The shear stress field is induced by the vibration loading. Depending on the ratio between normal stress and tangential stress, the contact state can range from total stick to partial sliding to mainly sliding [START_REF] Scheel | Challenging an experimental nonlinear modal analysis method with a new strongly friction-damped structure[END_REF]. When partial sliding occurs in the contact interfaces, the vibration energy is dissipated by friction [START_REF] Scheel | Challenging an experimental nonlinear modal analysis method with a new strongly friction-damped structure[END_REF][START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF]. The dissipation behavior is strongly nonlinear as it is linked to the threshold between the three contact states. The advantages of friction damping are the independence from the excitation frequency and the low dependence on temperature variations.

The design of structures that benefit from frictional dissipation remains challenging due to the nonlinearity. Depending on the system, a variation of the normal stress, also called contact pressure, can lead to a change of the modal stiffness [START_REF] Bouaziz | Vibration reduction of an assembly by control of the tightening load[END_REF] or the modal mass [START_REF] Vakilinejad | A comparison of robustness and performance of linear and nonlinear lanchester dampers[END_REF]. Brake [START_REF] Brake | The mechanics of jointed structures: recent research and open challenges for developing predictive models for structural dynamics[END_REF] and Teloli et al. [START_REF] De O. Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The orion beam[END_REF] observed that a decrease of the normal tightening load in structures that are assembled by bolted joints results in a higher frictional damping and also in softening effects (decrease of the contact stiffness) in vicinity of resonance frequencies. In addition, the contact pressure distribution is influenced by geometric irregularities on the macroscopic scale (e.g. uneven or wavy surfaces) and on the microscopic scale (roughness) [START_REF] Peyret | Dynamic damping in joints: multiscale model taking into account defects in a nominally plane surface[END_REF][START_REF] Brake | Observations of variability and repeatability in jointed structures[END_REF]. An interesting possibility to establish a constant contact pressure is the use of magnetically induced forces. Attractive forces between permanent magnets [START_REF] Dai | A new permanent magnetic friction damper device for passive energy dissipation[END_REF] or between permanent magnets and ferromagnetic materials [START_REF] Choi | Smart damper using the combination of magnetic friction and pre-compressed rubber springs[END_REF] have been used to establish the contact in joints of vibrating systems in the context of seismic damping. In the recent works of Amjadian et al. [START_REF] Amjadian | Modeling, design, and testing of a proofof-concept prototype damper with friction and eddy current damping effects[END_REF][START_REF] Amjadian | Planar arrangement of permanent magnets in design of a magneto-solid damper by finite element method[END_REF] the contact pressure is generated by repulsive forces between two permanent magnets. The force is transmitted to a friction pad and its intensity is controlled by the distance between the magnets.

In this work we present a new concept for friction damping in assembled structures. It is based on the use of permanently magnetized elastomer composites, also known as magnetically Hard MagnetoActive Elastomers (H-MAEs) [START_REF] Hermann | Magnetic and dynamic mechanical properties of a highly coercive mre based on ndfeb particles and a stiff matrix[END_REF]. H-MAEs are composed of magnetic particles and a viscoelastic matrix. Their magnetic behavior is comparable to the behavior of permanent magnets [START_REF] Hermann | Magnetic and dynamic mechanical properties of a highly coercive mre based on ndfeb particles and a stiff matrix[END_REF]. Their elastic modulus, however, is much lower (order of magnitude: 2 MPa to 10 MPa). The key idea of the concept is to use the permanent magnetization and the high flexibility of H-MAEs in the joints of an assembled structure to obtain an evenly distributed contact pressure that favors energy dissipation by friction. In the present study, the H-MAE is used in the core of a sandwich structure with ferromagnetic skins. The mechanical link is established by the magnetic forces.

The use of flexible magnets instead of stiff magnets is favorable for friction damping as the H-MAE can adapt its shape to the deformation of the skins. The contact is thus maintained even for a relatively high deformation of the structure and the assembly is less sensitive to geometric irregularities. The H-MAE is also less heavy than a permanent magnet with comparable volume. In addition, our concept offers the possibility to combine friction damping with viscoelastic damping.

In order to assess the vibration performances, a so-called Magnetofriction Sandwich Structure (MSS) is subjected to random vibrations and harmonic vibrations around the three first natural frequencies with different excitation amplitudes. The equivalent modal parameters of the test structure are evaluated with an original identification method and reveal a nonlinear behavior of the structure. A numerical model of the experiment is used to show the influence of the excitation amplitude on the local contact state between the core and the skins.

The paper is organized as follows: section 2 presents the concept of magnetically induced friction damping and describes the assembly of the experimental test sample. Section 3 contains the description of the test setup and details the experimental investigations and results. Section 4 contains a numerical study of the local contact parameters that are not accessible in the experiment. The last section consists of conclusions on magnetically induced friction damping and opens perspectives for future studies.

Magnetofriction in viscoelastic joints

At the beginning of this section, the concept that allows the combination of friction damping and viscoelastic damping is introduced. The experimental test sample which is used to proof the concept is presented subsequently.

Concept

The new concept is based on magnetically Hard MagnetoActive Elastomers (H-MAEs) which are used as core of a sandwich structure where they establish the mechanical linkage of the assembly. Fig. 1 illustrates the use of the H-MAE to achieve friction damping in sandwich structures. Fig. 1 (a) sketches an oscillating sandwich beam structure. It is composed of a magnetized H-MAE core which is enclosed by two ferromagnetic skins. The magnetic flux, generated by the magnetized H-MAE, propagates into the magnetically permeable skins (Fig. 1 (b)). The ferromagnetic skins are attracted by the magnetized H-MAE and the magneto-mechanical interactions maintain the skins clamped to the core without additional means of fixation (Fig. 1 (c)). The contact pressure is hence established by the magnetization. When the structure is subjected to transverse vibrations, the sandwich structure is deformed and shear stresses arise in the contact interfaces (Fig. 1 (d)). Normal and shear stresses are now present in the contact area. Their ratio determines the presence and the magnitude of sliding and therefore of friction damping. For small vibration amplitudes, the normal stresses prevail and the contact is maintained without sliding in the interface. In this case, the vibration damping is due to the viscoelasticity of the H-MAE composite. With an increasing vibration amplitude, the tangential stresses increase. Throughout the contact area, more and more parts begin to slide and friction damping occurs. The damping is maximized for a specific ratio of sticking and sliding which is determined by the vibration amplitude. By further increasing the amplitude, the tangential forces become predominant until the contact surfaces mainly slide. In this case, the friction damping is less effective [START_REF] Berthillier | Blades forced response analysis with friction dampers[END_REF][START_REF] Ye | Torsional friction damper optimization[END_REF]. The concept is called magnetofriction, since contact pressure is established magnetically.

Test sample

The concept is proven experimentally with the test sample sketched in Fig. 2 (a).

The assembly is called Magnetofriction Sandwich Structure (MSS) in the following. It comprises H-MAE patches, steel skins and an aluminum spacer.

The H-MAE is composed of a two-part silicone elastomer MED-4014, provided by NuSil T M , and of magnetically hard NdFeB particles (MQFP-14-12 from Magnequench T M ). Details about the elaboration process can be found in [START_REF] Hermann | Magnetic and dynamic mechanical properties of a highly coercive mre based on ndfeb particles and a stiff matrix[END_REF]. The H-MAE in the present study has a volume fraction of 36 % of NdFeB particles and shows a residual flux density of 0.3 T. The patches used in this study are magnetized in their direction of thickness. The two skins of the MSS consist of the steel alloy M400-50A which is ferromagnetic and magnetically soft. The magnetization curve of the material can be found in the Appendix D (Fig. D.1). The varying width of the sandwich minimizes the effects of the clamp on the experimental results in terms of damping and stiffness. The wider part has three drill holes of which two are used for the alignment of the sandwich in the clamping. The third drill hole is not used in the present study.

The MSS is assembled such as an aluminum spacer and the H-MAE patches form the core between the steel skins. A sample holder is used to ensure the alignment of the components during the assembly. The aluminum spacer is placed in the region of the MSS that is clamped in the experimental setup (Fig. 2 (c)). The stiffness of the spacer makes the clamping less sensitive to the clamping force. The spacer does not contribute to the generation of the contact pressure as it is magnetically neutral. The spacer is maintained in place by the clamping force that fixes the MSS in the experimental setup. The H-MAE core is positioned adjacent to the spacer in the narrow part of the sandwich. The core does not consist of one continuous H-MAE layer since the size of the magnetization chamber limits the patch size. The upper part of the core is therefore composed of five magnetized H-MAE patches. They are oriented such as the magnetization of two adjacent patches is opposed (Fig. 2 (b)). It is impossible to position patches with the same direction of magnetization next to each other in the MSS since their outer regions repel each other due to the magnetic forces in this configuration.

In the narrow part of the sandwich, the magnetic forces between the H-MAE patches and the metal sheets maintain the structure assembled. No glue or bolts are necessary for the assembly and the structural integrity. To highlight the benefits of the magnetic forces compared to more classical adhesion forces, the structure has also been tested in a configuration where the patches were partly glued. A glue point has been applied to the center of each patch on each side that is in contact with the skin. The amount of glue remained small and it was not applied on the entire patch surface (Appendix F, Fig. F.1). This approach was chosen on the basis of additional experiments that revealed an accumulation of the glue in the intersections between H-MAE patches. The glued configuration is called gMSS.

Vibration experiments 200

In this section, the experimental investigations on magnetically induced friction damping are presented. The experimental setup is explained at the beginning of the section. Afterwards, the random and harmonic excitation methods, which are used to characterize the dynamic behavior of the sandwich structure, are detailed. The two different excitation strategies are il-205 lustrated by the experimental results and first conclusions on the dynamic behavior are drawn. In the subsequent section, the identification of the modal parameters for the two different excitation methods is presented.

Experimental setup

In the experimental setup, the sandwich structure is clamped to an electromagnetic shaker (Tira TV 52125 200 N) as shown in Fig. 3 (a). The necessary power for the shaker is delivered by an amplifier. A triaxial B&K DeltaTron R accelerometer (Type 4520) is glued to the clamp to drive the base-excitation tests. The acceleration component a 2 (t) measured in direction e 2 is used as reference signal in the post-processing. The structure velocity is measured with a contactless Polytec vibrometer PSV-500 Xtra R . The velocity measurements are either performed on a single point with the position (x 1 , x 3 ) on the MSS or on multiple points. In the second case, a velocity field v 2 (x 1 , x 3 , t) is obtained in direction e 2 . For the random excitation tests, the data acquisition and the post-processing are performed by the software of the laser vibrometer (Fig. 3 (b)). During the harmonic oscillation, the excitation amplitude is feedback-controlled by in-house subroutines [START_REF] De O. Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The orion beam[END_REF]. The algorithms are implemented in Matlab R and transmitted via a NI9234 hardware.

Random excitation

The first part of the experiment consists of random vibration tests with the MSS and the gMSS. The goals of these tests are the identification of the first three vibration modes and the identification of the extended modal parameters, natural frequency f 0 and damping ξ, that depends on the vibration amplitude. The frequencies of the modes are used to define the frequency ranges for the harmonic excitation tests and the modal parameters are used to compare the MSS with the gMSS. The two different RMS acceleration levels 0.1 m/s 2 and 1 m/s 2 are used for the excitation. The structures are subjected to random vibrations in a frequency range between 20 Hz and 370 Hz. During each test, single-point measurements of the velocity are performed on one point which is represented in the schematic view of Fig. 4 (a). The frequency response functions (FRFs) are estimated for each frequency f on the measurement points using the H 1 estimator (Eq. ( 1)). The estimator is calculated from the average autopower S aa of the input signals â(f ) (acceleration) and the average crosspower S av of the same input signals and the output signals v(f ) (velocity) obtained from n repetitions of the measurement, yielding:

H 1 (f ) = S av (f ) S aa (f ) = 1 n n j=1 âj (f )v H j (f ) 1 n n j=1 âj (f )â H j (f ) , (1)
where (•) H denotes the complex conjugate transpose. The FRFs have been estimated for each of the four random excitation tests. The result of the MSS for the high excitation level is shown in Fig. 4 (b) as an example. The first three resonances are located in the measurement range and are well separated.

Harmonic excitation

In the second part of the experiment, the MSS is subjected to harmonic vibrations with different acceleration magnitudes. The goal of this study is to highlight the nonlinear behavior of the structure with respect to the excitation amplitude experimentally. A controlled step-sine procedure is used for the excitation. During this type of test, the MSS is excited in the steady-state regime with a fixed base-acceleration magnitude, whereas the excitation frequency is incrementally varied around each resonance frequency. Due to the acceleration and frequency control, this procedure is well adapted to identify the nonlinear behavior of vibrating structures. The procedure is described in detail in [START_REF] De O. Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The orion beam[END_REF]. In the present study, the three resonance modes of the MSS As an example for the step-sine tests, amplitude and phase of the FRFs of the third resonance mode for different acceleration amplitudes are shown in Fig. 5 (a) and Fig. 5 (b), respectively. The amplitude represents the ratio of the velocity measured on the structure and the acceleration measured on the base. The results indicate a nonlinear behavior of the MSS with respect to the excitation amplitude: the resonance frequency decreases as the base acceleration amplitude increases and the shape of the amplitude peaks is modified. Under these excitation conditions, note that the amplitude peaks of the FRFs at the resonance frequencies tend to decrease and incline to the left.

The so-called softening-effect takes place in assembled structures due to the presence of partial sliding that is responsible for reducing the contact stiffness.

This nonlinear behavior is related to the friction dissipation [START_REF] Berthillier | Blades forced response analysis with friction dampers[END_REF][START_REF] Ye | Torsional friction damper optimization[END_REF][START_REF] Vakilinejad | A comparison of robustness and performance of linear and nonlinear lanchester dampers[END_REF][START_REF] De O. Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The orion beam[END_REF]. The results shown in Fig. 5 prove the existence of friction in the MSS and hence the concept of magnetofriction. Limiting frequencies for sliding and adhesion have been obtained by numerical modeling of the experiment. The limits are represented in Fig. 5. Two resonance frequencies are not included in the interval due to the fact that the clamping was not exactly represented in the model. However, the experimental and numerical results are in good accordance and a varying contact state is confirmed. A comparison of the tests with increasing and decreasing frequency reveals jumps in vicinity of the resonances (Fig. 5). In fact, the assembled configuration of the MSS depends on the last sliding between the core and the skins. This sliding occurs around the resonance state, thus the assembled state is

Identification of modal parameters

For the random vibration tests, a conventional fitting procedure of the vibrometer software is used to identify the natural frequencies f 0 (i) and modal damping ξ (i) of the first three modes. The results are summarized in Tab. 2. All natural frequencies of the MSS decrease significantly as the excitation level increases while the natural frequencies of the gMSS decrease slightly. The evolution of f 0 (i) is related to the stiffness which depends mainly on the slip in the interfaces in our case. Once slip occurs, the stiffness decreases. The resulting decrease in the stiffness of the structure ("softening" effect) is responsible for the shift of f 0 (i) . Sliding is prevented by the glue in the gMSS;

conversely, the amount of sliding is higher in the MSS, which is required for friction damping and therefore, the assembly technique based on magnetic forces is more favorable for friction damping than the use of glue. In addition, the glue slightly increases the stiffness of the structure. As a result, the resonance frequencies are higher for the gMSS than for the MSS. The results are less obvious for damping because the effect of stiffness nonlinearity on structures under random excitation is to widen the frequency bandwidth [START_REF] Caughey | Equivalent linearization techniques[END_REF]. Damping has the same effect: it is impossible to separate the sources of damping sources, therefore it is difficult to conclude. This is not a problem since the harmonic tests allow us to evaluate the dependence of the damping on the amplitude of the vibration.

Sample Excitation Under harmonic excitation, the experimental results indicate that the stiffness and the damping evolve in a nonlinear manner. Thus, for identifying the equivalent modal parameters, each bending mode of the MSS is modeled based on the following nonlinear differential equation:

1 st mode 2 nd mode 3 rd mode a RM S [m/s 2 ] f 1 0 [Hz] ξ 1 [%] f 2 0 [Hz] ξ 2 [%] f 3 0 [Hz] ξ 3 [%] MSS 0.
m (i) q(t) + c (i) q(t) + k (i) q(t) + β (i) q(t)q 2 (t) + α (i) q(t) 3 = Ca sin(Ωt), (2)
where Ω = 2πf corresponds to the angular excitation frequency. The parameters m (i) , c (i) , k (i) correspond to modal values of mass, damping and stiffness, respectively, of the i th vibrating mode; β (i) and α (i) are the modal parameters that govern the nonlinear behavior of the structure. They are related to the nonlinear damping and nonlinear stiffness, respectively. The identification procedure performed to estimate the model parameters is formulated into a harmonic balance framework considering a single harmonic.

Further details related to this identification strategy are addressed in [START_REF] De O. Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The orion beam[END_REF]. It is important to point out that, in the present case, the acceleration and the velocity have been measured and the excitation force is not known. Thus, the method has been improved to process tests considering a base-acceleration excitation. A preliminary step is to derive the excitation force from the acceleration measured at the base.

Additional full-field measurements show that the mode shape, and hence the modal mass m (i) , does not change significantly with the acceleration magnitude (Appendix B, Fig. B.1). Moreover, modal damping is composed of a linear viscous part and a nonlinear viscous one that depends on the square of the displacement. The nonlinear stiffness force is composed of linear part and a cubic contribution. Based on the harmonic balance formulation, the resulting expressions of damping and frequency are given by:

f (i) 0 = 1 2π
k (i) + 3 4 α (i) Q 2 m (i) , ξ (i) = 1 πm (i) f (i) 0 c (i) + β (i) Q 2 4 , (3)
where Q is the fundamental harmonic amplitude of the system response.

Figure 6 depicts the results of the post-processed identification as the dependence of the natural frequency and damping on the vibration amplitude. These graphs are usually presented under the name "backbone" in the literature [START_REF] Peyret | A framework for backbone experimental tracking: Piezoelectric actuators, stop-sine signal and kalman filtering[END_REF].

The identification process has led to build the backbones corresponding to the first three modes of the structure. As for all the structures that implement a friction-induced damping solution, the natural frequencies of each mode decrease with the vibration amplitude [START_REF] Brake | The mechanics of jointed structures: recent research and open challenges for developing predictive models for structural dynamics[END_REF]. For this new concept, so-called magneto-friction structures, the dependence is stronger than in most of the cases studied [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF][START_REF] De O. Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The orion beam[END_REF][START_REF] Lacayo | Updating structural models containing nonlinear iwan joints using quasi-static modal analysis[END_REF]. This reflects the fact that a large amount of strain energy is mobilized in the interface. A Modal Strain Energy Fraction (MSEF) can be defined by analogy with the electro-mechanical coupling [START_REF] Chevallier | Piezoceramic shunted damping concept: testing, modelling and correlation[END_REF] or the viscoelastic induced damping criteria [START_REF] Hammami | Numerical design and test on an assembled structure of a bolted joint with viscoelastic damping[END_REF], which is given by:

K (i) = max(f (i) 0 ) -min(f (i) 0 ) min(f (i) 0 ) . ( 4
)
From the MSEF point of view, the setup benefits from the fact that the interface is distributed over the entire domain. Concerning the damping, the results are remarkable since more than 8% of damping was obtained on the first mode. Only the assembly studied by [START_REF] Scheel | Challenging an experimental nonlinear modal analysis method with a new strongly friction-damped structure[END_REF], the Rubbing Beam Resonator (RuBber) structure, produced better results. In both cases, MSS and RuBber, it is the ability to produce relatively uniform contact pressures that has allowed us to achieve these excellent performances. To illustrate how remarkable the performance obtained with the MSS structure is, Table 3 presents a comparison in terms of MSEF and damping between the MSS and other structures described in the literature. The RubBer experiment presents comparable performances, even slightly superior. It uses a very original device to manage the dry friction. But technologically compared to RubBer, the device used by MSS is really much simpler and easier to implement. Brake-ReußBeam (BRB) and ORION are using bolted joints. Their performance is definitely not as good. MAIAS [START_REF] Hammami | Numerical design and test on an assembled structure of a bolted joint with viscoelastic damping[END_REF] uses viscoelastic joints that allows to get similar performances in terms of damping but with a lower MSEF. Piezo-R uses a piezoelectric device shunted with a resistance that leads the same result than for a optimized bolted structure like ORION. Thus, this magneto-friction technology appears to be a promising technology for vibration damping.

Case MSS MSS MSS ORION

Numerical study of the local contact state

The experimental results show a dependence of the resonance modes on the acceleration amplitude which is characteristic for systems with friction damping ( [START_REF] Vakilinejad | A comparison of robustness and performance of linear and nonlinear lanchester dampers[END_REF][START_REF] Berthillier | Blades forced response analysis with friction dampers[END_REF][START_REF] Ye | Torsional friction damper optimization[END_REF]). This type of dissipation is caused by partial sliding in the contact interface at the positions where the shear stress is superior to the normal stress. It is extremely difficult to analyze the local contact state experimentally; hence a numerical model of the sandwich structure is implemented in FE code. The so-called digital twin is used to study the stress contributions that determine the contact state: the magnetically induced normal stress on the one hand and the vibration induced shear stress on the other hand.

Implementation of the digital twin

In the present study, the digital twin is a two-dimensional representation of the sandwich (Fig. 7). The digital twin has the same geometric and material properties as the experimental sample. The Maxwell stress tensor is commonly used to calculate the magnetic loading. Force calculations based on this tensor can be performed on material interfaces that separate magnetizable and non-magnetizable materials [START_REF] Hirsinger | Les tenseurs de contraintes électromagnétiques [electromagnetic stress tensors[END_REF]. For this reason, a gap is required between the sandwich core and each of the skins. The size of the gap set to 0.3 mm. According to the parametric study shown in Appendix E, this value is the best compromise between the calculation time and the accuracy of the local quantities obtained by the FE analysis. To establish the mechanical contact in the sandwich, however, a non-magnetic elastic material is situated in the gaps. The patches are arranged with alternating direction of magnetization in the model similar to the experiment. The components of the sandwich are meshed with a quadratic mesh (element size 0.1 mm -0.2 mm). A triangular mesh with an element size between 0.1 mm and 2 mm is used for the air surrounding the sandwich. The mesh is refined in vicinity of the sandwich to increase the accuracy of the local magnetic quantities.

The H-MAE patches are permanently magnetized and the variations of the magnetic field H remain small. Under these conditions, the composite material can be described by the linear magnetic material model shown in Eq. ( 5). This material model relates the magnetic flux density B to H using the magnetic constant µ 0 , the magnetic permeability µ r and the residual flux density B r :

B = µ 0 µ r H + B r . (5)
The nonlinear magnetic behavior of the steel skins is directly implemented as behavior law in the numerical model (cf. Appendix D, Fig. D.1). The magnetically neutral behavior of the aluminum spacer, the surrounding air and the non-magnetic elastic material can also be described by Eq. ( 5)

with B r = 0.

The mechanical strain remains small in the test cases; therefore a linear mechanical behavior of the skins and the core is implemented in the model. In the behavior law, the strain tensor ε is linked to the stress tensor σ by the fourth order stiffness tensor C:

σ = C : ε . (6)
The materials are modeled as mechanically isotropic and the corresponding stiffness tensors are hence determined by their Young's moduli E and their Poisson's ratios ν. The material in the gap between the core and the skins possesses the same mechanical properties as the H-MAEs. All magnetic and mechanical material parameters are listed in table 4.

In the numerical simulations, a structural magneto-mechanical coupling is implemented to calculate the magnetic forces between the components of the sandwich. The calculation of the local magnetic fields accounts for the influence of the shape of the different parts of the sandwich on their internal magnetic fields [START_REF] Ivaneyko | Mechanical properties of magneto-sensitive elastomers: unification of the continuum-mechanics and microscopic theoretical approaches[END_REF]. MAEs are generally known for their magneto-elastic material coupling represented by a dependence of their elastic modulus on externally applied magnetic fields [START_REF] Jolly | The magnetoviscoelastic response of elastomer composites consisting of ferrous particles embedded in a polymer matrix[END_REF] or magnetostriction [START_REF] Danas | Experiments and modeling of iron-particle-filled magnetorheological elastomers[END_REF]. coupling is strong in magnetically soft MAEs when the material is not magnetically saturated [START_REF] Galipeau | The effect of particle shape and distribution on the macroscopic behavior of magnetoelastic composites[END_REF]. In the present case, however, the material coupling can be neglected for several reasons. The MAE in the present study is magnetically hard and initially magnetized. The magnetization has no significant influence on the elastic modulus of the MAE used in the present study [START_REF] Hermann | Magnetic and dynamic mechanical properties of a highly coercive mre based on ndfeb particles and a stiff matrix[END_REF]. This behavior is due to the comparatively high stiffness of the elastomer matrix, which has a high influence on the intensity of magneto-elastic couplings in MAE [START_REF] Castañeda | Homogenization-based constitutive models for magnetorheological elastomers at finite strain[END_REF]. In addition, no external magnetic field is applied and the variation of the magnetic field in the sandwich remains small during the tests due to the small modal displacements. Furthermore, the magnetized MAE has a low magnetic permeability in the direction of magnetization. This state is comparable to the magnetic saturation regime of magnetically soft MAEs for which the magneto-elastic coupling effect is weaker than prior to the saturation [START_REF] Galipeau | The effect of particle shape and distribution on the macroscopic behavior of magnetoelastic composites[END_REF]. For these reasons, uncoupled magnetic and the mechanical behavior laws are used to model the MAE material.

Study cases

Two numerical study cases are performed to determine the magnetically induced normal stress and the shear stress induced by vibrations, respectively. In order to evaluate the stresses, the sliding in the contact interfaces is blocked. For the magnetic part of the simulations, performed in a quasistatic regime, the absence of electric fields and in-plane magnetic fields are assumed. The mechanical part of the simulations assumes in-plane mechanical strains and a stationary regime. The remaining magneto-mechanical governing equations are stated in the Appendix D.

The first study case concerns the magneto-mechanical loading that maintains the sandwich structure assembled. Fully coupled magneto-mechanical sim-ulations have been performed considering the attraction force between the skins and the H-MAE as the only loading. In the first step, the surface force density of magnetic origin, f mag S , is calculated on the boundaries of the skins and the core. The assumptions regarding the absence of currents and the quasi-static regime allow to calculate the forces by the help of the Maxwell stress tensor T in the non-magnetic materials surrounding the skin and the surface normal n on the boundaries of the magnetic components ( [START_REF] Hirsinger | Les tenseurs de contraintes électromagnétiques [electromagnetic stress tensors[END_REF]):

f mag S = T • n = 1 µ 0 B ⊗ B - B • B 2µ 0 I • n , (7)
with ⊗, • and I representing the vector product, the scalar product and the identity matrix, respectively. The force density is then applied as mechanical loading on the boundaries of the skins and the core. The resulting normal stress σ 22 in the contact interface, corresponds to the component of f mag S in the direction e 2 :

σ 22 = f mag S • e 2 . (8)
The load case of the second study is purely mechanical. The third bending mode is chosen for the numerical study in this paper since the full-field velocity was easier to measure for this mode than for the two other ones.

The experimental modal velocity field of the third bending mode has been obtained from the multi-point measurements with three different acceleration magnitudes presented in the Appendix B. Based on the experimental results obtained for a base acceleration of a = 1 m/s 2 , a displacement field function compatible with the numerical mesh is generated by interpolation The displacement fields are applied on the upper skin on the outside of the sandwich, where they have been measured in the experiment, in three different load cases as stationary mechanical loading. During the study, the central line in the spacer is not allowed to move in the e 1 direction and the bottom edge of the sandwich is fixed.

Results of normal and shear stress

The normal stress σ 22 in the interface between the upper skin and the nonmagnetic elastic layer is represented in Fig. 8. The normal stress is mainly negative which means that the skin is attracted by the H-MAE core. Attraction forces were also obtained for the second skin. In the middle region of the patches and the spacer, the stress magnitude is nearly constant. Stronger variations of the stress are observed in vicinity of geometric or material discontinuities, like corners and material interfaces. The distribution of the normal stress in the interface is directly related to the intensity and the distribution of the magnetic flux density (cf. Eq. ( 7)). The magnetic flux concentrates near the intersections of magnetic patches where it forms dense closed loops (cf. Appendix E, Fig. E.1). As the stress increases with the magnitude of the flux density, it is higher at an intersection of two patches than in the middle of a patch. Furthermore, the distribution of the magnetic flux generates a low normal stress in the clamping and an increased normal stress at the tip of the sandwich. Even though the contact pressure varies locally, it is well distributed in the contact interface compared to common mechanical assemblies, such as bolted joints (cf. [START_REF] Brake | Observations of variability and repeatability in jointed structures[END_REF] for example).

In addition, the pressure has a relatively weak magnitude which is of great technological interest as it offers the possibility to obtain partial sliding (and hence friction damping) at low vibration magnitudes.

The shear stress σ 12 in the contact interface of interest is shown in Fig. 8

Evaluation of the contact state

A Coulomb friction model (cf. [START_REF] Bowden | The friction and lubrication of solids[END_REF] for example) is used to study whether the contact surfaces show a relative movement or not (sliding or sticking). To analyze the contact state, the stress ratio |σ 12 |/|σ 22 | is determined for the three acceleration magnitudes. The result is shown in Fig. 9. In one of our studies on the friction behavior of the H-MAE, a friction coefficient of 1 has been measured experimentally. The value is highlighted in Fig. 9 by a red line and will be used to evaluate the contact state in the following. For the lowest excitation magnitude of 1 m/s 2 , the main part of the curve lies below c µ , thus the contact surfaces are mainly sticking together. For the intermediate excitation magnitude of 20 m/s 2 , the ratio |σ 12 |/|σ 22 | is of the same order of magnitude as c µ . Parts of the contact surfaces are sticking while others are sliding. Friction damping occurs principally for this contact state. For the highest excitation magnitude of 73 m/s 2 , the main part of the curve lies above c µ which means that there is mainly sliding in the contact interface.

The combination of experimental results, numerical simulations and the friction model allows to study the evolution of the contact in the digital twin of the sandwich structure. Since the concept of H-MAE based friction damping is fairly new, the model comprises several assumptions in its current state. The 2D representation and the assumption of in-plane strain and magnetic fields allow to reduce the CPU time but the physical concept of the present model can be implemented in the same manner in a 3D model. Simulations with 3D models will be helpful to get more detailed information, for example on the role of boundary effects, which are likely to influence the threshold at which friction damping is activated. The model furthermore assumes a negligible convection effect, i.e. transport terms, for the magnetic quantities in the stationary vibration regime. This assumption can be verified in future studies as well as the effect of non-stationary electro-magnetic effects, such as eddy currents [START_REF] Amjadian | Modeling, design, and testing of a proofof-concept prototype damper with friction and eddy current damping effects[END_REF]. However, in its current state it has already provided useful information on the normal stress distribution and the evolution of the contact state

Conclusions and perspectives

In the present work, we introduce the concept of magnetofriction for vibration damping. In this concept, magnetic forces are used to generate the contact pressure while vibrations induce the tangential forces in the same contact interface.

The concept is demonstrated by the so-called Magnetofriction Sandwich Structure (MSS). The MSS comprises a core that consists of a permanently magnetized MagnetoActive Elastomer (MAE) and ferromagnetic steel skins. Due to the magnetic attraction forces between core and skins, the MSS can be assembled without glue or other additional mechanical means of fixation.

Compared to other assembly strategies (e.g. use of bolts, screws, magnets), the flexible MAE in the core of the MSS adapts its shape to the deformations of the vibrating structure. This shape adaptability assures a distributed contact interface which enhances friction damping. Random and harmonic vibrations are performed with the MSS to evaluate its amplitude-dependent characteristics. The backbones are identified into an original harmonic balance framework. The MSS shows a strongly nonlinear behavior for the first three modes. With increasing vibration magnitude, the natural frequency decreases, whereas the modal damping increases. This behavior is characteristic for friction damping. It is also shown that the Modal Strain Energy Fraction is very high compared to many other damping technologies. Moreover, the damping ratio is one of the highest obtained among all the studies dedicated to friction induced damping. It is remarkable that the best value has been obtained on the first eigenmode. The setup benefits from the fact that the contact pressure is low (about 2.4 kPa compared to 1 MPa for the ORION beam) and almost uniformly distributed over the entire domain.

A numerical model of the experiments is implemented for the study of the local contact state. The simulations show that the magnetic contact pressure is almost uniformly distributed in the contact interface which is favorable for friction damping. The comparison of normal and tangential forces in the contact interface confirms the dependence of the contact states on the vibration amplitude: from total stick for a small vibration magnitude to mainly sliding for high vibration magnitudes.

The present work can be classified as a proof of concept for magnetofriction. A topic for future works is the separation of the damping phenomena: viscoelastic damping, friction damping and eddy current damping. Another perspective for further inquisition lies in evaluating wear and damage over the entire life of the structure. The current sandwich is quite heavy because of the NeFeB particles. Instead of using a full H-MAE core, the core could be positioned where the structure deforms the most. The results of these studies will allow the design of magnetofriction devices with optimum dissipation capacities. from 1, but the variation is very small. In conclusion, it is stated that the nonlinearity does not affect the mode shape in a significant manner. The figures below reveal the comparison of the FRFs measured experimentally and computed with the Harmonic Balance Method from the identification process described in section 3.4. The amplitude represents the ratio of the displacements measured on the structure and on the base [mm/mm].

Figure 1 :

 1 Figure 1: Notational diagram of a part of an oscillating sandwich beam structure (a) and a zoom on one part of the structure sketching the magnetic state (b), the clamping force (c) and the contact situation (d) in the beam.

Figure 2 :

 2 Figure 2: Components of the MSS (a), schematic side view of the assembled MSS (b) and schematic view of the clamping (exploded view) (c).

Figure 3 :

 3 Figure 3: Experimental setup with MSS specimen mounted on the electromagnetic shaker in side view: photograph (a) and schematic view (b).

Figure 4 :

 4 Figure 4: Schematic front view of the clamped MSS (a) and estimation of the FRFs by H 1 (b) for the MSS with an RMS acceleration magnitude of 1 m/s 2 .

Figure 5 :

 5 Figure 5: Results of harmonic vibration test: amplitudes and phases of the FRFs obtained for eight different acceleration amplitudes with either increasing -△-or decreasing -▽frequency during the measurement.

  modified slightly during each step-sine test. The effect is present and visible between two tests with the same base-acceleration magnitude and different frequency sweeping directions. The results of the two first vibration modes can be found in the Appendix A (Fig. A.1 and Fig. A.2). The identification of the modal parameters f 0 and ξ of the harmonic excitation tests is presented in the second part of section 3.4.

Figure 6 :

 6 Figure6: Backbones calculated for the resonance frequency f (i) and the damping ξ (i) for the first three resonance mode of the MSS. The displacement denotes the vibration magnitude.

Figure 7 :

 7 Figure 7: Numerical model of the sandwich structure, axis ratio e 2 :e 1 = 10:1 in graphic representation for better visibility.

(

  cf. Appendix D, Fig. D.2). The displacement functions of the two additional vibration amplitudes are obtained by multiplying the displacement function by the corresponding vibration level. This method can be used as mode shapes are similar regardless of the vibration magnitude (cf. Fig. B.1).

  (b) for the three studies with different base acceleration magnitudes. The qualitative evolution of the shear stress in the interface is similar for the three cases since the mechanical loading is generated by one displacement field function with different amplitudes. The shear stress is low near the fixation, where the modal displacements are small and the structure is hence not deformed by the loading. The direction of the shear stress changes at the points where the displacement of the structure is at a local maximum or minimum (cf. Appendix D, Fig. D.2).

  According to the Coulomb friction model, sliding occurs when the ratio of tangential and normal stress are equal to a friction coefficient: |σ 12 |/|σ 22 | = c µ . In the numerical model, however, the relative movements are blocked. In this case, sliding is represented by a stress ratio of |σ 12 |/|σ 22 | that is equal to or greater than c µ : Contact state : |σ 12 | < c µ |σ 22 | → sticking |σ 12 | ≥ c µ |σ 22 | → sliding (

Figure 8 :

 8 Figure 8: Distribution of the magnetic flux density B, direction of H-MAE magnetization (red arrows) and evaluation interface (green line) in the model (a), magnetically induced normal stress σ 22 (b) and mechanically induced shear stresses for the three tested excitation magnitudes (c) in the evaluation interface.

Figure 9 :

 9 Figure 9: Distribution of the stress magnitude ratio |σ 12 |/|σ 22 | for different acceleration amplitudes and the friction coefficient c µ = 1 in the evaluation interface of the sandwich structure.

Figures A. 1

 1 Figures A.1 and A.2 depict the FRFs estimated for the second and third 580

Figure A. 1 :

 1 Figure A.1: Experimental results of harmonic vibration test for the first resonance mode: amplitudes and phases of the FRFs obtained for six different acceleration amplitudes with increasing -△-or decreasing -▽-frequency during the measurement.

Figure A. 2 :

 2 Figure A.2: Experimental results of harmonic vibration test for the second resonance mode: amplitudes and phases of the FRFs obtained for seven different acceleration amplitudes with decreasing -▽-frequency during the measurement.

Figure B. 1 :

 1 Figure B.1: Normalized mode shapes of the structure for the three excitation amplitudes (a) and modal assurance criterion (MAC) for the three configurations (b).

Figure C. 1 :

 1 Figure C.1: Results of the identification process ("Simulated") for the harmonic vibrations of the first mode in comparison with experimental results.

Figure C. 2 :

 2 Figure C.2: Results of the identification process ("Simulated") for the harmonic vibrations of the second mode in comparison with experimental results.

Figure C. 3 :

 3 Figure C.3: Results of the identification process ("Simulated") for the harmonic vibrations of the third mode in comparison with experimental results.

Figure E. 1 :

 1 Figure E.1: Local magnetic flux density B in the sandwich structure (cones, color corresponding to magnitude ||B||) and orientation of the residual flux density B r of the magnetic patches (red arrows).

Table 1 :

 1 Summary

of the acceleration magnitudes âk and frequency range [f 1 -f n ] used the step-sine tests for each resonance mode. x -applied during test

Table 2 :

 2 Resonance frequencies f 0 of the first three vibration modes obtained for the random excitation with different excitation levels with the MSS the gMSS.

Table 3 :

 3 Bench-marking between several results from the literature. Modal Strain Energy Fraction indicates the involved energy. Modal damping indicates the target performance.

	RubBeR BRB MAIAS Piezo-R

Table 4 :

 4 Material properties of the different materials modeled in the digital twin.

	The material

* cf. Appendix D, Fig. D.1

Acknowledgments

This work has been funded by the French company Moving Magnet Technologies SA and EIPHI Graduate School, ANR-17-EURE-0002.

CRediT authorship contribution statement

Declaration of conflicting interests

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. In a second analysis, the velocity fields, measured for the different configurations, are compared point by point. For each configuration, the field contains N velocities measured on the different points j which, at the resonance, form an approximation of the k th eigenvector Φ k as shown in Eq. ((B.1)), where the exponent T indicates the transpose of the vector:

To compare the vectors between the three different configurations, the coefficients of a Modal Assurance Criterion (MAC) matrix are calculated by

where k and l correspond to one of the three configurations. A coefficient M AC kl of the matrix represents the correlation between velocity fields of the configuration k and l. If the coefficient is equal to 1, the eigenvectors are co-linear and that the velocity fields are identical. A coefficient equal to 0 means, that the eigenvectors are oriented in an orthogonal manner and that the velocity fields are very different. The 3x3 MAC matrix for the three tested configurations is shown in

Appendix D.2. Mechanical equations

The mechanical balance equations are implemented for stationary cases for which Cauchy's equation of motion is described by

where σ is the Cauchy stress tensor.

Appendix D.3. Prescribed displacement in the mechanical simulations

The prescribed displacement field is computed from the experimental data obtained with laser full field measurement and reported in Appendix C.

Appendix E.1. Distribution of the magnetic flux density in the sample

The magnetization of the H-MAE patches generates a magnetic flux that propagates in the sandwich structure. The orientation and intensity of local magnetic flux density are represented in Fig. E.1 by the orientation and the color of the cones. The orientation of the residual flux density B r of the patches is represented by red arrows in addition. The magnitude of the magnetic flux density is very high in the skins compared to the other parts and reaches 2.15 T at the maximum. The high magnetic permeability of the ferromagnetic material causes the flux to concentrate in the skins. The color gradient shows that the intensity of the flux density is not evenly distributed in the skins: the magnetic flux is concentrated at the intersection between two patches while the remaining section of the patch is only weakly magnetized (0.1 T on average). The weak magnetization is due to the material properties and the geometry of the patches, while the concentration at the intersection is caused by the magnetic circuit. At the intersection of two oppositely magnetized patches, closed loops of the magnetic flux form easily and the flux lines are dense.

Appendix E.2. Convergence study of the normal load

The contact pressure is governed by the computation of the magnetic forces. As mentioned in section 4, an artificial gap has been created as the conventional FE code is able to calculate the magnetically induced pressure only 655 between magnetic and non-magnetic materials. The sensitivity analysis in Fig. E.2 shows that the local pressure decreases with the air-gap, which is comprehensible as the magnetic field decreases with the distance towards the magnetized MAE. The choice of 0.3mm was based on the pressure in the central regions of the 660 five patches as the pressure peaks arise near singularity points of the geometry. The zoom on these regions in Fig. E.3 shows that the pressure decreases slightly stronger than before at distances above 0.3mm. The computation time, on the other hand, increased strongly with the reduction of the air gap choice of 0.3mm is a compromise between the calculation time accuracy of the force calculation.


domain: info math spi stat

text:
Paul Saves
email: paul.saves@isae-supaero.fr

Youssef Diouane
email: youssef.diouane@polymtl.ca

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Thierry Lefebvre
email: thierry.lefebvre@onera.fr

J Morlier
email: joseph.morlier@isae-supaero.fr

A general square exponential kernel

I. Nomenclature

𝑛

= number of continuous variables 𝑚 = number of integer variables 𝑙 = number of categorical variables Ω ∈ R 𝑛 = continuous space 𝑆 ∈ Z 𝑚 = integer space F 𝑙 = categorical space 𝐿 𝑖 , 𝑖 ∈ {1, . . . , 𝑙} = number of levels for the 𝑖 𝑡 ℎ categorical variable 𝜃 𝑐𝑜𝑛𝑡 = vector of hyperparameters for the continuous part of the Gaussian process model 𝑘 = correlation kernel 𝜃 𝑐𝑜𝑛𝑡 𝑗 , 𝑗 ∈ {1, . . . , 𝑛 + 𝑚} = hyperparameter for the 𝑗 𝑡 ℎ continuous or integer variable 𝑅 𝑐𝑜𝑛𝑡 = correlation matrix for continuous and integer inputs Θ 𝑐𝑎𝑡 = hyperparameters for the categorical part of the Gaussian process model 𝐾 𝑖 = categorical kernel for the 𝑖 𝑡 ℎ categorical variable Θ 𝑖 = matrix of hyperparameters for the 𝑖 𝑡 ℎ categorical variable 𝑅 𝑐𝑎𝑡 = correlation matrix for categorical inputs Θ = [Θ 𝑐𝑎𝑡 , 𝜃 𝑐𝑜𝑛𝑡 ] = hyperparameters for the Gaussian process model 𝑅 = 𝑅 𝑐𝑜𝑛𝑡 𝑅 𝑐𝑎𝑡 = correlation matrix for mixed integer inputs

II. Introduction

N ew aircraft configurations with a lower footprint on the environment (also known as Eco-aircraft design) have seen a resurgence of interest [START_REF] Duriez | HALE multidisciplinary design optimization with a focus on eco-material selection[END_REF][START_REF] Priem | An efficient application of Bayesian optimization to an industrial MDO framework for aircraft design[END_REF][START_REF] Ciampa | The AGILE Paradigm: the next generation of collaborative MDO[END_REF]. In this context, one targets to minimize the footprint on the environment of the aircraft using a Multidisciplinary Design Analysis (MDA) [START_REF] Lambe | A unified description of MDO architectures[END_REF][START_REF] Lambe | Extensions to the design structure matrix for the description of multidisciplinary design, analysis, and optimization processes[END_REF][START_REF] Martins | Multidisciplinary Design Optimization: A Survey of Architectures[END_REF]. This is an example of an expensive-to-evaluate without derivative problem that could be encountered on industry. Therefore, it could be useful to use a surrogate model that simplifies by a lot an expensive model and gives a good approximation from a small data set of known configurations. An example of an industrial application of surrogate model in the context of aircraft design is given in Fig. 1. Nevertheless, in this context, the process generally involves mixed continuous-categorical design variables. For instance, the size of aircraft structural parts can be described using continuous variables; in case of thin-sheet stiffened sizing, they represent panel thicknesses and stiffening cross-sectional areas. The set of discrete variables can encompass design variables such as the number of panels, the list of cross sectional areas or the material choices. In this work, we target to construct an inexpensive surrogate model f for a black-box simulation function of the form

𝑓 : Ω × 𝑆 × F 𝑙 → R. (1)
The function 𝑓 is typically expensive-to-evaluate simulations with no exploitable derivative information. Ω ⊂ R 𝑛 represents the bounded continuous design set for the 𝑛 continuous variables. 𝑆 ⊂ Z 𝑚 represents the bounded integer set where 𝐿 1 , . . . , 𝐿 𝑚 are the numbers of levels of the 𝑚 quantitative integer variables on which we can define an order relation and F 𝑙 = {1, . . . , 𝐿 1 } × {1, . . . , 𝐿 2 } × . . . × {1, . . . , 𝐿 𝑙 } is the design space for 𝑙 categorical qualitative variables with their respective 𝐿 1 , . . . , 𝐿 𝑙 levels. In this context, Gaussian processes (GP) [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] Rasmussen | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF][START_REF] Forrester | Engineering Design via Surrogate Modelling: A Practical Guide[END_REF][START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF], also called Kriging models, are known to be a good modelling strategy to define response surface models. Namely, we will consider that our unknown black-box function 𝑓 is a realization of an underlying GP of mean 𝜇 𝑓 and of standard deviation 𝑠 𝑓 , i.e.,

𝑓 ∼ f = GP 𝜇 𝑓 , [𝑠 𝑓 ] 2 .
(

) 2
For a general problem involving categorical or integer variables, several modeling strategies to build a mixedcategorical GP have been proposed [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF][START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF][START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF][START_REF] Rufato | Creating Recommender Systems for Industrial Engineering Problems Using a Mixed Categorical-Continuous Data-Driven Method[END_REF][START_REF] Cuesta-Ramirez | A comparison of mixed-variables Bayesian optimization approaches[END_REF]. Compared to a continuous GP, the major changes are in the estimation of the correlation matrix: the latter is essential for building estimates of 𝜇 𝑓 and 𝑠 𝑓 . Similarly to the process of constructing a GP with continuous inputs, relaxation techniques [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF] and Gower distance based models [START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF] use a kernel-based approach to estimate the correlation matrix. Other recent approaches try to estimate the correlation matrix directly independently of a kernel choice [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF][START_REF] Cuesta-Ramirez | A comparison of mixed-variables Bayesian optimization approaches[END_REF] which shows good results as these methods model completely the correlations. However, the direct estimation of the correlation matrix as proposed in [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF] is not adapted for high-dimensional problems as it is very expensive to compute all the required hyperparameters. In fact, dimension reduction methods such as principal components analysis (known as 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][START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by Partial Least Squares dimension reduction[END_REF], Kriging model with Partial Least Squares) require the construction of the correlation matrix via a kernel function. KPLS models are used to reduce to number of hyperparameters and to handle a large number of mixed inputs [START_REF] Saves | Bayesian optimization for mixed variables using an adaptive dimension reduction process: applications to aircraft design[END_REF].

In this work, we target to extend the classical paradigm used for continuous inputs to cover the mixed categorical case. We will present a kernel-based approach that will lead to a unified approach for existing approximation methods [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF][START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF][START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF]. A similar process for the estimation of the hyperparameters could be applied to both continuous and categorical inputs. The good potential of the proposed approach is shown over a set of analytical test cases. In this paper, a particular attention will be given to the Gaussian kernel, but our proposed approach can be straightforwardly extended to other existing kernels.

The reminder of this paper is as follows. In Section III, a detailed review of the GP model for continuous and for categorical inputs is given. The extended kernel-based approach for constructing the correlation matrix is presented in Section IV. Section V presents academical tests as well as the obtained results. Conclusions and perspectives are finally drawn in Section VI.

III. Gaussian process meta-models for mixed categorical inputs

In general, a GP model is used to fit a response surface model from an initial set of points, known as the Design of Experiments (DoE) [START_REF] Rasmussen | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF][START_REF] Forrester | Engineering Design via Surrogate Modelling: A Practical Guide[END_REF][START_REF] Kim | Surrogate-Based Optimization for Mixed-Integer Nonlinear Problems[END_REF]. The GP provides a mean response hypersurface as well as a pointwise estimation of its variance. In what comes next, let 𝑛 𝑡 be the size of DoE data set (𝑥, 𝑦 𝑓 ).

A. GP meta-models for continuous inputs

In this subsection, we will only consider that all the design variables are continuous in problem (1): namely, the design space will be restricted to Ω ⊂ R 𝑛 . In this case, we assume that 𝑓 : R 𝑛 → R is defined only over the continuous design space: 𝑛 𝑡 is the number of already evaluated points in R 𝑛 of the deterministic function 𝑓 and ∀𝑟 ∈ {1, .., 𝑛 𝑡 }. Let 𝑥 𝑟 = (𝑥 𝑟 1 , ..., 𝑥 𝑟 𝑛 ) ∈ R 𝑛 be the 𝑟 𝑡 ℎ point with its respective 𝑛 continuous variable values and 𝑦 𝑓 𝑟 ∈ R be the associated values of 𝑓 (𝑥 𝑟 ) and denote the DoE as (𝑥, 𝑦 𝑓 ). The stochastic model [START_REF] Duvenaud | Automatic model construction with Gaussian processes[END_REF] writes as: f (𝑥) = 𝜇(𝑥) + 𝜖 (𝑥) ∈ R with 𝜖 the error term between 𝑓 and the model approximation 𝜇(𝑥). The error terms are considered as independent and identically distributed random variables of variance 𝜎 2 . Let 𝑅 𝑐𝑜𝑛𝑡 be the correlation matrix between the input points whose elements are defined by

[𝑅 𝑐𝑜𝑛𝑡 ] 𝑟 ,𝑠 = 𝐶𝑜𝑟𝑟 (𝜖 (𝑥 𝑟 ), 𝜖 (𝑥 𝑠 ))
The correlation function 𝐶𝑜𝑟𝑟 is computed using a kernel function 𝑘 that relies on 𝑛 hyperparameters 𝜃 𝑐𝑜𝑛𝑡 estimated typically using maximum likelihood estimator (MLE) [START_REF] Rossi | Mathematical statistics: an introduction to likelihood based inference[END_REF]: 𝐶𝑜𝑟𝑟 (., .) = 𝑘 (., ., 𝜃 𝑐𝑜𝑛𝑡 )

Let 𝑟 𝑐𝑜𝑛𝑡 (𝑥 * ) = (𝐶𝑜𝑟𝑟 (𝜖 (𝑥 * ), 𝜖 (𝑥 1 )), . . . , 𝐶𝑜𝑟𝑟 (𝜖 (𝑥 * ), 𝜖 (𝑥 𝑛 𝑡 ))) for a given 𝑥 * and 1 be the 𝑛 𝑡 vector of ones, then, we have:

𝜇 𝑓 (𝑥 * ) = μ 𝑓 + 𝑟 𝑐𝑜𝑛𝑡 (𝑥 * ) 𝑇 [𝑅 𝑐𝑜𝑛𝑡 ] -1 (y 𝑓 -1 μ 𝑓 ), (3)
and

[𝑠 𝑓 ] 2 (𝑥 * ) = [ σ 𝑓 ] 2 1 -𝑟 𝑐𝑜𝑛𝑡 (𝑥 * ) 𝑇 [𝑅 𝑐𝑜𝑛𝑡 ] -1 𝑟 𝑐𝑜𝑛𝑡 (𝑥 * ) + (1 -1 𝑇 [𝑅 𝑐𝑜𝑛𝑡 ] -1 𝑟 𝑐𝑜𝑛𝑡 (𝑥 * )) 2 1 𝑇 [𝑅 𝑐𝑜𝑛𝑡 ] -1 1 , (4)
where μ 𝑓 and σ 𝑓 , respectively, are the MLE of 𝜇 and 𝜎 with respect to 𝜃 𝑐𝑜𝑛𝑡 given the DoE data set (𝑥, 𝑦 𝑓 ). In these formulae, f and [𝑠 𝑓 ] 2 both depend on 𝑅 and 𝑟 which are characterized by the correlation kernel 𝑘 (., ., 𝜃 𝑐𝑜𝑛𝑡 ). For two continuous inputs 𝑥 𝑟 and 𝑥 𝑠 , the Gaussian kernel is defined as:

𝑘 (𝑥 𝑟 , 𝑥 𝑠 , 𝜃 𝑐𝑜𝑛𝑡 ) = 𝑛 𝑗=1 exp -𝜃 𝑐𝑜𝑛𝑡 𝑗 𝑥 𝑟 𝑗 -𝑥 𝑠 𝑗 2 = 𝑛 𝑗=1 exp -|𝑥 𝑟 𝑗 -𝑥 𝑠 𝑗 |𝜃 𝑐𝑜𝑛𝑡 𝑗 |𝑥 𝑟 𝑗 -𝑥 𝑠 𝑗 | (5)
Other kernels can be used like the Matérn 3/2 kernel [START_REF] Sacks | Designs for Computer Experiments[END_REF]:

𝑘 (𝑥 𝑟 , 𝑥 𝑠 , 𝜃 𝑐𝑜𝑛𝑡 ) = 𝑛 𝑗=1 1 + √ 3𝜃 𝑐𝑜𝑛𝑡 𝑗 𝑥 𝑟 𝑗 -𝑥 𝑠 𝑗 exp - √ 3𝜃 𝑐𝑜𝑛𝑡 𝑗 𝑥 𝑟 𝑗 -𝑥 𝑠 𝑗 (6)
As the hyperparameters are always multiplied by the distance between two points, they can be interpreted as being the inverse correlation length. A GP whose kernel is stationary and based on a distance between two points is well adapted to a continuous context and its extension to categorical or integer variables is not straightforward. In the next part, for the general mixed-categorical case, we will consider both continous and categorical hyperparameters Θ.

B. GP meta-models for mixed-categorical inputs

In this subsection, we are considering the case where the design variables could be categorical or integer. Namely, we assume that 𝑓 : R 𝑛 × Z 𝑚 × F 𝑙 ↦ → R. The set F 𝑙 = {1, . . . , 𝐿 1 } × {1, . . . , 𝐿 2 } × . . . × {1, . . . , 𝐿 𝑙 } is the design space for 𝑙 categorical qualitative variables with their respective 𝐿 1 , ..., 𝐿 𝑙 levels. Our goal is to build a GP surrogate model for 𝑓 . In this case, the GP model will be constructed following the same methodology used for continuous design space (see Eq. ( 3) and Eq. ( 4)). The only changes are related to the construction of the correlation matrix 𝑅. In fact, for a given couple (𝑟, 𝑠) ∈ ({1, . . . , 𝑛 𝑡 }) 2 , let 𝑤 𝑟 = (𝑥 𝑟 , 𝑧 𝑟 , 𝑐 𝑟 ) ∈ R 𝑛 × Z 𝑚 × F 𝑙 and 𝑤 𝑠 = (𝑥 𝑠 , 𝑧 𝑠 , 𝑐 𝑠 ) ∈ R 𝑛 × Z 𝑚 × F 𝑙 two points from the design space. In this case, the correlation kernel [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF] is given by the product of continuous and categorical kernels as:

𝑘 (𝑤 𝑟 , 𝑤 𝑠 , Θ) = 𝑘 𝑐𝑎𝑡 𝑐 𝑟 , 𝑐 𝑠 , Θ 𝑐𝑎𝑡 𝑘 𝑐𝑜𝑛𝑡 (𝑥 𝑟 , 𝑧 𝑟 ), (𝑥 𝑠 , 𝑧 𝑠 ), 𝜃 𝑐𝑜𝑛𝑡 , (7)
where Θ = [Θ 𝑐𝑎𝑡 , 𝜃 𝑐𝑜𝑛𝑡 ], the kernel 𝑘 𝑐𝑜𝑛𝑡 ((., .), (., .), 𝜃 𝑐𝑜𝑛𝑡 ) is constructed efficiently as before (with the continuous relaxation of the integer inputs 𝑧) and the term 𝑘 𝑐𝑎𝑡 (., ., Θ 𝑐𝑎𝑡 ) is a categorical kernel [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF] that depends on a matrix of hyperparameters Θ 𝑐𝑎𝑡 . Using Eq. ( 7), we have

[𝑅] 𝑟 ,𝑠 (Θ) = 𝐶𝑜𝑟𝑟 (𝑤 𝑟 , 𝑤 𝑠 ) = [𝑅] 𝑐𝑎𝑡 𝑟 ,𝑠 (Θ 𝑐𝑎𝑡 ) [𝑅] 𝑐𝑜𝑛𝑡 𝑟 ,𝑠 (𝜃 𝑐𝑜𝑛𝑡 ), (8)
where

𝑅 𝑐𝑜𝑛𝑡 𝑟 ,𝑠 (𝜃 𝑐𝑜𝑛𝑡 ) = 𝑘 𝑐𝑜𝑛𝑡 ((𝑥 𝑟 , 𝑧 𝑟 ), (𝑥 𝑠 , 𝑧 𝑠 ), 𝜃 𝑐𝑜𝑛𝑡 ) and 𝑅 𝑐𝑎𝑡 𝑟 ,𝑠 (Θ 𝑐𝑎𝑡 ) = 𝑘 𝑐𝑎𝑡 (𝑐 𝑟 , 𝑐 𝑠 , Θ 𝑐𝑎𝑡 ).
In the general setting, the categorical kernel 𝑘 𝑐𝑎𝑡 needs to be chosen such that the correlation matrix 𝑅 𝑐𝑎𝑡 is symmetric positive definite (SPD) [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF]. For the categorical inputs, the kernel 𝑘 𝑐𝑎𝑡 (𝑐 𝑟 , 𝑐 𝑠 , Θ 𝑐𝑎𝑡 ) is constructed on the following way. ∀𝑖 ∈ {1, . . . , 𝑙}, 𝑐 𝑟 𝑖 is the level taken by the 𝑖 𝑡 ℎ component of the input 𝑐 𝑟 . As in Pelamatti et al. [START_REF] Pelamatti | Overview and Comparison of Gaussian Process-Based Surrogate Models for Mixed Continuous and Discrete Variables: Application on Aerospace Design Problems[END_REF], let 𝑘 𝑐𝑎𝑡 be formulated level-wise as:

𝑘 𝑐𝑎𝑡 (𝑐 𝑟 , 𝑐 𝑠 , Θ 𝑐𝑎𝑡 ) = 𝑙 𝑖=1 𝐾 𝑖 𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖
where every sub-kernel 𝐾 𝑖 is associated with a correlation matrix 𝑅 𝑖 that contains the correlations between the various levels of the categorical variable 𝑖. Namely, we have

𝑅 𝑐𝑎𝑡 𝑟 ,𝑠 (Θ 𝑐𝑎𝑡 ) = 𝑙 𝑖=1 [𝑅 𝑖 ] 𝑟 ,𝑠 (Θ 𝑖 ).
Thus, the hyperparameters Θ 𝑐𝑎𝑡 can be seen as a concatenation of the set of matrices Θ 1 , . . . , Θ 𝑙 , i.e., Θ 𝑐𝑎𝑡 = [Θ 1 , . . . , Θ 𝑙 ]. The full set of hyperparameters Θ will be estimated using the DoE data set (𝑥, 𝑦 𝑓 ) via an MLE approach on the following way

Θ * = arg max Θ - 1 2
𝑦 𝑓 𝑇 [𝑅(Θ)] -1 𝑦 𝑓 - 1 2 log | [𝑅(Θ)] | - 𝑛 𝑡 2 log 2𝜋 , (9)
where 𝑅(Θ) is computed using Eq. [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF].

IV. Towards a general correlation matrix representation for a Gaussian kernel

We propose a novel approach that tackles the problem of extending correlation kernels to categorical variables by replacing the distance between input points with the kernel that depends only on the hyperparameters. In this section, we will present first the mathematical framework for the Gaussian kernel which fits well with classical dimension reduction techniques such as 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]. The Gaussian kernel has a natural extension to the use of hyperparameters in a matrix form that is an usual form when handling categorical design variables. The proposed extended model will lead to a generalization of both continuous relaxation and Gower distance based methods.

For our purposes, the treatment of continuous inputs will not bring any additional difficulty. Thus, without loss of generality, we will consider only categorical inputs. Hence, in the following, we assume Θ = Θ 𝑐𝑎𝑡 and 𝑅 = 𝑅 𝑐𝑎𝑡 .

A. An extended correlation matrix approach for Gaussian kernel

In the case of Gaussian kernel, a natural extension of hyperparameters to the matrix form can be as follows. Starting from (5) and replacing the vector 𝜃 𝑐𝑜𝑛𝑡 by a given symmetric matrix Θ, we obtain the following kernel:

𝑘 (𝑥 𝑟 , 𝑥 𝑠 , Θ) = 𝑛 𝑗=1 𝑛 𝑗 ′ =1 exp -|𝑥 𝑟 𝑗 -𝑥 𝑠 𝑗 | [Θ] 𝑗, 𝑗 ′ |𝑥 𝑟 𝑗 ′ -𝑥 𝑠 𝑗 ′ | (10)
If Θ is a diagonal matrix with only positive values, then the resulting kernel would be the exact Gaussian kernel for continuous inputs with the hyperparameters being the inverse correlation length. This kernel, given by Eq. ( 10) generalizes the continuous Gaussian kernel for a 2D correlation matrix. Now, for a given 𝑖 ∈ {1, . . . , 𝑙}, let 𝑐 𝑖 be a categorical variable characterized by 𝐿 𝑖 levels. The mapping to a 𝐿 𝑖 -dimensional Hilbert space defined in such a way that the only non-zero coordinate of the image is 1 in the dimension associated to the mapped level is the so-called one-hot encoding [START_REF] Golovin | Google Vizier: A Service for Black-Box Optimization[END_REF]. Let 𝑒 𝑐 𝑖 be the one-hot encoding of 𝑐 𝑖 that takes value 0 everywhere and value 1 on the dimension corresponding to the level taken by the category 𝑐 on the variable 𝑖, 𝑒 𝑐 𝑖 ∈ {0, 1} 𝐿 𝑖 . For example, if 𝑐 takes the 𝑗 𝑡 ℎ level on the variable 𝑖, 𝑒 𝑐 𝑖 𝑗 = 1 and 𝑒 𝑐 𝑖 = [0, . . . , 0, 1, 0, . . . , 0] ∈ R 𝐿 𝑖 . Inspired by Eq. [START_REF] Rasmussen | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF], it is now possible to define a 𝐿 𝑖 × 𝐿 𝑖 symmetric matrix Θ 𝑖 using positive correlation values. This leads to the following general formulation which is the natural extension of Gaussian kernel when dealing with a matrix of hyperparameters:

𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 ) = 𝐿 𝑖 𝑗=1 𝐿 𝑖 𝑗 ′ =1 exp -(𝑒 𝑐 𝑟 𝑖 ) 𝑗 -(𝑒 𝑐 𝑠 𝑖 ) 𝑗 [Θ 𝑖 ] 𝑗, 𝑗 ′ (𝑒 𝑐 𝑟 𝑖 ) 𝑗 ′ -(𝑒 𝑐 𝑠 𝑖 ) 𝑗 ′ (11)
Hence by the definition of 𝑒 𝑐 𝑟 𝑖 and

𝑒 𝑐 𝑠 𝑖 , if 𝑐 𝑟 𝑖 = 𝑐 𝑠 𝑖 , one deduces that 𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑟 𝑖 , Θ 𝑖 ) = exp(0) = 1. Otherwise, if 𝑐 𝑟 𝑖 ≠ 𝑐 𝑠 𝑖 , one gets 𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 ) = exp - 𝐿 𝑖 ∑ 𝑗=1 𝐿 𝑖 ∑ 𝑗 ′ =1 (𝑒 𝑐 𝑟 𝑖 ) 𝑗 -(𝑒 𝑐 𝑠 𝑖 ) 𝑗 [Θ 𝑖 ] 𝑗, 𝑗 ′ (𝑒 𝑐 𝑟 𝑖 ) 𝑗 ′ -(𝑒 𝑐 𝑠 𝑖 ) 𝑗 ′ = exp -[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑟 𝑖 + [Θ 𝑖 ] 𝑐 𝑠 𝑖 ,𝑐 𝑠 𝑖 + [Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖 + [Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖 = exp -[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑟 𝑖 + [Θ 𝑖 ] 𝑐 𝑠 𝑖 ,𝑐 𝑠 𝑖 exp -2[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖 . (12)
where

[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠
𝑖 is the coefficient characterizing the correlation between the two discrete categorical levels taken by 𝑐 𝑟 and 𝑐 𝑠 in the 𝑖 𝑡 ℎ categorical component.

In the following, as far as the matrices Θ 𝑖 respect a specific parameterization, we will show that our approach will guarantee to the correlation matrix 𝑅 to be SPD with unit diagonal and off-diagonal term values in [0, 1] [START_REF] Qian | Gaussian Process Models for Computer Experiments With Qualitative and Quantitative Factors[END_REF]. The latter properties are needed to avoid numerical issues during the computations, see Eq. (3) and Eq. (4). In fact, for all 𝑖 ∈ {1, . . . , 𝑙}, we propose to use the following parameterization for the hyperparameter matrix Θ 𝑖 :

[Θ 𝑖 ] 𝑗, 𝑗 ≥ 0 [Θ 𝑖 ] 𝑗, 𝑗 ′ = log 𝜖 2
( [𝐶 𝑖 𝐶 𝑇 𝑖 ] 𝑗, 𝑗 ′ -1) if 𝑗 ≠ 𝑗 ′ (13)
where for every categorical variable 𝑖, 𝐶 𝑖 is a Cholesky lower triangular matrix that relies on 𝐿 𝑖 (𝐿 𝑖 -1)/2 elements in [0, 𝜋 2 ] that represent the coordinates of a point on the surface of a sphere with a unit radius as in [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF] (see Appendix VII.A for a proof in this context). The parameter 𝜖 is a small tolerance (1 > 𝜖 > 0). Equation ( 13) is chosen such that the elements of 𝑅 𝑖 will be in [0, 1], thanks to the hypersphere decomposition [START_REF] Zhou | A Simple Approach to Emulation for Computer Models With Qualitative and Quantitative Factors[END_REF][START_REF] Rebonato | The Most General Methodology to Create a Valid Correlation Matrix for Risk Management and Option Pricing Purposes[END_REF]. In the following, we will show how this model generalizes all the previous works and proposes a new framework and new matrix configurations.

B. Equivalence with other categorical kernels

The proposed parameterization of the matrix 𝑅 𝑖 using Eq. ( 12) can be seen as the product of the continuous relaxation kernel [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF] and the Gaussian homoscedastic hypersphere kernel, for a total of 𝐿 𝑖 (𝐿 𝑖 +1) 2 hyperparameters per categorical variable. Our proposed approach guarantees that the correlation matrix is SPD and that all its elements are in [0, 1] as in the continuous case. It follows that our proposed kernel is also equivalent to the Gaussian homoscedastic hypersphere model alone.

Note that the hypersphere decomposition, as proposed by Pelamatti et al. [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF], allows negative correlation. Indeed, for categorical variables, values between -1 and 0 are considered because two levels can correspond to opposite effects whereas for continuous GP, the furthest the points, the smaller the correlation. In this case, it follows that our proposed Gaussian homoscedastic hypersphere model

𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 ) = exp -2[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖
is not completely equivalent to the model proposed in [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF] that does not use a positive kernel (like the Gaussian one). A negative correlation would be close to 0 through Gaussian kernel. Indeed, it is equivalent if and only if the correlations are strictly positive, i.e on [0, 𝜋 2 ](Appendix VII.A).

In what comes next, this GP model using our kernel-based approach will be called the homogeneous full model (defined in Eq. ( 12)). It allows to generalize existing approaches in the following way:

• If ∀𝑖 ∈ {1, . . . , 𝑙} the matrices Θ 𝑖 are set to be diagonal with only positive values, then the obtained kernel will correspond to the continuous relaxation method one [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF] with 𝐿 𝑖 hyperparameters per categorical variable, i.e.,

𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 ) = 𝐿 𝑖 𝑗=1 exp -[Θ 𝑖 ] 𝑗, 𝑗 (𝑒 𝑐 𝑟 𝑖 ) 𝑗 -(𝑒 𝑐 𝑠 𝑖 ) 𝑗 2 , ∀𝑐 𝑟 𝑖 ≠ 𝑐 𝑠 𝑖 .
• If ∀𝑖 ∈ {1, . . . , 𝑙} the matrices Θ 𝑖 have all the diagonal terms equal to zero, then, the obtained kernel will be reduced to the Gaussian homoscedastic hypersphere one [START_REF] Pelamatti | Overview and Comparison of Gaussian Process-Based Surrogate Models for Mixed Continuous and Discrete Variables: Application on Aerospace Design Problems[END_REF] with 𝐿 𝑖 (𝐿 𝑖 -1)

2 hyperparameters per categorical variable, i.e., 𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 ) = exp -2[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖 (14)
• If ∀𝑖 ∈ {1, . . . , 𝑙} the matrices Θ 𝑖 have all the off-diagonal terms equal to [Θ 𝑖 ] 𝑐𝑜𝑣 < 0 and zero on the diagonal. In this case, the obtained kernel will be reduced to Gower distance based model [START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF] with only 1 hyperparameter, i.e.,

𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 ) = exp (-2[Θ 𝑖 ] 𝑐𝑜𝑣 ) , ∀𝑐 𝑟 𝑖 ≠ 𝑐 𝑠 𝑖 . ( 15
)
As one can see, our proposed framework generalizes different approaches by considering different matrix structures of the same model. In what follows, the full parameterization is called 𝑚𝑎𝑡_𝐹𝑈 𝐿𝐿, the diagonal representation (which is equivalent to continuous relaxation) will be called 𝑚𝑎𝑡_𝐶 𝑅, the off-diagonal representation will be is called 𝑚𝑎𝑡_𝐺𝐻𝐻, as it is the Gaussian homoscedastic hypersphere model. Last, when we will consider only one covariance, as in Gower distance, we will denote the model 𝑚𝑎𝑡_𝐺𝑂𝑊 𝐸 𝑅. Table 1 summarizes how to obtain the several existing approaches using our proposed framework.

          [Θ 𝑖 ] 1,1 9
Sym.9

[Θ 𝑖 ] 1,2 [Θ 𝑖 ] 2,2 9 . . . . . . . . . 9 [Θ 𝑖 ] 1,𝐿𝑖 . . . [Θ 𝑖 ] 𝐿𝑖-1,𝐿𝑖 [Θ 𝑖 ] 𝐿𝑖 ,𝐿𝑖           exp -[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑟 𝑖 + [Θ 𝑖 ] 𝑐 𝑠 𝑖 ,𝑐 𝑠 𝑖 exp -2[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖 1 2 𝐿 𝑖 (𝐿 𝑖 + 1)
Continuous relaxation [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF] (i.e., Θ 𝑖 = 𝑚𝑎𝑡_𝐶 𝑅)

          [Θ 𝑖 ] 1,1 9 9 Sym. 0 [Θ 𝑖 ] 2,2 9 . . . . . . . . . 9 0 . . . 0 [Θ 𝑖 ] 𝐿𝑖 ,𝐿𝑖           exp -[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑟 𝑖 + [Θ 𝑖 ] 𝑐 𝑠 𝑖 ,𝑐 𝑠 𝑖 𝐿 𝑖
Gaussian homoscedastic hypersphere [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF] (i.e., Θ 𝑖 = 𝑚𝑎𝑡_𝐺𝐻𝐻)

          0 9 Sym.9 [Θ 𝑖 ] 1,2 0 9 . . . . . . . . . 9 [Θ 𝑖 ] 1,𝐿𝑖 . . . [Θ 𝑖 ] 𝐿𝑖-1,𝐿𝑖 0           exp -2[Θ 𝑖 ] 𝑐 𝑟 𝑖 ,𝑐 𝑠 𝑖 1 2 𝐿 𝑖 (𝐿 𝑖 -1)
Gower distance [START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF] (i.e., Θ 𝑖 = 𝑚𝑎𝑡_𝐺𝑂𝑊 𝐸 𝑅)

[Θ 𝑖 ] 𝑐𝑜𝑣           0 9
Sym.9

1 0 9 . . . . . . . . . 9 1 . . . 1 0           exp (-2[Θ 𝑖 ] 𝑐𝑜𝑣 ) 1
Finally, we note that, although our methodology is detailed only for Gaussian kernels, it is possible to extend our proposed approach to other kernels. In the next section, we will see how these models behave on different test cases. in [0, 1] and on a categorical variable with 13 levels. Let 𝑤 = (𝑥, 𝑐) be a given point with 𝑥 being the continuous variable and 𝑐 being the categorical variable, 𝑐 ∈ {1, . . . , 13}. There are two groups of curves corresponding to levels 1 to 9 and levels 10 to 13 with strong within-group correlations, and strong negative between-group correlations (see Appendix VII.D for a detailed description of the function).

The number of relaxed dimensions for continuous relaxation is 14. We draw a 14 × 7=98 points DoE by Latin Hypercube Sampling (LHS) [START_REF] Jin | An efficient algorithm for constructing optimal design of computer experiments[END_REF] and plot the mean Gaussian process models on Fig. 3 for Gower distance, continuous relaxation and Gaussian homoscedastic hypersphere. The number of hyperparameters to optimize is therefore 2 for Gower, 14 for continuous relaxation and 79 for Gaussian homoscedastic hypersphere. Then, from these models, we compute the root mean square error (RMSE) as RMSE =

√ 𝑛 𝑖=1 1 𝑛 f (𝑤 𝑖 ) -𝑓 (𝑤 𝑖 )
where 𝑛 is the size of the validation set, f (𝑤 𝑖 ) is the prediction of our model at 𝑤 𝑖 and 𝑓 (𝑤 𝑖 ) is the true value.

As expected, the more the hyperparameters, the better the model, as it can be seen on the decreasing RMSE. However, Gower distance takes 1.4 seconds to compute, continuous relaxation takes 24.5 seconds and Gaussian homoscedastic hypersphere takes 514.5 seconds to compute which motivates the use of a reduced order model. We can also consider the full homogeneous model with 92 hyperparameters; this model takes 642 seconds to compute and is worse than the Gaussian homoscedastic hypersphere. This model is not to consider for practical use cases. The estimated correlation matrix 𝑅 𝑖 = 𝑅 1 is shown in Fig. 4. For two given levels {𝑟, 𝑠}, the correlation [𝑅 1 ] 𝑟 ,𝑠 is in blue if the correlation is close to 1 and in red if the correlation is close to 0. We can see on the figure that the correlation between a level and itself is always 1. For Gower distance, there is only one estimated "mean correlation" as in Fig. 4a. For continuous relaxation (see Fig. 4b), we have

[𝑅 1 ] 𝑟 ,𝑠 = exp(-[Θ 1 ] 𝑟 ,𝑟 + [Θ 1 ] 𝑠,𝑠
), therefore for the most important levels (1 to 9) are strongly correlated (in blue) with one another and the other levels (10 to 13) that should also have been also correlated are badly estimated because of the model limitation that neglected them. In contrast, the Gaussian parametrization (see Fig. 4c) of the Gaussian homoscedastic hypersphere decomposition model gives a good approximation of the real correlation and we see that there are two groups of highly correlated levels. The levels 1 to 9 are strongly similar, the levels 10 to 13 are strongly similar and the two groups are less similar. The latter correlations should have been different but we do not allow negative values through the Gaussian kernel.

As previously mentioned, the full homogeneous model (see Fig. 4d) that combines both continuous relaxation and Gaussian homoscedastic hypersphere adds irrelevant parameters making it more hard to optimize numerically while being equivalent to the Gaussian homoscedastic hypersphere model. For this reason, this model should not be considered for real applications. We compare in Fig. 5 our Gaussian hypersphere decompositon model (see Fig. 5a) with the homoscedastic hypersphere model [START_REF] Pelamatti | Overview and Comparison of Gaussian Process-Based Surrogate Models for Mixed Continuous and Discrete Variables: Application on Aerospace Design Problems[END_REF] that allows negative values (see Fig. 5b). By imposing the correlation to be the matrix correlation of Fig. 5b in the case of Homoscedastic hypersphere, we obtain a likelihood of around 210 against 138 for the Gaussian homoscedastic hypersphere model with only positive values. Nevertheless, even if the Homoscedastic hyperphere model is more general than our model, the two RMSE are of the same order of magnitude, indicating similar performances on that particular test case.

Categorical Branin function

Let the function 𝑓 to model be the modified categorical Branin function [START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF]. This problem has 3 variables: two continuous variables in [0, 1] and one categorical variable with three levels and the two first levels are totally correlated. We draw a DoE of 60 points by LHS to compare the given GP models and compute the error terms. To begin with, we start by plotting the models built with the different methods in Fig. 6 and compute their respective RMSE to compare them with the original formulation of the methods. The obtained values are given in Fig. 6 from a validation base of size 30603 that corresponds to 101 points from 0 to 1 in every continuous direction for every level (see Appendix VII.C for a detailed description of the function). In this test case, we clearly show that the full model, the continuous relaxation model and the Gaussian homoscedastic one are the same. However, we still have numerical instabilities but using the Gaussian homoscedastic kernel as proposed in this paper instead of using the raw matrix leads to a better estimation of the hyperparameters. For homoscedastic hypersphere, the RMSE that we found is 63.021. The Gower distance model is the only one that differs visually and that differs consequently in error from the others (77.8 instead of 60.7). As mentioned in Section V.A, with 3 levels, Gaussian homoscedastic hypersphere and continuous relaxation are equivalent methods, so these results were theoretically expected. However, when comparing with the original homoscedastic hypersphere, the Gaussian model differs as the third level is negatively correlated from the two others and because the Gaussian kernel returns only positive value. In this case, the estimated correlation is 0.03 with the Gaussian kernel against -0.15 with the original one.

VI. Conclusion and perspectives

In this work, we have proposed a kernel-based class of models that extends continuous Gaussian kernels to handle mixed-categorical variables. We have observed on test cases that our proposed models generalizes Gower distance and continuous relaxation based models. Numerical illustrations on analytical problems showed the good potential of the proposed models. Further works will consider to include dimension reduction techniques (such as KPLS) to improve the computational efficiency of our model and be able to tackle higher dimensional problems. A rigorous classification between the proposed models as well as a proof of the SPD nature of the correlation matrices within our proposed model shall be investigated.

Then, we take 𝐶 = 𝐵𝐵 𝑇 as a Cholesky decomposition and this is the so-called "hypersphere decomposition" [START_REF] Zhou | A Simple Approach to Emulation for Computer Models With Qualitative and Quantitative Factors[END_REF]. By restricting our coordinates 𝜃 to be in [0, 𝜋 2 ], we still have a bijection, by construction. This bijection takes values between 𝜃 = [0, 𝜋 2 ] and 𝜉 = 𝐹 ( [0, 𝜋 2 ]) = [0, 1]. This is the so-called first quadrant as in Fig. 7. For a given 𝜉 * > 𝜖 that corresponds to the optimal correlation according to the maximum of likelihood estimator, with the homoscedastic hypersphere method, ∃! 𝜃 * ∈ [0, 𝜋 2 ] : 𝐹 (𝜃 * ) = 𝜉 * . For the same 𝜉 * > 𝜖 that corresponds to the optimal correlation, with the Gaussian homoscedastic hypersphere method, ∃! 𝜃 ′ * ∈ [0, 𝜋 2 ] : 𝐹 (𝑔(𝜃 ′ * )) = 𝜉 * .

Therefore, the two methods are equivalent over [𝜖, 1]. By choosing 𝜖 sufficiently small, the two methods are equivalent to model positive correlations. Note that we obtain the same optimal correlation and the same likelihood but the values of 𝜃 * and 𝜃 ′ * may differ even if being uniquely defined on [0, 𝜋 2 ].

B. 2D blue/red test case

This test case has one categorical variable with two levels: blue or red and one continuous variable in [0, 4]. The blue DoE is the following: x= {0, 1, 4}, y={0, 9, 16} The red DoE is the following: x= {0, 1, 2, 3}, y={0, 1, 8, 27} Therefore, we have a DoE consisting of 7 points either blue or red, with continuous value ranging between 0 and 4 and taking value between 0 and 27.

C. Categorical Branin function case

This test case has one categorical variable with three levels and two continuous variables in [0, 1]. Let 𝑤 = (𝑥 1 , 𝑥 2 , 𝑐) be a given point with 𝑥 1 and 𝑥 2 being the continuous variables and 𝑐 being the categorical variable, 𝑐 ∈ {0, 1, 2}. The DoE is given by a LHS of 60 points. Our validation set is a evenly spaced grid of 101 points in 𝑥 1 ranging from 0.01 to 0.99, 101 points in 𝑥 2 ranging from 0.01 to 0.99 for every of the three categorical levels (0, 1, 2) for a total of 30603 points.

𝑓 (𝑤

D. Categorical cosine case

This test case has one categorical variable with 13 levels and one continuous variable in The DoE is given by a LHS of 98 points. Our validation set is a evenly spaced grid of 101 points in 𝑥 1 ranging from 0.01 to 0.99, 101 points in 𝑥 2 ranging from 0.01 to 0.99 for every of the three categorical levels (0, 1, 2) for a total of 30603 points.

Fig. 1

 1 Fig. 1 Drag polar and aerodynamic properties for an efficient supersonic air vehicle obtained from a surrogate model for different Mach speed and sweep angles [7, Figure 3].

( a )

 a Our model (with Θ 1 = 𝑚𝑎𝑡_𝐺𝑂𝑊 𝐸 𝑅): 2 hyperparameters, CPU time =1.4s, RMSE=30.079. (b) Our model (with Θ 1 = 𝑚𝑎𝑡_𝐶 𝑅): 14 hyperparameters, CPU time =24.5s, RMSE=22.347. (c) Our model (with Θ 1 = 𝑚𝑎𝑡_𝐺𝐻𝐻): 79 hyperparameters, CPU time =514.5s, RMSE=2.941. (d) Our model (with Θ 1 = 𝑚𝑎𝑡_𝐹𝑈 𝐿𝐿): 92 hyperparameters, CPU time =642s, RMSE=22.61.

Fig. 3

 3 Fig. 3 Mean predictions for the cosine problem with a 98 points DoE for our 4 proposed models.

  (a) The matrix 𝑅 𝑐𝑎𝑡 with Θ 1 = 𝑚𝑎𝑡_𝐺𝑂𝑊 𝐸 𝑅. (b) The matrix 𝑅 𝑐𝑎𝑡 with Θ 1 = 𝑚𝑎𝑡_𝐶 𝑅. (c) The matrix 𝑅 𝑐𝑎𝑡 with Θ 1 = 𝑚𝑎𝑡_𝐺𝐻𝐻. (d) The matrix 𝑅 𝑐𝑎𝑡 with Θ 1 = 𝑚𝑎𝑡_𝐹𝑈 𝐿𝐿.

Fig. 4

 4 Fig. 4 Obtained correlation matrices for the cosine problem using a DoE of 98 points.

( a )Fig. 5

 a5 Fig. 5 Comparison results between our model and Homoscedastic hypersphere [27] on the cosine problem using a DoE of 98 points.

( a )Fig. 6

 a6 Fig. 6 Mean predictions (over the three levels) for the categorical Branin problem using a DoE of 60 points (in red in the curve plots).

Fig. 7 2 𝑛: 2 𝑛:

 722 Fig. 7 First quadrant on the sphere S 1 .

15 𝑐 - 𝑐 20 , 2 𝑥 - 𝑐 20 ,

 1520220 [0, 1]. Let 𝑤 = (𝑥, 𝑐) be a given point with 𝑥 being the continuous variable and 𝑐 being the categorical variable, 𝑐 ∈ {1, . . . , 13}. 𝑓 (𝑤) = cos 7𝜋 2 𝑥 + 0.4𝜋 + 𝜋 if c ∈ {10, . . . , 9} 𝑓 (𝑤) = cos 7𝜋 if c ∈ {10, . . . , 13} The reference landscapes of the objective function (with respect to the categorical choices) are drawn on Fig. 8.

Fig. 8

 8 Fig. 8 Landscape of the cosine test case from [26].

Table 1 A summary of how our proposed kernel-based approach generalizes existing categorical models.

 1

	Models	Θ 𝑖	𝐾 𝑖 (𝑐 𝑟 𝑖 , 𝑐 𝑠 𝑖 , Θ 𝑖 )	# of Hyperparam.
	Our full model
	(i.e., Θ 𝑖 = 𝑚𝑎𝑡_𝐹𝑈 𝐿𝐿)

Acknowledgements

This work is part of the activities of ONERA -ISAE -ENAC joint research group. The research presented in this paper has been performed in the framework of the AGILE 4.0 project (Towards Cyber-physical Collaborative Aircraft Development) and has received funding from the European Union Horizon 2020 Programme under grant agreement n • 815122. The authors are grateful to the partners of the AGILE 4.0 consortium for their contribution and feedback.

V. Results

In this section, we propose several illustrations and comparisons on three different test cases (with 1 or 2 continuous variables, 1 categorical variable up to 13 levels) to show the interest of our method and the equivalence with other models from the literature. The optimization of the likelihood as a function of the hyperparameters needs a performing algorithm, in this work, we are using COBYLA [START_REF] Powell | A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation[END_REF] to maximize this quantity. We note that the implementation of our proposed method has been released in the toolbox SMT v1.2 * [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF] and further developments are to appear in the next release.

A. Hyperparameters equivalence

We start with a 2D test case with one continuous variable in [0, 4] and one categorical variable with two levels (blue or red). We consider a DoE of 3 blue points and 4 red points (see Appendix VII.B for a detailed description of the test case). The mixed integer GP model is shown in Fig. 2 with the two associated levels (blue and red) and is obtained from continuous relaxation. As any other method leads to the same GP model of Fig. 2, we only report the optimal values of the hyperparameters in Tab. 2 to illustrate this equivalence. In this case, with only two levels, it is easy to prove that the Gower distance kernel and the homoscedastic hypersphere one are equivalent: there is only one categorical hyperparameter and we can check in Tab. 2 that exp (-0.23015) = 0.7944. On this particular test case, all the correlations are positive and so, the Gaussian homoscedastic hypersphere and the original homoscedastic hypersphere models are equivalent as they can be restricted to angles in [0, 𝜋 2 ] (see Appendix VII.A for a proof). A negative correlation would have been model by a close to 0 correlation through Gaussian kernel.

VII. Appendix

A. Positive homoscedastic hypersphere Theorem 1. We have equivalence between homoscedastic hypersphere [START_REF] Pelamatti | Mixed-variable Bayesian optimization : application to aerospace system design[END_REF] and Gaussian homoscedastic hypersphere if and only if the correlations are positive (> 𝜖) or by restricting the angles to be on [0, 𝜋 2 ] that is the first quadrant. Proof. Let 𝐹 be the polyspherical change of coordinates [START_REF] Klimyk | Representations of Lie groups and special functions[END_REF], such that 𝐹 (𝜃 1 , . . . , 𝜃 𝑛-1 ) = 𝜉 1 , . . . , 𝜉 𝑛 with 𝜉 = (𝜉 1 , . . . , 𝜉 𝑛 ) ∈ S 𝑛-1 and 𝜃 = (𝜃 1 , . . . , 𝜃 𝑛-1 ) ∈ [0, 2𝜋] 𝑛-1 , for a given integer 𝑛.

The formula for the so-called "polyspherical change of coordinates" given by Vilenkin et al. [START_REF] Klimyk | Representations of Lie groups and special functions[END_REF], restricted to real values, is:

The function 𝐹 is bijective as a change of coordinates [START_REF] Klimyk | Representations of Lie groups and special functions[END_REF]. Correlations are values in [-1, 1] so we consider that we have "unit radius" [START_REF] Rebonato | The Most General Methodology to Create a Valid Correlation Matrix for Risk Management and Option Pricing Purposes[END_REF] (𝑟 = 1). 𝜉 are the so-called "cartesian coordinates" and 𝜃 are the so-called "spherical coordinates".

This leads to the formulation of Rebonato [START_REF] Rebonato | The Most General Methodology to Create a Valid Correlation Matrix for Risk Management and Option Pricing Purposes[END_REF] for a 𝑚 × 𝑚 matrix. Let 𝑏 𝑖 𝑗 be the elements of a triangular superior matrix 𝐵 of size 𝑚 × 𝑚. The 𝑚 × (𝑚 -1) corresponding 𝑏 𝑖 𝑗 elements are determined from the 𝑚 × (𝑚 -1) 𝜃 𝑖 𝑗 as: