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

Eric Nguyen Van
email: eric.nguyen_van@onera.fr

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Thierry Lefebvre
email: thierry.lefebvre@onera.fr

Christophe David
email: christophe.david@onera.fr

Sébastien Defoort
email: sebastien.defoort@onera.fr

Youssef Diouane
email: youssef.diouane@isae-supaero.fr

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

Bayesian optimization for mixed variables using an adaptive dimension reduction process: applications to aircraft design

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. Introduction

New 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]. 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]. The process of finding the best configuration, known as Multidisciplinary Design Optimization (MDO), is formulated as a minimization constrained problem where the objective and the constraint functions are typically expensive-to-evaluate and their derivatives are not available.

In the context of aircraft design, the MDO 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. Thus, the regarded optimization problem is of the following form: min =( , , ) ∈Ω× ×F { ( ) s.t. ( ) ≤ 0 and ℎ( ) = 0}

: R × Z × F ↦ → R is the objective function, : R × Z × F ↦ → R gives the inequality constraints, and ℎ : R × Z × F ↦ → R returns the equality constraints. Ω ⊂ 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. The functions , , and ℎ are typically expensive-to-evaluate simulations with no exploitable derivative information.

When only continuous design variables are optimized (i.e., F is reduced to a single categorical choice), Bayesian optimization (BO) is shown to be powerful strategy for solving problem [START_REF] Duriez | HALE multidisciplinary design optimization with a focus on eco-material selection[END_REF] [START_REF] Mockus | Bayesian approach to global optimization: theory and applications[END_REF]. BO uses Gaussian processes (GPs) [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 denoted by Kriging [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF]) to define response surface models, and the sequential enrichment is performed by maximizing a given acquisition function [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]. The latter is meant to model a compromise between exploration of new zones in the design space and exploitation (i.e. minimization) of the GPs. For general mixed integer problems, several modeling strategies to build GPs 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]. Based on these GP models, a trade-off acquisition function was adapted for optimization [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF]. Some other modelling strategies consist in computing a continuous model for each category [START_REF] Roy | A Mixed Integer Efficient Global Optimization Algorithm with Multiple Infill Strategy -Applied to a Wing Topology Optimization Problem[END_REF], either by continuously relaxing the design variables [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF], by using a multi-armed bandit strategy to handle the categorical choices [START_REF] Nguyen | Bayesian Optimization for Categorical and Category-Specific Continuous Inputs[END_REF] or by considering a Gower distance to model simultaneously the proximity over categorical and continuous variables [START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF]. Recently, a continuous relaxation BO based method [START_REF] Karlsson | Continuous surrogate-based optimization algorithms are well-suited for expensive discrete problems[END_REF] to tackle mixed integer variables has been shown to solve efficiently expensive-to-evaluate optimization problems. In fact, using continuous relaxation within BO leads to better results. However, the relaxation of the categorical design variables increases the number of the hyperparameters needed (to be tuned) associated with the GP model. This in particular constrained the method in [START_REF] Karlsson | Continuous surrogate-based optimization algorithms are well-suited for expensive discrete problems[END_REF] to be used only for small dimensional optimization problems. Since, the construction of the GP model may not be scalable to practical applications involving a large number of mixed variables, some specific strategies have to be investigated to handle this complexity: the Kriging model with Partial Least Squares (KPLS) is one of the most used technique [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by Partial Least Squares dimension reduction[END_REF][START_REF] Bouhlel | An Improved Approach for Estimating the Hyperparameters of the Kriging Model for High-Dimensional Problems through the Partial Least Squares Method[END_REF].

In this work, we show how to reduce the computational cost related to the construction of the mixed categorical GP model as proposed in [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF]. Our proposed approach allows in particular to solve higher dimensional mixed integer MDO problems. The method relies on combining continuous relaxation and the use of KPLS. Inspired by the Wold's R criterion [START_REF] Li | Model selection for partial least squares regression[END_REF][START_REF] Tarpey | A Note on the Prediction Sum of Squares Statistic for Restricted Least Squares[END_REF][START_REF] Wold | Cross-Validatory Estimation of the Number of Components in Factor and Principal Components Models[END_REF], we propose also an adaptive mechanism for choosing the number of hyperparameter to use within KPLS. The good potential of the proposed approach is first showed over a set of analytical test cases.

The performance of the proposed approach is then confirmed on two MDO applications from the FAST-OAD framework [START_REF] David | From FAST to FAST-OAD: An open source framework for rapid Overall Aircraft Design[END_REF]: "CERAS" and "DRAGON". FAST-OAD is an open-source Python framework that provides a flexible way to build and solve the Overall Aircraft Design problems by assembling discipline models from various sources: it currently comes with some bundled, quick and simple, models dedicated to commercial aircraft. In this paper, FAST-OAD will resolve an MDA problem that mainly: (a) sizes the geometry of main aircraft components, (b) computes mass and centers of gravity of aircraft parts, (c) estimates the aerodynamics and propulsion along the computed mission, and (d) returns the fuel consumption related to the mission. These estimated quantities will be used to define the objective and the constraints of our two optimization problems.

The outline of the paper is as follows. In Section II, a detailed review of the Bayesian optimization framework is given. The continuous relaxation as well as the use of the KPLS technique are detailed in Section III. Section IV presents our academical tests as well as the obtained results on the two regarded MDO test cases. Conclusions and perspectives are finally drawn in Section V.

II. Continuous constrained Bayesian optimization

In this section, we will only consider that all the design variables are continuous in problem [START_REF] Duriez | HALE multidisciplinary design optimization with a focus on eco-material selection[END_REF]. Namely, in this section, the design space will be restricted to Ω ⊂ R ; hence = in the optimization problem [START_REF] Duriez | HALE multidisciplinary design optimization with a focus on eco-material selection[END_REF]. In the context of unconstrained optimization, the Efficient Global Optimization (EGO) framework [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] suggests to use the GP model to find the global minimum of an expensive-to-evaluate black-box function (based on the preliminary results of Močkus [START_REF] Močkus | On bayesian methods for seeking the extremum[END_REF]). In general, a GP 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. Thereafter, we will consider that our unknown black-box objective function is a realization of an underlying distribution of mean ˆ and of standard deviation such that (.) ∼ ( ˆ (.), [ ] 2 (.)). Let be 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. The stochastic model [START_REF] Duvenaud | Automatic model construction with Gaussian processes[END_REF] writes as: ( ) = ( ) + ∈ R with the error term between and the model approximation ( ). The errors terms are considered iid of variance 2 . Let be the error correlation matrix between the inputs points = ( ( ), ( )). 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 ( * ) = ( ( * ), ( )) for a given * , then

ˆ ( * ) = ˆ + ( * ) -1 (y -1 ˆ ), (2)
and

[ ] 2 ( * ) = [ ˆ ] 2 1 -( * ) -1 ( * ) + (1 -1 -1 ( * )) 2 1 -1 1 , (3)
where ˆ and ˆ , respectively, are the MLE of and with respect to given the data set (x, y ) [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF].

Within EGO, at a given iteration , a GP surrogate model, referred as Kriging, is computed based on the current DoE to approximate the black-box . Henceforth, one wants to estimate the best new point to evaluate, as it is costly, by taking into account the model information to converge as fast as possible to the real optimum of the black-box. The point that we will evaluate next is the one that gives the best improvement a priori according to an acquisition function like the Expected Improvement (EI) defined over the model. The objective value at this new point will then be evaluated and incorporated into the next surrogate model. The hyperparameters that characterize and define the model are thus updated at each iteration until convergence. The Bayesian optimization process is made from these GPs in an iterative manner.

To tackle constrained Bayesian optimization, EGO was extended to Super-Efficient Global Optimization (SEGO) method [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF]. SEGO uses surrogate models of the constraints to give an estimation of the search space Ω through a given criterion in which the function is optimized. The latter was enhanced to tackle multi-modal and equality constraints with the Upper Trust Bound (UTB) criterion [START_REF] Priem | An efficient application of Bayesian optimization to an industrial MDO framework for aircraft design[END_REF][START_REF] Priem | Upper trust bound feasibility criterion for mixed constrained Bayesian optimization with application to aircraft design[END_REF][START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF]. The acquisition function that we use is the WB2s (Watson and Barnes 2 criterion with scaling) [START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF] that is known to be more robust than the Expected improvement (EI) criterion, especially in high dimension [START_REF] Jesus | Surrogate based MDO of a canard configuration aircraft[END_REF]. WB2s can be seen as smooth version compared to the WB2 [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF] criterion and is less multimodal compared to EI. Algorithm 1 details the SEGO optimization procedure.

Algorithm 1: SEGO for continuous inputs.

Result: Solution of the problem (1) over the continuous design space Ω. Inputs: Initial DoE D 0 and set = 0. Search space Ω. while the stopping criterion is not satisfied do 1) Build the surrogate model of the objective function to obtain the mean and standard deviation prediction at a given point: ( ˆ , ) from the DoE D . 2) Build the surrogate models for every constraint ( ˆ , ), ( ĥ , ℎ ) from the DoE D to compute an estimation of the search space Ω .

3) Construct the acquisition function 2 (.) = ˆ (.), (.) from the objective model.

end

III. Mixed categorical constrained Bayesian optimization

To handle mixed categorical design variables, we propose to use the continuous relaxation method that has been recently shown to be well-suited for expensive discrete problems [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF][START_REF] Karlsson | Continuous surrogate-based optimization algorithms are well-suited for expensive discrete problems[END_REF]. The main drawback of such method is that it enlarges the dimension of the design space according to the size of the categorical space. To overcome such issue, we propose to combine continuous relaxation with the well-known partial least squares procedure [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by Partial Least Squares dimension reduction[END_REF] to reduce the number of the GP hyperparameters. For mixed categorical design variables, the proposed treatment relies on continuous relaxation. The design space Ω × × F is relaxed to a continuous space Ω constructed on the following way:

• ∀ ∈ {1, . . . , }, the variable is relaxed within its bounds and treated as continuous.

• ∀ ∈ {1, . . . , }, we use a relaxed one-hot encoding [START_REF] Golovin | Google Vizier: A Service for Black-Box Optimization[END_REF] for and add new continuous dimensions into Ω . Therefore, we get, after relaxation, a new design space Ω ⊆ R where = + + =1 > + + . The nature of the variables should be respected when evaluating a point in the relaxed space so we define the inverse operator that projects a point ∈ Ω to its closer point in Ω × × F . Namely, rounds the value of an integer variable to the closer value among its levels and, for a categorical variable , selects the level which corresponding dimension value is the highest.

In this work, when building the Kriging model, the error correlation will be estimated using a squared exponential (or Gaussian) correlation kernel over the relaxed design space. We denote by = ( ) the relaxation in Ω of a point ∈ Ω × × F . The mixed categorical kernel is

( , , ) = =1 exp - - 2 , ∈ R + (4)
This kernel relies on hyperparameters estimated by maximum of likelihood such that the more the number of variables for the problem, the more the number of hyperparameters to optimize. Reducing this number leads to a better estimation for the hyperparameters, a more convenient optimization of the likelihood and makes the model scalable for high-dimensional problems. To do so, the Partial Least Squares (PLS) method [START_REF] Wold | Soft Modelling by Latent Variables: The Non-Linear Iterative Partial Least Squares (NIPALS) Approach[END_REF] searches the direction that maximizes the variance between the input and output variables. This is done by a projection into a smaller space spanned by the so-called principal components. The number of principal components that corresponds to the new number of hyperparameters for KPLS is chosen to be much lower than . The resulting PLS squared exponential kernel is given by

( , , ˆ ) = =1 =1 exp -ˆ * - * 2 , ∀ ∈ R + (5)
where [ * ] , are scalars that measure the influence of the input variables on the output . Combining this model construction with SEGO gives the method described in Algorithm 2.

Algorithm 2: SEGO for mixed categorical inputs using an adaptive KPLS.

Result

end

KPLS is an efficient method that can tackle high-dimensional problems by reducing the number of effective dimensions to a small number and it is currently used in the community [START_REF] Bhosekar | Advances in surrogate based modeling, feasibility analysis, and optimization: A review[END_REF][START_REF] Zuhal | On dimensionality reduction via partial least squares for Kriging-based reliability analysis with active learning[END_REF][START_REF] Li | An Improved High-Dimensional Kriging Surrogate Modeling Method through Principal Component Dimension Reduction[END_REF]. The number of principal components being a key point that must be chosen for the PLS technique, we propose a strategy to choose this number in an adaptive process through the optimization iterations. As a rule of the thumb, taking 2 to 5 active components is, in general, efficient for most problems. However, when we do not have any a-priori information, an intuitive idea is to learn from GP the number of active dimension. Bouhlel et al. [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by Partial Least Squares dimension reduction[END_REF] proposed to use the leave-one-out strategy to find the best number of hyperparameters for KPLS. However, this strategy is not efficient because at every iteration, it implies to compute a large number of reduced size surrogate models to select the best one. As the size of the DoE increases, both the model computation and the leave-one-out criterion become more and more costly. In [START_REF] Li | Model selection for partial least squares regression[END_REF], they propose to use the adjusted Wold's R criterion [START_REF] Wold | Cross-Validatory Estimation of the Number of Components in Factor and Principal Components Models[END_REF] for dimensionality reduction purposes. Wold's R criterion is based on cross-validation over -folds [START_REF] Li | Model selection for partial least squares regression[END_REF] and consists in measuring the ratio of the PRedicted Error Sum of Squares (PRESS) [START_REF] Tarpey | A Note on the Prediction Sum of Squares Statistic for Restricted Least Squares[END_REF]. If for a given number of components , the criterion defined as

( ) = PRESS( + 1) PRESS( ) (6)
and should be smaller than a given threshold. The ratio is typically 1, which means we add components until the approximate error stops decreasing. However a threshold of 0.9 or 0.95 could be used in order to be more selective and add fewer components.

Initialization:

• For all subsamples ∈ {1, . . . , }, build a KPLS model with min components over all data except the inputs in the subsample .

• For all ∈ 1, . . . , , let ˆ ,-be the prediction of the KPLS model build without the point in the ℎ fold at the corresponding . This fold is therefore used as a validation set to compute the -fold cross validation PRESS as

PRESS( min ) = =1 =1 -ˆ ,- 2 .
• Set Ψ ←min . while Ψ < max do 1) Divide the DoE D into subsamples.

2) For all subsamples ∈ {1, . . . , }, build a KPLS model with Ψ components over all data except the inputs in the subsample . 3) For all ∈ 1, . . . , , let ˆ ,-be the prediction of the KPLS model build without the point in the ℎ fold at the corresponding . This fold is therefore used as a validation set to compute the -fold cross validation PRESS as PRESS(

Ψ + 1) = =1 =1 -ˆ ,- 2 . Set ( ) = PRESS( Ψ +1) PRESS( Ψ ) . 4) if ( Ψ ) ≥ then STOP; 5) Increment Ψ , i.e., Ψ ← -Ψ + 1 .
end Return KPLS model with Ψ components; Note that when min = max , Algo. 3 is equivalent to KPLS with a constant number of dimensions. The code implementation of our proposed method has been released in the toolbox SMT v1.1 * [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF].

IV. Results and discussion

In this section, we carry out experiments for unconstrained and constrained test cases, with several number of variables and an increasing complexity. We optimize analytical test cases as a benchmark study and then we solve aircraft design optimization problems as a validation and application of the present work.

A. Implementation choices

In order to compare, we used several optimization algorithms detailed hereafter: Bandit-BO, NSGA2, SEGO with Kriging, SEGO with Gower distance and the proposed method, SEGO coupling Kriging and PLS.

The Bandit-BO implementation used is the one by Nguyen et al. [START_REF] Nguyen | Bayesian Optimization for Categorical and Category-Specific Continuous Inputs[END_REF], we are not considering parallelization or batch evaluations. The NSGA2 [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] algorithm used is the implementation from the toolbox pymoo [START_REF] Blank | Pymoo: Multi-Objective Optimization in Python[END_REF] with the default parameters (probability of crossover = 1, eta = 3). Fronts are not relevant in our study as we are considering single-objective optimization. The optimization with SEGO is made from SEGOMOE [START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF] for both constrained and unconstrained cases. For SEGO using Gower distance [START_REF] Halstrup | Black-Box Optimization of Mixed Discrete-Continuous Optimization Problems[END_REF] (denoted by SEGO+GD), we are considering the implementation of the Surrogate Modeling Toolbox (SMT) [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF], an open-source python toolbox where some computations associated to the present work have been done and where automatic PLS was implemented. The same holds * https://smt.readthedocs.io/en/latest/ for SEGO and Kriging coupled with PLS to reduce to the number of hyperparameters (denoted by SEGO+KPLS(d= ). As the PLS could potentially lead to numerical instabilities, we are using the homoscedastic noise that maximizes the likelihood as a so-called nugget. For automatic PLS (denoted by SEGO+KPLS(auto), the -fold cross-validation uses = 4, the experiments [START_REF] Li | Model selection for partial least squares regression[END_REF] suggest to use 4 to 6 fold cross-validation.

For SEGO using Kriging (denoted by SEGO+KRG), we also use the implementation from the toolbox SMT. For the constrained analytical test cases, we are using the criterion [START_REF] Priem | Upper trust bound feasibility criterion for mixed constrained Bayesian optimization with application to aircraft design[END_REF]. Some adaptions have been done to Bandit-BO and NSGA2 to consider both integer and categorical variables. As NSGA2 can only consider integer variables, categorical variables are treated as integer ones. Contrarily, Bandit-BO can treat only categorical variables so integer variables are treated as categorical ones. Bandit-BO creates a GP model for each arm, so it requires at least 2 × initial points, being the number of categorical possibilities for the problem inputs. For unconstrained test cases, these 2 points by categorical possibility are sampled randomly. If we are not using Bandit-BO, for constrained optimization, we are using a continuous relaxed Latin hypercube sampling and then we project the output points to obtain the mixed integer DoE.

For Kriging, Kriging with PLS or Gower distance, the hyperparameters are optimized with COBYLA [START_REF] Powell | A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation[END_REF] and the chosen model regression is constant. When optimizing with SEGO, the acquisition function is maximized using ISRES [START_REF] Runarsson | Search biases in constrained evolutionary optimization[END_REF] to find some interesting starting points and SNOPT [START_REF] Gill | SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization[END_REF] to finalize the process based on these starting points. The squared exponential kernel is the only one considered for these methods.

In order to compute some statistical data (median and variance), we are doing 20 repetitions of the optimization process for a given method and an initial DoE size. We consider that a constraint is respected if the constraint violation is smaller than the threshold value 10 -4 .

B. Benchmark test cases including unconstrained and constrained optimization problems

In this section, a benchmark of different problems is proposed in order to compare the efficiency of the proposed algorithm with some state-of-the-art methods. The first ten are analytical cases, with or without constraints in order to provide some data profiles and the last two concern some more complex applications for aircraft conceptual design.

Unconstrained optimization

To begin with, we validate our method on unconstrained problems up to 14 dimensions. The first analytical test case is a modified Branin function [START_REF] Roy | A Mixed Integer Efficient Global Optimization Algorithm with Multiple Infill Strategy -Applied to a Wing Topology Optimization Problem[END_REF], denoted by "Branin 5", where the first variable is an integer 1 ∈ {-5, -4, . . . , 9, 10} and the second one 2 ∈ [0, 10] is a continuous variable. As this problem is only 2-dimensional, SEGO-KRG is considered without the coupling with the PLS technique. For Bandit-BO, we represent 1 as a bandit with 16 arms associated to the 16 integer values from 0 to 15. Therefore, an initial DoE of at least 32 points is required for Bandit-BO. A smaller initial DoE is also considered for NSGA2 and SEGO-KRG in order to compare the convergence according to the DoE size: 5 points or 32 points. For a given DoE size (5 or 32), 20 different initial DoE are obtained via Latin hypercube sampling in order to compute some statistical data (median and variance) about the convergence results. In Fig. 1a, the medians and the associated quartiles (25% and 75%) of the 20 runs are illustrated for Bandit-BO, NSGA2 and SEGO with Kriging. The initial DoE is shown before the black dotted line. For the DoE with 32 initial points, we are doing 50 iterations of the methods Bandit-BO and SEGO-KRG, for a total of 82 evaluations. For NSGA2, 200 iterations are performed. When the DoE size is reduced to 5 points, 200 iterations of SEGO-KRG are made and compared to 300 evaluations for NSGA2. These comparisons show that SEGO-KRG outperforms BANDIT-BO and NSGA2 in terms of convergence. Moreover, we find that a smaller initial DoE performs better for SEGO, that is a known result [START_REF] Riche | Revisiting Bayesian Optimization in the light of the COCO benchmark[END_REF].

To compare the dispersion of the best results over the 20 runs, Fig. 1b shows the boxplots for SEGO-KRG, Bandit-BO and NSGA2 after 50 iterations associated to a initial DoE of 32 points. These plots show that, not only SEGO-KRG converges better in median but also it is the only method to converge for every of the 20 initial DoE. However, Bandit-BO converges better than NSGA2 on this test case. The second analytical test case is a toy function that consists in the choice between a set of 10 1-dimensional continuous functions [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF], denoted by "Set 1", where the first variable is a categorical variable 1 ∈ {0, 1, . . . , 9} and the second one 2 ∈ [0, 1] is a continuous variable. This problem only has 2 variables but, as the relaxed space is in dimension 11, SEGO with KPLS is then considered to compare with SEGO with Kriging such that we can see how reducing the number of spatial correlations from 11 to 5 in the model can impact the optimization process. For Bandit-BO, we represent 1 as a bandit with 10 arms from 0 to 9 and an initial DoE with a minimum of 20 points is required. As previously for NSGA2 and SEGO-KRG, two initial DoE sizes are considered with 5 points and 20 points. For a given DoE size, 20 different initial DoE are obtained to compute statistics on the convergence results.

In Fig. 2, the medians and the associated quartiles (25% and 75%) of the 20 runs are illustrated for each of the four algorithms. The initial DoE is shown before the black dotted line. For the DoE with 20 initial points in Fig. 2b, 50 iterations of the methods Bandit-BO, SEGO+KPLS(d=5) and, SEGO-KRG are performed for a total of 70 evaluations and for NSGA2 200 iterations are done. When the DoE size is reduced to 5 points in Fig. 2a, 200 iterations of SEGO-KRG are done and compared to 300 evaluations for NSGA2. These convergence plots show that the smaller the DoE, the faster the convergence. Also, KPLS slows the convergence at the start but this dimension reduction does not change the convergence overall and the incumbent is even better at the end with KPLS than without. On this test case, Bandit-BO method does not perform well. To compare the dispersion of the best results over the 20 runs, Fig. 2c shows the boxplots for SEGO-KRG, Bandit-BO, NSGA2 and, SEGO+KPLS(d=5) after 50 iterations associated to a initial DoE of 20 points. These plots show that SEGO+KPLS does not converge properly because 2 runs are outliers that did not have converged whereas with NSGA2 and SEGO-KRG, there is only one outlier. Nevertheless, NSGA2 is not as precise as SEGO-KRG on average.

Constrained optimization

For constrained optimization, Bandit-BO can not be considered as it can not deal with constraints, so SEGO is only compared to NSGA2. The first constrained test case is a modified Branin function [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF] with one constraint, denoted by "Branin 3", where the two first variables are categorical variables with 2 levels each such that there is 4 possible Branin function variations and the two last variables are the continuous ones. This problem has 4 variables in the initial space and 6 in the relaxed one, so SEGO with Kriging is applied without any dimension reduction technique. According to the previous experiments, a small initial DoE with 5 points is considered in order to obtain better results for a given number of evaluations. In Fig. 3a, the medians and the associated quartiles (25% and 75%) of the 20 runs are illustrated for NSGA2 and SEGO-KRG. The initial DoE is shown before the black dotted line, then 50 iterations of SEGO-KRG are performed and compared with 200 iterations of NSGA2. To compare the dispersion of the best results over the 20 runs, Fig. 3b shows the boxplots for SEGO-KRG and NSGA2 after 50 iterations. On this low-dimension constrained case, the mixed integer version of SEGO with Kriging is shown to perform well and be adapted to the constrained mixed optimization. The second constrained test case is an augmented dimension modified Branin function [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF] with one constraint, denoted by "Branin 4". The first two variables are categorical with 2 levels each and the 10 last variables are continuous. Therefore, we can see the test case as a set of 4 continuous augmented 10-dimensional Branin functions. This problem has 12 variables in the initial space and 14 in the relaxed space. SEGO with KPLS is considered with an automatic number of principal components to compare the dimension reduction with SEGO+KRG (equivalent to SEGO+KPLS(d=14)) and with SEGO+KPLS(d=2). As previously, an initial DoE of 5 points is chosen, shown before the dotted line in Fig. 4a. The medians and the associated quartiles (25% and 75%) of the 20 runs are illustrated for NSGA2 and SEGO. On this graph, 50 iterations of SEGO+KRG are compared with 200 iterations of NSGA2. To illustrate the dispersion of the best results over the 20 runs, Fig. 4b shows the boxplots for SEGO+KPLS(d=14), SEGO+KPLS(d=2), SEGO+KPLS(auto) and NSGA2 after 50 iterations. On this high-dimension constrained case, our method is still efficient. SEGO+KPLS(auto) and SEGO+KPLS(d=2) give similar results. To see how the different PLS models behave on the first run, we plotted on Fig. 5 the number of components through the 50 iterations. The adaptive method gives a mean number of 1.9 components for the objective model and of 1.25 components for the constraint model: the dimension reduction is quite efficient compared to KPLS(d=2) even in a few iterations more components are necessary (up to 5 for the objective function). For both models, on Fig. 6, the estimated PRESS error reduction (based on Eq. ( 6)) is reported through the iterations in which PLS components are added. We can see that this error reduction increases with the size of the DoE. Due to the fact that a small DoE leads to a lot of uncertainties, adding components gives marginal gains with sparse data. Contrarily, when the number of known points is larger, adding components to the model increases the quality of its prediction (up to 10%). These preliminary results on the 'Branin 4" test case motivate the interest for the adaptive method. Fig. 6 Error gained when choosing more than 1 component.

These four analytical cases have shown that SEGO performs better than both Bandit-BO for mixed integer Bayesian Optimization and NSGA2 for constrained optimization. It has been shown that the PLS technique allows the method to be scalable for high-dimension and even to give better results than Kriging by carrying favor to the most relevant space.

Overall benchmark results

For 10 different test cases (3 constrained and 7 unconstrained) we are considering 20 runs with different DoE sampled by Latin hypercube sampling for a total of 200 instances. The resulting percentages of instances that have converged after a given budget for every method are plotted on Fig. 7. For unconstrained test cases, in order to compare with Bandit-BO, the size of the initial DoE is given by min{5, 2 × } where is the number of categorical possibilities and for constrained test cases, we took 5 points for the initial DoE size. This allow us to compare Bandit-BO, SEGO, SEGO with Gower Distance, SEGO with KPLS and NSGA2. Some tests are in dimension 2, so in order to compare with the same number of hyperparameters for all tests, we had to choose between 1 or 2 principal components for KPLS. As the number of points increases, the projected points are really closed to each other, so, to insure a certain stability, we choose KPLS with 2 principal components and noise evaluation, denoted by KPLS(d=2) in the following.

These 10 test cases were extracted from several state-of-the-art papers [13-16, 19, 45, 46] from dimension 2 to 12, with integer, continuous and/or categorical variables. For constrained case, we keep only the inputs that give a constraint violation smaller than 10 -4 . A problem is considered solved if the error to the known solution is smaller than 2% on Fig. 7a and smaller than 0.5% on Fig. 7b. The mean error after 50 iterations can be found in Tab. 1. SEGO with PLS gives the smaller errors on constrained test cases but, for unconstrained ones, SEGO-KRG performs the best. As we can see on the data profiles, SEGO-KRG and SEGO+KPLS(d=2) are similar and outperform the three other methods. However, SEGO with PLS being less efficient on unconstrained test cases, SEGO-KRG gives better results over all tests. These preliminary results are promising and some more realistic applications are considered in the next section.

C. Aircraft design test cases

For the core MDO application, we apply FAST-OAD on two different aircraft design problems. The first one named "CERAS" is a classic well-known aircraft and the second "DRAGON" is a more innovative new one that is currently under development. This time, as the evaluations are expensive, we are doing only 10 runs instead of 20. For each run, we draw a random starting DoE of 5 points. As we have no equality constraint and a small budget, we will not force the constraints to be as large as possible and we will use only the mean prediction of the constraint models and not the criterion [START_REF] Priem | Upper trust bound feasibility criterion for mixed constrained Bayesian optimization with application to aircraft design[END_REF]. Also, to have realistic results, the constraints violation will be forced to be less than 10 -3 . On the following, let MAC denote the Mean Average Chord, VT be the Vertical Tail, HT be the Horizontal Tail and TOFL be the Take-off Field Length.

"CERAS" configuration

To validate the results of Section IV.B on a real test case, we are considering the data from the CEntral Reference Aircraft System ("CERAS") based on an Airbus A320 aircraft. The problem to solve is a constrained optimization problem with 6 continuous design variables, 2 categorical variables and 2 integer ones, for a total of 10 design variables. The presented version of SEGO (with or without the PLS technique) is used as an optimizer in a Multidisciplinary feasible (MDF) architecture where the MDA is computed with FAST-OAD. The optimization problem is described in Tab. 2 where the total number of variables is reported using the relaxation technique. We are testing three methods, NSGA2, SEGO with Kriging and SEGO with KPLS over 10 different DoE and 200 iterations. As for the analytical test cases, Fig. 8a shows that PLS method with 4 hyperparameters speeds up the optimization process at the start. However, the boxplots of the final results on Fig. 8b show that, at the end, the KPLS(d=4) and Kriging versions give similar results while KPLS uses 3 times less correlation length hyperparameters. The Kriging version starts slower but catches up with KPLS quickly and for every given budget Kriging gives better results than NSGA2. The best configuration is obtained after a SEGO+KPLS(d=4) optimization, the optimal result is given in Tab. 3. The optimal aircraft geometries obtained using each of the three methods and the baseline are plotted in Fig. 9. As SEGO with Kriging is almost the same as the best, SEGO+KRG and SEGO+KPLS are grouped on the geometries. With KPLS we obtain a constraint margin of 0.500 and an objective of 16722.73 kg. From an aircraft design point of view, the obtained optimized values can be surprising for some parameters, but are logical given the limits of current models in FAST-OAD:

• as the flight ceilings (buffeting, aerodynamics, ...) are not computed, the algorithm logically goes toward high altitudes, • the mass of a strongly swept wing is not properly computed, and does not bring enough penalty on such configuration, • similarly, the models do not compute the additional mass that should be required for a T-tail with such a large aspect ratio on vertical tail. However, the capabilities of the proposed algorithm are very promising for a conceptual design stage where a lot of architecture choices are still undetermined, leading to a large combinatorial problem. In that perspective, the ability of the method to capture the right trends regarding number of engines (integer) and their position (categorical) is a good perspective, that needs to be confirmed with an even more complex design case using many categorical and integer variables. The "DRAGON" hybrid aircraft case described below provides this feature. Fig. 9 "CERAS" best configuration geometries. Comparisons between the initial configuration, the NSGA2, the SEGO+KRG and the SEGO+KPLS(d=4) results.

"DRAGON" configuration

The "DRAGON" aircraft concept in Fig. 10 has been introduced by ONERA in 2019 [START_REF] Schmollgruber | Multidisciplinary Exploration of DRAGON: an ONERA Hybrid Electric Distributed Propulsion Concept[END_REF] within the scope of the European CleanSky 2 program † which sets the objective of 30% reduction of CO2 emissions by 2035 with respect to 2014 state of the art. A first publication in SciTech 2019 [START_REF] Schmollgruber | Multidisciplinary Exploration of DRAGON: an ONERA Hybrid Electric Distributed Propulsion Concept[END_REF] was followed by an up-to-date report in SciTech 2020 [START_REF] Schmollgruber | Multidisciplinary Design and performance of the ONERA Hybrid Electric Distributed Propulsion concept (DRAGON)[END_REF]. ONERA responded to this objective by proposing a distributed electric propulsion aircraft that improves the aircraft fuel consumption essentially by increasing the propulsive efficiency. Such efficiency increase is obtained through improvement of the bypass ratio by distributing a large number of small electric fans on the pressure side on the wing rather than having large diameter turbofans. This design choice avoids the problems associated with large under-wing turbofans and at the same time allows the aircraft to travel at transonic speed. Thus the design mission set for "DRAGON" is 150 passengers over 2750 Nautical Miles at Mach 0.78. The employment of a distributed propulsion comes at a certain cost; a turbo-electric propulsive chain is necessary to power the electric fans which brings additional complexity and weight. Typically, turboshafts coupled to electric generators are generating the electrical power on board of the aircraft. The power is then carried to the electric fans through an electric architecture sized to ensure robustness to single component failure. This safety feature is obtained with redundant components as depicted in Fig. 11. The baseline configuration is two turboshafts, four generators, four propulsion buses with cross-feed and forty fans. This configuration was selected for the initial study as it satisfies the safety criterion. However it was not designed to optimize aircraft weight. The turboelectric propulsive chain being an important weight penalty, it is of particular interest to optimize the chain and particularly the number and type of each component, characterized by some discrete or particular values. In Aerobest [START_REF] Saves | Constrained Bayesian optimization over mixed categorical variables, with application to aircraft design[END_REF], we solved a constrained optimization problem with 8 continuous design variables and 4 integer ones, for a total of 12 design variables. Like for "CERAS", we found that the best configuration with SEGO and Kriging with and without PLS regression were almost the same while PLS reduced the number of variables from 12 to 4. These first results were promising but no categorical variable was considered. Here the motivation is to add some of these variables in order to increase the number of considered possibilities. Moreover, the MDA fixed-point algorithm is now replaced to be more flexible and general as the "DRAGON" test case is treated with Overall Aircraft Design method in FAST-OAD [START_REF] David | From FAST to FAST-OAD: An open source framework for rapid Overall Aircraft Design[END_REF].

In [START_REF] Saves | Constrained Bayesian optimization over mixed categorical variables, with application to aircraft design[END_REF], the observation of the optimization results from an aircraft point of view led to surprising conclusions, ultimately resulting in questioning the validity of the models utilized for the design. Among the problems raised in this previous work, the sizing rules of the propulsive chain, based on the default architecture, proved to be insufficiently general to explore more diversified and valid architectures. The results indicated that an architecture comprising 40 motors connected to 6 cores of equal power would be optimal, although it is violating the design rule stating that each motor should be connected to a single core.

In this paper, we propose an up-to-date optimization problem for "DRAGON" with the most recent enhancements. The consideration of the variables related to architecture was revised to take full advantage of mixed variables optimization. Three configurations with different number of electric components were considered, each with their own sizing rules. The default configuration is conserved for the analysis and two new configurations, one with low distribution, the other with high distribution, were created and analysed to establish the sizing rules. The choice of architecture is hence represented by a categorical variable. The number of motors was to remain an optimization variable as it is an important driver of the propulsive and aerodynamic efficiency of the aircraft but it is constrained by the type of architecture. In this analysis, the number of motors could only be multiple of 4, 8 or 12 depending of the selected architecture. The solution consisted in expending the levels of the categorical variable representing the architecture and assign a specific and valid number of motors to each level.

Additionally, it was noticed that the layout of "DRAGON", with the turbogenerators at the rear of the fuselage was disadvantageous for two reasons:

• The added weight at the rear forces the wing to be more aft, reducing the level arm with the tail surfaces. To maintain the static stability, the tail surfaces have to be increased, adding weight and friction drag penalty. • The electric cables running from the rear generators to the center part of the wing are heavy, accounting for ∼10% of the total propulsion system weight. For these two reasons, it would be advantageous to locate the turbogenerator below the wing. However, this would restrict the space available for the electric fans, which would in turn limit the maximum propulsive efficiency. To make the trade off between lighter propulsion and better propulsive efficiency the layout was included into the optimization problem as a categorical variable.

Five constraints are added among which, the take off field length, climb duration and top of climb slope angle are strong drivers for the sizing of the hybrid electric propulsion. A portion of the wing trailing edge at the wingtip, should be left free for the ailerons hence limiting the space available for the electric fans. Finally, the wingspan limit is imposed by airport regulation.

To know how optimizing the fuel mass will impact the aircraft design, we are considering the new optimization problem described in Tab. 4. We can now solve a constrained optimization problem with 10 continuous design variables and 2 categorical variables with 17 and 2 levels respectively, for a total of 12 design variables. For the optimization, this new problem is a hard test case involving 29 relaxed variables and 5 constraints. The definition of the architecture variable is given in Tab. 5 and the definition of the turboshaft layout is given in Tab. 6. To validate our method, we are comparing two variants of SEGO with automatic KPLS or with its more expensive version using Kriging, and NSGA2 in Fig. 12 where 10 runs are performed in order to plot the mean and the associated quartiles. The best method after 10 iterations is the proposed one involving automatic PLS regression as shown in the boxplots in Fig. 12b. In Fig. 12a, after 100 iterations, we still find that SEGO is better than NSGA2 and that the PLS technique helps for the convergence. As for "CERAS", we find that the best configurations found with SEGO and Kriging with and without PLS regression are almost the same. The description of this best configuration is given in Tab. 7 and the best geometries for NSGA2, SEGO+KRG and SEGO+KPLS(auto) are given in Fig. 13. SEGO+KRG and SEGO+KPLS(auto) are grouped as they give similar results. The optimal configuration was found with Kriging for 10816 kg of fuel and KPLS found a similar point at 10819 kg. The best configuration found is the configuration 10, with the smaller number of 8 motors but with 4 cores and electric generators.

On the point of view of the aircraft design, the number of motors was restricted to the minimum available and the two architectures allowing 8 motors (1 and 10) are performing almost identically, showing a low influence of this variable over the number of motors. In this case, it is the climb constraints that are sizing for the whole propulsion chain. Therefore the architecture sizing laws are not at play and no difference can be seen.

Despite an important space still available at the wing trailing edge, a small number of motors is looked for by the optimizer to lower the wetted surface area of the fans. This is probably oversimplified as the integration of the fans at the trailing edge is coarsely modeled. The fact that this direction is selected during the optimization tells us that this model should be updated to further improve the design.

The most favorable layout is found to be with the turbo-generators at the rear. The level arm between the wing and the horizontal tail being actually larger due to the maximum sweep angle employed for the horizontal tail. As for the "CERAS" test case, the combination of high sweep and high aspect ratio is too badly taken into account on a structural point of view and leads to unrealistic weight for the horizontal stabiliser, which probably favour too much this layout. Nevertheless, the trade-off found by optimization is correct given the models used in FAST-OAD.

V. Conclusion and perspectives

To conclude, we have observed on both analytical and industrial cases than SEGO with KPLS is well-suited and efficient for a mixed integer high-dimensional constrained efficient global optimization problem. We have seen that continuous relaxation allows straightforward use of continuous GP but can be impractical as it increases the computational effort required to build the surrogate model. By using the PLS regression, it is possible to reduce the computational cost and makes the continuous relaxation affordable in practical contexts. Using the Wold's R criterion, it is possible to automatize the PLS regression for high-dimensional problems. This method has been applied to aircraft design optimization and contributed to the development of a MDA tool designing future aircraft configurations.

4 )

 4 Maximize the acquisition function 2 over Ω : = arg max ∈Ω 2 ( ). 5) Add , ( ), ( ), ℎ( ) to the DoE D +1 . Increment .

: 2 ) 4 ) 5 )

 245 Solution of the problem (1) over the mixed categorical design space Ω × × F . Inputs: Initial DoE D 0 and set = 0. The search space Ω × × F . while the stopping criterion is not satisfied do 1) Relax continuously integer and categorical input variables to a real bounded space Ω of dimension = + + =1 . Namely, we continuously relax the mixed categorical DoE D = { } ∈ (Ω × × F ) to a continuous DoE D = { } ∈ (Ω ) using the relaxation procedure . Build the GP models for the objective function and the constraints and ℎ related to the DoE D with PLS to reduce the number of the hyperparameters from to , and ℎ as detailed in Alg. 3. 3) Build an estimation of the feasible domain Ω ⊂ Ω with the criterion and construct the acquisition function 2 . Maximize the acquisition function 2 over Ω : = arg max ∈Ω 2 ( ). Project the obtained continuous solution over Ω × × F : = ( ). 6) Add , ( ), ( ), ℎ( ) to the DoE D +1 . Increment .

Algorithm 3 :

 3 Adaptive dimension reduction for KPLS models. Result: Optimal model of reduced order Inputs: A DoE of size , i.e, D = {( 1 , 1 ), . . . , ( , )} associated with a given function Ψ (Ψ( ) = , ∀ = 1, . . . , ). The maximal and minimal number of components max and min , respectively. A threshold ∈ [0, 1]. A chosen number of folds . ∀ ∈ {1, . . . , } denote the corresponding known points: {( 1 , 1 ), . . . , ( , )}.

  with a DoE of 32 points.

Fig. 1 "

 1 Fig. 1 "Branin 5" obtained optimization results. The Boxplots are generated, after 50 iterations, using the 20 best points.

Fig. 2 "

 2 Fig. 2 "Set 1" obtained optimization results.

Fig. 3 "

 3 Fig. 3 "Branin 3" optimization results.

Fig. 4 "Fig. 5 "

 45 Fig. 4 "Branin 4" optimization results.

Fig. 7

 7 Fig. 7 Data profiles generated using 10 analytical test cases.

Fig. 8 "

 8 Fig. 8 "CERAS" optimization results using a DoE of 5 points. The Boxplots are generated, after 200 iterations, using the 10 best points.

Fig. 10 "

 10 Fig. 10 "DRAGON" aircraft mock-up.

Fig. 11 turboelectric

 11 Fig. 11 turboelectric propulsive architecture.

Fig. 12 "

 12 Fig. 12 "DRAGON" optimization results using a DoE of 5 points over 10 runs. The Boxplots are generated, after 100 iterations, using the 10 best points.

Fig. 13 "

 13 Fig. 13 "DRAGON" best configuration geometry. Comparisons between the initial configuration, the NSGA2, the SEGO+KRG and the SEGO+KPLS(auto) results.

Table 1 Mean errors of each method after 50 iterations.

 1

	ERRORS	NSGA2 Bandit-BO SEGO+KRG SEGO+KPLS(d=2) SEGO+GD
	Unconstrained test cases 14.29%	5.84%	2.27%	5.74%	8.41%
	Constrained test cases	58.27%	-	27.61%	25.41%	47.83%

Table 2 Definition of the "CERAS" optimization problem.

 2

		Function/variable	Nature	Quantity	Range
	Minimize	Fuel mass	cont	1
	with respect to	x position of MAC	cont	1	[16., 18.] ( )
		Wing aspect ratio	cont	1	[5., 11.]
		Vertical tail aspect ratio	cont	1	[1.5, 6.]
		Horizontal tail aspect ratio cont	1	[1.5, 6.]
		Wing taper aspect ratio	cont	1	[0., 1.]
		Angle for swept wing	cont	1	[20., 30.] ( • )
		Total continuous variables	6
		Cruise altitude	discrete	1 {30k,32k,34k,36k} ( )
		Number of engines	discrete	1	{2,3,4}
		Total discrete variables	2
		Tail geometry	cat	2 levels	{T-tail, no T-tail}
		Engine position	cat	2 levels	{front or rear engines}
		Total categorical variables	2
		Total relaxed variables	12
	subject to	0.05 < Static margin < 0.1 cont	2
		Total constraints		2

Table 3 "CERAS" Optimal aircraft configuration.

 3

	Name	Nature	Value
	Fuel mass	cont	16722.55 kg
	static margin	cont	0.0495
	x position of MAC	cont	16.2825 m
	Wing aspect ratio	cont	11
	VT aspect ratio	cont	6
	HT aspect ratio	cont	6
	Wing taper ratio	cont	0.5099
	Wing sweep angle	cont	30.0 •
	Cruise altitude	discrete	36,000 ft
	Tail geometry	cat	T-tail
	Engine position	cat	Front engine
	Number of engines discrete	2

Table 4 Definition of the "DRAGON" optimization problem.

 4

		Function/variable		Nature Quantity	Range
	Minimize	Fuel mass		cont	1
	with respect to Fan operating pressure ratio		cont	1	[1.05, 1.3]
		Wing aspect ratio		cont	1	[8, 12]
		Angle for swept wing		cont	1	[15, 40] ( • )
		Wing taper ratio		cont	1	[0.2, 0.5]
		HT aspect ratio		cont	1	[3, 6]
		Angle for swept HT		cont	1	[20, 40] ( • )
		HT taper ratio		cont	1	[0.3, 0.5]
		TOFL for sizing		cont	1	[1800., 2500.] ( )
		Top of climb vertical speed for sizing	cont	1 [300., 800.]( /	)
		Start of climb slope angle		cont	1	[0.075., 0.15.](	)
		Total continuous variables			10
		Architecture		cat	17 levels {1,2,3, . . . ,15,16,17}
		Turboshaft layout		cat	2 levels	{1,2}
		Total categorical variables			2
		Total relaxed variables			29
	subject to	Wing span < 36 ( )		cont	1
		TOFL < 2200 ( )		cont	1
		Wing trailing edge occupied by fans < 14.4 ( ) cont	1
		Climb duration < 1740 ( )		cont	1
		Top of climb slope > 0.0108 (	)	cont	1
		Total constraints			5

Table 5 Definition of the architecture variable and its 17 associated levels.

 5

	Architecture number number of motors number of cores number of generators
	1	8	2	2
	2	12	2	2
	3	16	2	2
	4	20	2	2
	5	24	2	2
	6	28	2	2
	7	32	2	2
	8	36	2	2
	9	40	2	2
	10	8	4	4
	11	16	4	4
	12	24	4	4
	13	32	4	4
	14	40	4	4
	15	12	6	6
	16	24	6	6
	17	36	6	6

Table 6 Definition of the turboshaft layout variable and its 2 associated levels

 6

	layout	position	y ratio	tail	VT aspect ratio VT taper ratio
	1	under wing	0.25	without T-tail	1.8	0.3
	2	behind	0.34	with T-tail	1.2	0.85

Table 7 "DRAGON" Optimal aircraft configuration.

 7

	Name	Nature	Value
	Fuel mass	cont	10816 kg
	Wing span	cont	36 m
	TOFL	cont	1722.7 m
	Wing trailing edge occupied by fan	cont	10.65 m
	Climb duration	cont	1735.3 s
	Top of climb slope	cont	0.0108 rad
	Fan operating pressure ratio	cont	1.09
	Wing aspect ratio	cont	10.9
	Angle for swept wing	cont	32.2 •
	Wing taper ratio	cont	0.235
	HT aspect ratio	cont	6
	Angle for swept HT	cont	40 •
	HT taper ratio	cont	0.3
	TOFL for sizing	cont	1803 m
	Top of climb vertical speed for sizing	cont	494 ft/min
	Start of climb slope angle	cont	0.104 rad
	Architecture	cat	10
	Turboshaft layout	cat	2

† https://www.cleansky.eu/technology-evaluator

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 consortium for their contribution and feedback. We also thanks the people who worked on DRAGON and made this work possible: Carsten Döll, David Donjat, Michael Ridel, Jean Hermetz, Bernard Paluch and Olivier Altinault.


Luc Lesoil

Hugo Martin

Mathieu Acher

Arnaud Blouin

Jean-Marc Jézéquel

Transferring Performance between Distinct Configurable Systems : A Case Study

Many research studies predict the performance of configurable software using machine learning techniques, thus requiring large amounts of data. Transfer learning aims to reduce the amount of data needed to train these models and has been successfully applied on different executing environments (hardware) or software versions. In this paper we investigate for the first time the idea of applying transfer learning between distinct configurable systems. We design a study involving two video encoders (namely x264 and x265) coming from different code bases. Our results are encouraging since transfer learning outperforms traditional learning for two performance properties (out of three). We discuss the open challenges to overcome for a more general application.

INTRODUCTION

Software offers more and more options that users can (de)select to customize the system for their specific needs. With the exploding number of options e.g., +2041 options in three years for the Linux kernel [START_REF] Martin | Transfer Learning Across Variants and Versions : The Case of Linux Kernel Size[END_REF], it becomes complex to accurately estimate the individual impact of options, difficult to predict software performance and unthinkable to measure exhaustively the configuration space of real-world systems. To overcome this problem, related work has proposed to train machine learning models a.k.a. performance models [START_REF] Guo | Variability-Aware Performance Prediction: A Statistical Learning Approach[END_REF] taking software configurations as input, training the model to link a configuration to a performance value and finally capable of predicting the performance of a new configuration. But this method has a cost [START_REF] Nair | Transfer Learning with Bellwethers to find Good Configurations[END_REF]: since we feed the model with configurations and performance, we have to measure numerous configurations. Transfer learning techniques [START_REF] Martin | Transfer Learning Across Variants and Versions : The Case of Linux Kernel Size[END_REF][START_REF] Moradi | Performance Prediction in Dynamic Clouds using Transfer Learning[END_REF] aim at reducing this measurement cost: if a performance model has already been trained on one executing environment, we can reuse this model to reduce the number of measurements needed for a second executing environment [START_REF] Jamshidi | Learning to Sample: Exploiting Similarities across Environments to Learn Performance Models for Configurable Systems[END_REF], assuming that both environments are similar (e.g., same operating system). We propose to apply transfer learning to distinct software systems performing the same task, such as compilers (e.g., gcc and llvm), container managers (e.g., podman and docker), etc. Intuitively, the model trained on one software could be used -at least partially-to train the other performance model, as depicted in Figure 1 for x264 and x265. This paper presents the first minimal example showing that under certain conditions, it is possible to transfer performance models across software systems. We also discuss the limitations of our work and highlight the open challenges to face when scaling to other software systems. Data1 and code2 are publicly available.

TRANSFER LEARNING ACROSS SOFTWARE SYSTEMS: A PROOF OF CONCEPT 2.1 Research Questions

To get started with this topic, we first study the differences of performance between different software systems executing the same task. Are the distribution of performance consistent across similar software systems? If they are, we can probably use similarities between their performance distributions when training the model. We then address the main issue of paper : Can we transfer performance from one software system to another? How useful is it compared to non-transfer technique? How much do we gain, both in terms of measurements and accuracy? More specifically, we want to ensure that transfer learning does indeed outperform simple learning and avoid any instance of negative transfer.

To answer these questions, we design the following study. 4 . We exploit these common features and make their values vary in order to generate a set of 3125 configurations working for both systems. We check the uniformity of the resulting distribution of options' values with a Kolmogorov-Smirnov test [START_REF] Frank | The Kolmogorov-Smirnov test for goodness of fit[END_REF] 5 .

Study Design

Input Data. We select eight input videos [START_REF] Lesoil | The Interaction between Inputs and Configurations fed to Software Systems: an Empirical Study[END_REF] extracted from the Youtube UGC Dataset [START_REF] Wang | YouTube UGC Dataset for Video Compression Research[END_REF], well-known in the community of video compression. For this selection, we vary the content (LiveMusic, Sports) and the resolution (360P, 480P) of videos.

Performance Properties. We then use x264 and x265 to transcode these eight videos from the mkv to the mp4 format. During each execution, we measure the percentage of cpu usage, etime the elapsed time in seconds and the file size of the resulting video in bytes.

Executing Environment. We measure all performance sequentially on a dedicated (and warmed-up) server -model Intel(R) Xeon(R) CPU D-1520 @ 2.20GHz, running Ubuntu 20.04 LTS.

Compare Performance.

We analyze the differences between the distribution of performance properties of x264 and x265. As a measure of (dis)similarity, we compute their Spearman rankbased correlation [START_REF] George | Rank correlation methods[END_REF]. It is suited for our case since all performance properties are quantitative variables relative to the same configurations. If both encoders obtain the same rankings in terms of performance, the correlation is close to 1, and there is a good chance of getting good results with transfer learning. If they react differently to the same configurations, the correlation is close to 0 and the transfer might be challenging to achieve.

Transfer

Performance. Third, and as displayed in Figure 1, we try to transfer the performance from x264 (i.e., source software) to x265 (i.e., target software). We use Model Shift, a simple and stateof-the-art transfer learning approach defined by Valov et al. [START_REF] Valov | Transferring Performance Prediction Models Across Different Hardware Platforms[END_REF]. The protocol should be read following Figure 2 : 1. First, it trains a shifting function, mapping the performance distribution of the source on the target software's performance distribution, 2. Then, it trains a performance model on the source software, 3. Finally, it predicts the performance distribution of the 4 The list of selected configuration options can be consulted here : https://anonymous. 4open.science/r/TL_cross_soft-855B/replication/x26x/README.md 5 Results can be consulted at: https://anonymous.4open.science/r/TL_cross_soft-855B/ replication/x26x/x264_x265_options.ipynb source software and applies the shifting function to the predictions, in order to estimate the performance of the target software system. We compare Model Shift to a simple baseline acting as a control approach, training a performance model directly on the target software, without using any measurement of the source. We call this baseline No Transfer. For both, we used a Random Forest algorithm [START_REF] Thais | How many trees in a random forest?[END_REF] to predict software performance, without tuning its hyperparameters. We separate the dataset of the target into training and test, varying the size of the training set. In the evaluation, we compare and display the Mean Absolute Percentage Error (MAPE) [START_REF] Molnar | Interpretable Machine Learning[END_REF] between the predicted values (i.e., predicted by the approaches) and the real values (i.e., measured on the test set of the target software). To reduce the variance induced by machine learning randomness, we repeated the process five times and display the average MAPE for the test sets. We rely on the python library scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. 62.3 Evaluation 2.3.1 Compare Performance. Among the three performance properties, we can distinguish three cases : between x264 and x265, the cpu consumption has an average a correlation close to 0 (lower than 0.1), the elapsed time etime is overall positively correlated (about 0.75), and there is almost no differences between the sizes of the resulting mp4 videos (correlations close to 1). Transferring sizes is likely to be easy and transferring cpu consumption will probably be too difficult. The elapsed time etime is in between. equivalent to the baseline whatever the number of configurations and the video. This can be explained by the high correlations observed in Table 1. At first, the transfer is really outperforming the baseline e.g., in Figure 3h, for 5 configurations, the transfer has an error of about 106 ± 64 % while the baseline amounts to 1071 ± 556 %.

Conclusion.

The effectiveness of the transfer varies with the performance property we consider. Overall, it is possible to outperform the No Transfer baseline for 2 performance properties (out of 3), especially when the budget (i.e., the number of configurations in the training set) is low. As noticed by [START_REF] Jamshidi | Transfer Learning for Performance Modeling of Configurable Systems: An Exploratory Analysis[END_REF], the greater the correlation between performance distributions, the more accurate the transfer. As a concrete advice for developers, this correlation could be a cheap indicator to estimate a priori whether transfer techniques are adapted between two systems.

THREATS TO VALIDITY

Due to the cost of measurements (58 days of system time), we did not measure performance more than once. Therefore, the performance distribution of x264 and x265 could change with new measurements. To address this threat, we check the results for eight different input videos; given their consistency, we are confident that similar conclusions could be drawn by reproducing the experiment. Another threat to validity is related to the use of machine learning algorithms, which leads to non-deterministic results. To mitigate this threat, we display the average result of five runs of the model. In Section 2.2, we only select the configuration options common to x264 and x265, ignoring a large majority of features. This represents a potential threat to validity when generalising the transfer to the whole configuration space.

DISCUSSION

Although our experiment is encouraging to further explore this research direction, it does not cover all possible cases of transfer learning between software systems. In this section, we identify three lines of research and for each [START_REF] Mathieu Acher | Composing Feature Models[END_REF] we point out the limitations of our study; [START_REF] Mathieu Acher | Slicing Feature Models[END_REF] we discuss the open challenges to overcome when transferring performance between distinct software systems; (3) we describe the potential approaches and solutions.

Find Transferable Software Systems

4.1.1 Limitation. A threat to our study is the choice of x264 and x265. There might be other pairs of systems for which nonfunctional properties are dissimilar e.g., with a correlation of zero, the transfer learning might perform poorly. For x264 and x265 our experiments show that two non-functional properties out of three are positively and strongly correlated.

Open challenge.

We cannot guarantee the effectiveness of transfer learning for every pair of software systems. What is difficult is to know whether transfer will work for a given pair of software systems.

Possible Solution.

For now, the only reasonable assumption we can make is to choose software systems within the same domain e.g., compilers like gcc, llvm or clang, container managers like podman or docker, learning libraries like theano, pytorch or tensorflow, text editors like emacs, gedit or vim, video encoders like vp9, x264, or x265 Measures on how performance distributions differ (e.g., with correlations) across systems can provide a first indicator on whether transferring is worth.

When and How to Transfer?

4.2.1 Limitation. The results of Section 2 show some limits of transfer learning: after a given training size, the interest in adding noisy measures -such as source data -decreases. This learning size is a switching point, as shown by Figure 3a: before this point, it is preferable to use transfer learning and after this point, we should switch to simple learning. This point seems difficult to estimate a priori i.e., without any measurement. This may depend on the complexity of the configuration space: the more complex the configuration space, the more configurations are needed to make an accurate prediction. And the more configurations needed to be accurate, the larger the switching point.

Open Challenge.

It is challenging to know when to apply transfer learning and when to switch to non-transfer learning. These results could also be better (or worse) with other transfer learning techniques (see Section 5). The challenge here is to determine the best learning approach to use for a pair of software systems.

Possible Solution.

Empirical studies comparing different approaches on representative pairs of systems would be of great help in addressing this challenge. Specialized learning algorithms, capable of handling distribution shifts across software systems, are also expected.

What to Transfer?

4.3.1 Limitation. Our choice of configuration options is another limitation of our study. We have deliberately considered a favorable case: since x264 and x265 have been developed by the same team, the two systems share a common set of configuration options with the same range of values.

Open Challenge.

In general, the configuration spaces will often be very different between the source and the target systems, making the sampling of configurations difficult in practice. The challenge is to map the configuration space of the source to the configuration space of the target. We illustrate these differences between configuration spaces with vp9, x264, and x2657 :

• The same feature is implemented in the source and the target with different names e.g., -level for x264 and -level-idc for x265; • The same feature is implemented in the source and the target, but a value is only implemented in one software system e.g., unlike for x265, the motion-estimator feature -me of x264 does not implement the 'star' pattern search; • The feature of one system encapsulates one feature (or more) of the other e.g., activating -fullrange in x264 is equivalent to choose -range full for x265; • The feature is only implemented in one software system e.g., the feature -rc-grain of x265 does not exist for x264; • The feature does not have the same default value for the source and the target e.g., -qpmax is set to 51 by default for x264 and set to 69 by default for x265;

• Both software systems do not have the same requirements or feature interactions e.g., for x265, passing a input video in the yuv format does not work unless you specify -input-res, while it does for x264; • The same feature is implemented in the source and the target, but the scale of the values differ between the source and the target.

For instance, the constant rate factor -crf goes from 0 to 63 for vp9 but from to 0 to 51 for x264 and x265, which is problematic when comparing a value (e.g., -crf = 35) that will not have the same meaning for all systems. We could even imagine a situation where the values of a configuration option are increasing for one system and decreasing for the other.

In addition, and as acknowledged by [START_REF] Krüger | Towards a Better Understanding of Software Features and Their Characteristics: A Case Study of Marlin[END_REF], the study of configuration options also requires domain knowledge, which is also true when mapping the configuration spaces of different software systems.

Possible Solutions.

Mapping configuration spaces could be envisioned using recent advances in the variability community. We propose a very simple protocol: (i) Align features [START_REF] Linsbauer | Variability extraction and modeling for product variants[END_REF][START_REF] Schulze | Aligning Coevolving Artifacts Between Software Product Lines and Products[END_REF] between systems; (ii) Meticulously model both configuration spaces e.g., with feature models [START_REF] Benavides | Automated Reasoning on Feature Models[END_REF][START_REF] Kyo | Feature-Oriented Domain Analysis (FODA) feasibility study[END_REF]; (iii) Analyse [START_REF] Di | Feature Diagrams as Package Dependencies[END_REF] and instrument [START_REF] Mathieu Acher | Composing Feature Models[END_REF][START_REF] Mathieu Acher | Slicing Feature Models[END_REF] them in order to create a resulting feature model generating configurations accepted by both systems. If the objective is purely to predict performance, heterogeneous transfer learning may be a possible black-box alternative (see Section 5) that infers the existing relationships between distinct configuration spaces.

RELATED WORK

Performance Comparison. In many areas, related work compares the performance of different tools performing the same task [START_REF] Boyko | Performance Evaluation and Comparison of Software for Face Recognition, Based on Dlib and Opencv Library[END_REF][START_REF] Rasmus Emilsson | Container performance benchmark between Docker, LXD, Podman &amp[END_REF][START_REF] Li | Performance comparison and evaluation of software tools for microRNA deep-sequencing data analysis[END_REF][START_REF] Park | Performance comparison of GCC and LLVM on the EISC processor[END_REF][START_REF] Pereira | Software Variability Management: An Exploratory Study with Two Feature Modeling Tools[END_REF][START_REF] Pohl | A Performance Comparison of Contemporary Algorithmic Approaches for Automated Analysis Operations on Feature Models[END_REF]. These empirical evidences could be relevant; if two software systems have similar performance distributions, they are good candidates for the transfer.

Transfer Learning. We list here few other transfer learning techniques that could be applied across software systems. Jamshidi et al. define Learning to Sample (L2S) [START_REF] Jamshidi | Learning to Sample: Exploiting Similarities across Environments to Learn Performance Models for Configurable Systems[END_REF]: it combines an exploitation of the source and an exploration on the target to sample a list of configurations, improving the overall accuracy of the learning algorithm. As many other transfer learning works [START_REF] Ballesteros | Transfer Learning for Multiobjective Optimization Algorithms Supporting Dynamic Software Product Lines[END_REF][START_REF] Moradi | Performance Prediction in Dynamic Clouds using Transfer Learning[END_REF], it is applied to transfer performance of executing environments. Martin et al. develop TEAMs [START_REF] Martin | Transfer Learning Across Variants and Versions : The Case of Linux Kernel Size[END_REF], a transfer learning approach predicting the performance distribution of the Linux kernel using the measurements of its previous releases. Between two releases, related the same system but distant in time, one could consider that it is a simple case of transfer across systems. Krishna et al. implement BEETLE [START_REF] Krishna | Whence to Learn? Transferring Knowledge in Configurable Systems using BEETLE[END_REF], that we could use to find one bellwhether software i.e., a source software that lead to better transfer results whatever the target. João et al. propose Weighted Multisource Tradaboost [START_REF] Antunes | Weighted Multisource Tradaboost[END_REF]. Applied to our case, it would exploit multiple source software systems and give them different weight values based on their similarity with the target system. This approach can be useful to apprehend the wide diversity of existing systems.

Heterogeneous Transfer Learning. Heterogeneous Transfer Learning [START_REF] Day | A survey on heterogeneous transfer learning[END_REF][START_REF] He | Heterogeneous transfer learning with weighted instance-correspondence data[END_REF][START_REF] Wu | Online Transfer Learning with Multiple Homogeneous or Heterogeneous Sources[END_REF] (HTL) is an extension of homogeneous transfer learning handling the differences between the source and the target feature spaces (a.k.a., configuration spaces when studying software variability). It creates a representation of the feature space, in between the source and the target, and finally transforms both feature spaces so they fit in this representation. Applied to the transfer across software systems, it would handle the changes of features between the configuration space of the source software and the configuration space of the target software. Moving to HTL would generalise our work to a more realistic case, where we can decide to select any feature of the source and the target systems.

CONCLUSION

This short paper shows it is challenging yet possible to transfer performance between distinct configurable software systems, using two video encoders (namely x264 and x265). We also discuss the limitations of our work and highlight the open challenges to overcome when generalising the applicability of transfer learning across distinct systems.

Figure 1 :

 1 Figure 1: Can x264 be used to predict x265 performance?

Figure 2 :

 2 Figure 2: Model Shift, an approach of transfer learning

Figure 3 :

 3 Figure 3: Mean Average Percentage Error (y-axis, lower is better) when transferring x264 to x265 performance depending on the training size (x-axis, log scale), for eight input videos and three performance properties • The transfer of cpu consumptions is almost always negative. For example, after 15 configurations in the training set for the Cov-erSong video, in Figure 3d, the baseline is always more accurate than transfer learning. For the Lecture video, in Figure 3e, it is even always a negative transfer. • For etime, the transfer is cost-effective at first, until a given training size e.g., 65 for Figure 3c. After this threshold, the error of the baseline keeps dropping while the errors stabilise for the transfer. For the Animation video, in Figure 3a, and for a budget of 500 configurations, No Transfer decreases to a MAPE of 35±1.3 % and Model Shift stays at 220 ± 21 %.• Finally, for the encoded sizes of videos, Model Shift is at least equivalent to the baseline whatever the number of configurations and the video. This can be explained by the high correlations observed in Table1. At first, the transfer is really outperforming the baseline e.g., in Figure3h, for 5 configurations, the transfer has an error of about 106 ± 64 % while the baseline amounts to 1071 ± 556 %.

  We select x264 and x265 3 , two video encoders. x264 and x265 realize the same task, but with different compression standards (resp. H.264 and H.265). It has a profound impact on the visual quality algorithms internally implemented. Though x265 has the ambition to borrow heavily from x264's features, x265 is not directly based on x264 source code. x265 is fully developed in C++ (and assembly) whereas x264 is written in C. x265 also implements novel algorithms such as CU-Tree the successor to x264's macroblock-tree. Overall, x264 and x265 are two distinct software projects (i.e., x264 is not a version or a fork of x265).

	2.2.1 Measure Performance. We gather measurements about the performance of two software systems. Software Systems. Importantly, both are developed by VideoLAN, which makes it easier to find similarities between them (same options, same conventions, etc.). From this respect, x264 and x265 can be seen as a favorable yet challenging case of transfer learning across systems. Configuration Options. We search for common configuration options in their documentation. For instance, and according to their documentation, both x264 and x265 implement the features -ref and -preset. -ref could be set to 1, 8 or 16 while -preset can be fast or slow. Possible resulting configurations like (slow, 1) or (fast, 16) are accepted and valid for x264 and x265. In the end, we keep 35 configuration options common to x264 and x265 (out of resp. 118 and 168)

Table 1 :

 1 Correlations between the performance distributions of x264 and x265 for eight input videos

	Video	cpu etime size
	Animation	0.05	0.74 0.98
	CoverSong	0.0	0.73 0.98
	Gaming	-0.02	0.81 0.99
	Lecture	0.07	0.75 0.98
	LiveMusic	0.01	0.77 0.99
	LyricVideo	0.14	0.73 0.96
	MusicVideo 0.07	0.75 0.99
	Sports	-0.0	0.78 0.98

2.3.2 Transfer Performance.

In Figure

3

, we depict the results varying for different inputs and performance properties:

Dataset available at : https://zenodo.org/record/5662589

Code available at : https://anonymous.4open.science/r/TL_cross_soft-855B/

See x264 and x265 webpages : https://www.videolan.org/developers/x264.html and https://www.videolan.org/developers/x265.html

A description of our python environment can be consulted at: https://anonymous. 4open.science/r/TL_cross_soft-855B/replication/requirements.txt

One can verify our illustrations with the lists of features, available at: https://cinelerra-gg.org/download/CinelerraGG_Manual/VP9_parameters.html for vp9, http://www.chaneru.com/Roku/HLS/X264_Settings.htm for x264, https://x265.readthedocs.io/en/2.5/cli.html for x265.

This research was funded by the ANR-17-CE25-0010-01 VaryVary project.


Kévin Hérissé

Benoit Larras

Bruno Stefanelli

Antoine Frappé

Andreas Kaiser

Mixed-Signal In-Memory Multi-bit Matrix-Vector Multiplication

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Mixed

  Multiplication (MVM) is used in all State-of-the-Art (SoA) neural network implementations. MVM performs Multiplication and Accumulation (MAC) of input and weight vectors. Single MAC operation in a Von Neumann architecture [Murmann ESSCIRC 2020] : • Data is fetched 4 times • Register access cost = 50fJ/byte • Total access cost = 200fJ/byte • System maximum efficiency = 10TOPS/W • The processing elements are integrated into the memory. • The energy efficiency upper limit is set by the ADC. • The energy cost of the ADC is shared with multiple MACs across an Accumulation Line (AL). • ADC-limited efficiency of 200TOPS/W (arrays > 100 MACs). • Multi-bit operations up to 16 states of conductance (4 bits) [Hsieh EDL 2021] • Process variations prevent reaching a higher number of bits. • Writing necessitates high current spikes (around 6µA) [Dao & Koch ICEIC 2020] • Inference cost of approximately 250fJ/MAC [Murmann TVLSI 2021] • A trade-off exists between energy consumption and variability. • Needs additional digital circuitry to combine all the line's results. Only adequate for a reduced number of lines (<5) but dominates the power consumption for higher numbers of bits [Murmann TVLSI 2021] • Take advantage of the time to charge the capacitor with ultra-small currents (thanks to recent CMOS technology) • Envisioned core energy consumption of 0.5fJ/MAC -4000TOPS/W (dominated by the ADC 10fJ/MAC -200 TOPS/W for an 8bit SoA ADC) • PWM period and the current level can be tuned to address many applications. Switched current sources offer great opportunities for configurable, lowenergy multi-bit MVM. Enabled by ultra-low currents in advanced CMOS technologies, this contribution proposes to trade off time to reduce energy, being specifically recommended for applications with throughput lower than 50GOPS, such as Key Word Spotting (KWS), Voice Activity Detection (VAD), or sound recognition. It allows reaching the theoretical ADC-limited consumption of 200TOPS/W for applications with low throughput requirements. This work was supported in part by the French National Research Agency under Grant ANR-18-CE24-0006-01 LEOPAR Comparison of computation methods for multi-bit MVM


Kévin Hérissé
email: kevin.herisse@junia.com

Benoit Larras

Bruno Stefanelli

Antoine Frappé

Andreas Kaiser

Mixed-Signal In-Memory Multi-bit Matrix-Vector Multiplication

The applications for artificial intelligence are wide and cover multiple domains including industry, health, home automation, consumer electronics, automotive, and smart cities. Application-specific integrated circuits performing tiny machine learning at ultra-low power and high accuracy are needed. The Von Neumann wall forces us to shift the processing elements closer to the memory to prevent data movement and therefore reduce energy consumption. Matrix-Vector Multiplication (MVM) can be achieved with many approaches that perform well with binary weights but not with multi-bit multiplications. This paper tries to highlight the advantages of using current sources to perform in-memory computing, improving further the energy consumption to perform multi-bit MVM.

Introduction

Matrix-Vector Multiplication (MVM) is used in all

State-of-the-Art (SoA) neural network implementations, such as Convolutional Neural Network (CNN), Recurrent Neural Network (RNN), or Long Short-Term Memory (LSTM).

MVM performs Multiplication and Accumulation (MAC) of input and weight vectors. In a Von Neumann architecture [START_REF] Murmann | Mixed-Signal Processing Opportunities for AI[END_REF], we can see that data is fetched 4 times from the memory to perform a MAC operation, as shown in Fig. 1. By estimating the register access cost at 50fJ/byte, the total access cost reaches 200fJ/byte and therefore clamps the system to maximum efficiency of 10TOPS/W (1 MAC = 2 Operations). This is known as the Von Neumann memory wall. In-Memory Computing (IMC) proposes to shift the processing elements inside the memory and therefore gain energy by reducing memory access. Mixed-signal IMC has been studied extensively, but a majority of research focuses on binary (or ternary) solutions. This paper addresses mixedsignal architectures for multi-bit operation.

Section 2 provides an overview of SoA mixed-signal multi-bit solutions for IMC. Section 3 evaluates the potential of switched current sources and their advantages, and Section 4 concludes the document. In this configuration, the ADC usually sets the upper limit in energy efficiency. However, the energy cost of the ADC is shared with multiple MACs across an AL. ADClimited efficiencies of 200TOPS/W are envisioned [START_REF] Murmann | Mixed-Signal Processing Opportunities for AI[END_REF] for large arrays (>100 MACs). This efficiency value can only be reached if the process of multiplication of input and weight vectors can be realized at a fraction of the ADC energy cost. The next subsections detail and compare several approaches for this multiplication operation.

Overview of Mixed-Signal IMC for MVM

ReRAM approach

In the resistive approach, memristors are used to perform the multiplication. The digital input vector is sent as an analog voltage across all the rows of the array. The weights are stored as conductance values of the memristors, modulating the current through the AL. This approach requires DACs with large output currents, and a current comparator, as highlighted in [START_REF] Cosemans | Towards 10000TOPS/W DNN Inference with Analog in-Memory Computing -A Circuit Blueprint, Device Options and Requirements[END_REF]. It is possible to perform multi-bit operations with this approach by using up to 16 states of conductance [START_REF] Hsieh | Four-Bits-Per-Memory One-Transistorand-Eight-Resistive-Random-Access-Memory (1T8R) Array[END_REF]. However, process variations prevent reaching a higher number of bits. Furthermore, in terms of energy dissipation, writing necessitates high current spikes (around 6µA in [START_REF] Dao | Memristor-based Reconfigurable Circuits: Challenges in Implementation[END_REF]), while inference costs approximately 250fJ/MAC, according to [START_REF] Murmann | Mixed-Signal Computing for Deep Neural Network Inference[END_REF]. It is to note that a trade-off exists between energy consumption and variability.

Switched Capacitor approach

In the capacitive approach, the multiplication is based on sharing charges or currents on a capacitive line, usually using XNOR gates for binary multiplication [START_REF] Knag | A 617-TOPS/W All-Digital Binary Neural Network Accelerator in 10-nm FinFET CMOS[END_REF]. Unit capacitors are charged according to each binary multiplication and charges are redistributed across all capacitors on an AL, resulting in a voltage to be converted by the ADC [START_REF] Bankman | An 8-bit[END_REF]. To perform a multi-bit operation with switched capacitors, [START_REF] Murmann | Mixed-Signal Computing for Deep Neural Network Inference[END_REF] shows a topology similar to a digital multiplier using one AL for each bit. However, this method needs additional circuitry to combine all the line's results, which is adequate for a reduced number of lines (<5) but dominates the power consumption for higher numbers of bits.

Switched Current Sources advantages

By using current sources to charge or discharge a capacitive line, multi-bit operations can be performed without additional circuits at low energy consumption. Current sources are distributed over the AL and are switched by the activation lines with Pulse Width Modulated (PWM) signals. Recent CMOS technologies allow driving small currents precisely enough to perform MVM with good accuracy. The main principle is to take advantage of the time to charge the capacitor with ultrasmall currents. This solution is specifically recommended for applications with a throughput lower than 50GOPS, such as Key Word Spotting (KWS), Voice Activity Detection (VAD), or sound recognition. For example, for an AL composed of 100 multi-bit MACs, a unitary 100pA current, and an AL capacitance equivalent to 100fF, approximately 5µs are needed to charge a capacitor to 0.5V (allowing for signed operation in a 1V supply). The equivalent energy consumption for the computation of the MVM is then calculated to 0.5fJ/MAC (4000 TOPS/W). Using this configuration, the system's core consumption is small enough so that the total energy will be eventually dominated by the ADC and reaches 10fJ/MAC (200TOPS/W) for an 8-bit SoA ADC [START_REF] Murmann | Mixed-Signal Processing Opportunities for AI[END_REF]. Furthermore, another advantage of switched current sources is the configurability, since the PWM period or the current level can be tuned to address many applications. Changing the current value modifies the gain of the multiplication so that the number of MACs per AL can be adjusted. [START_REF] Murmann | Mixed-Signal Computing for Deep Neural Network Inference[END_REF] also showed that mixed-signal IMC is suitable for bit resolutions up to 8 bits and that fully digital solutions should be considered for wider bit width.

Memristor

Conclusion

Switched current sources offer great opportunities for configurable, low-energy multi-bit MVM. Enabled by ultra-low currents in advanced CMOS technologies, this contribution proposes to trade-off time to reduce energy. It allows reaching the theoretical ADC-limited consumption of 200TOPS/W for applications with low throughput requirements.

Fig. 2

 2 Fig. 2 highlights the principle of mixed-signal IMC for MVM. A digital input vector is broadcast across multiple accumulation lines (AL) where it is multiplied by the weight vectors stored in the memory array. The accumulation is done in the current or voltage domains and an ADC converts the result at the bottom of each AL.

Figure 1 -Figure 2 -

 12 Figure 1 -MAC operation in a Von Neumann architecture

Table 1 -Comparison of computation methods for multi-bit MVM

 1

			Switched	Switched Current
			Capacitor	Source
	Configurability	-	-	+
	Technology	Non-standard (CMOS + RRAM)	Standard CMOS Standard CMOS
	Typical energy efficiency	21.9 TOPS/W [8] 10 TOPS/W [5]	200 TOPS/W
	Max number of bits	3 -4	8	8