Multi-objective Multidisciplinary Design Optimization approach for partially Reusable Launch Vehicle design

stage may be investigated. One option that may be explored [START_REF] Vila | Weighting options for the next generation of Ariane launchers[END_REF] is to design the launch vehicle first stage for both reusable and expendable uses, in order to increase the flexibility and adaptability to different target missions.

Following this strategy, this paper proposes a multi-level Multidisciplinary Design Optimization (MDO) approach to design aerospace vehicles addressing multi-mission problems. This category of design problems can be formulated as multi-objective MDO problems which involve numerous evaluations of MultiDisciplinary Analysis (MDA) used to assess the performance of the vehicles. This MDA involves computationally expensive coupled models (aerodynamics, propulsion, structural design, trajectory, etc.) [START_REF] Balesdent | Advanced space vehicle design taking into account multidisciplinary couplings and mixed epistemic / aleatory uncertainties[END_REF].

Several multi-objective MDO approaches for launch vehicle design have been investigated in the literature [START_REF] Castellini | Comparative analysis of global techniques for performance and design optimization of launchers[END_REF][START_REF] Fazeley | A multi-objective, multidisciplinary design optimization methodology for the conceptual design of a spacecraft bi-propellant propulsion system[END_REF][START_REF] Fujikawa | Multi-objective, multidisciplinary design optimization of TSTO space planes with RBCC engines[END_REF].

Castellini et al. [START_REF] Castellini | Comparative analysis of global techniques for performance and design optimization of launchers[END_REF] proposed a comparison of seven population-based algorithms (e.g. NSGA-II [START_REF] Deb | A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II[END_REF], MOPSO [START_REF] Coello Coello | MOPSO: a proposal for multiple objective particle swarm optimization[END_REF], PAES [START_REF] Knowles | Approximating the nondominated front using the Pareto archived evolution strategy[END_REF]) dealing with multi-objective problems and compared them on expendable launch vehicle design problems. Two types of problem have been considered, ascent trajectory optimization with a fixed launcher design and launcher design optimization (with optimization of architectures: number of boosters, type of propulsion, type of engine cycle, etc.).

In this work, a MultiDiscipline Feasible (MDF) formulation [START_REF] Balling | Optimization of coupled systems-A critical overview of approaches[END_REF] is carried out to couple the different disciplines and conceptual design models are used. The two considered objectives are the sum of square errors on semiaxis, eccentricity and inclination and the payload mass (or Gross Lift-Off Weight -GLOW). Fazeley et al. [START_REF] Fazeley | A multi-objective, multidisciplinary design optimization methodology for the conceptual design of a spacecraft bi-propellant propulsion system[END_REF] proposed to compare MDF and Collaborative Optimization (CO) for the multi-objective optimization of a bi-propellant space propulsion system design for an expendable launcher. The two objectives for the design of the propulsion systems are the total wet mass and the total impulse. NSGA-II is used to perform the multi-objective optimization. In comparison, for this specific problem, the authors concluded that MDF requires less calls to the disciplines than CO. Fujikawa et al. [START_REF] Fujikawa | Multi-objective, multidisciplinary design optimization of TSTO space planes with RBCC engines[END_REF] performed a conceptual design study for a Two-Stage-To-Orbit space plane with ethanol-fueled rocket-based combined cycle engine using a multi-objective MDO approach. An All-At-Once formulation is employed in the proposed process. Three objectives are considered for the problem: the payload mass, the GLOW and the take-off velocity. At each iteration, the multi-objective problem is transformed into its relevant single-objective problem via min-max goal programming, and is subsequently solved using a gradient-based optimization method (SQP) [START_REF] Fujikawa | Multi-Objective, High-Accuracy Trajectory Optimization Using SQP and Sampling on the Manifold[END_REF]. Kosugi et al. [START_REF] Kosugi | Multidisciplinary and multi-objective design exploration methodology for conceptual design of a hybrid rocket[END_REF] used a multi-objective genetic algorithm to design a hybrid rocket. The two considered objectives are the GLOW and the maximal reached altitude by the sounding rocket and a MDF formulation is carried out.

To the best of the author knowledge, most of the existing papers in the literature focus on multi-objective problems for a single launch vehicle using population-based algorithms. In contrast, this paper is focused on the design of a family of launchers for different missions sharing design commonalities (one reusable mission dedicated to Sun Synchronous Orbit and one expendable mission to Geostationary Transfer Orbit) using Bayesian Optimization to account for the computational cost associated with the discipline simulations.

In general, due to the repeated evaluations of the MDA, multi-objective MDO solving using a classical optimization algorithm such as evolutionary approaches (e.g., NSGA-II [START_REF] Deb | A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II[END_REF], OMOPSO [START_REF] Sierra | Improving PSO-based multi-objective optimization using crowding, mutation and epsilon-dominance[END_REF], PAES [START_REF] Knowles | Approximating the nondominated front using the Pareto archived evolution strategy[END_REF]) is intractable. In this paper, an alternative technique is used based on Bayesian multi-objective optimization. Bayesian optimization algorithms are employed to deal with computationally intensive black-box function problems. They rely on surrogate models to build a statistical relationship between the input design variables and the quantities of interest (objective and constraint functions given by the MDA), to predict these latters using a limited number of exact MDA simulations. The evaluation cost of the surrogate models is much cheaper compared to the exact MDA. Bayesian optimization methods have been extended to solve such multi-objective problems [START_REF] Beume | SMS-EMOA: Multiobjective selection based on dominated hypervolume[END_REF]. Specific infill criteria such as the Expected HyperVolume Improvement [START_REF] Svenson | Multiobjective optimization of expensive-to-evaluate deterministic computer simulator models[END_REF][START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF] can be used to determine the design candidates for which the exact MDA needs to be evaluated.

These candidates are added to the dataset and the surrogate models are updated. This process is repeated until a given criterion is reached. In this paper, a surrogate model-based MDO approach using a multi-level decomposition and relying on Gaussian Processes (GP) is proposed to estimate with a reduced computational cost the optimal Pareto Front for multi-mission aerospace vehicle design problems. The novelty of the proposed approach relies in the decomposition of the multi-mission problem into a multi-level MDO formulation involving efficient optimization techniques.

The proposed multi-level MDO formulation is illustrated on the design of a family of conceptual Two-Stage-To-Orbit launch vehicles. The multi-mission problem considers two missions: a reusable configuration for a Sun Synchronous Orbit (SSO) with a medium payload range and recovery of the first stage using a gliding-back strategy; and an expendable configuration for a payload injected into a Geostationary Transfer Orbit (GTO). Gliding-back strategy (also known as deadleaf) belongs to the family of Vertical Take-Off and Horizontal Landing approaches and is an alternative strategy to toss-back configuration by taking advantage of aeronautical technologies (e.g. addition of lifting surfaces).

The following of the paper is organized as follows. In Section 2, an overview of Multi-objective Bayesian optimization is presented along with Gaussian Process principles and infill criterion description. In Section 3, the application of the proposed MDO methodology to the design of multi-mission reusable launch vehicles is described with a presentation of the expendable and reusable missions along with the vehicle concept and design requirements. Then, the multi-level MDO formulation for the multi-mission problems is introduced and the disciplinary models are presented with the MultiDisciplinary Analysis. In Section 4, the results of the proposed multi-objective MDO process on the launch vehicle design problem are presented and compared to reference approaches.

III. Multi-objective Bayesian optimization

The design of a reusable launch vehicle for multi-mission can be described as a multi-objective MDO problem. A general multi-objective problem characterized by n objective functions and q constraints may be formulated as:

min f(z) = [ f 1 (z), • • • , f n (z)] (1)
w.r.t. z ∈ R m s.t. g i (z) ≤ 0, i = 1, . . . , q (2)
z min ≤ z ≤ z max (3)
with z the vector of design variables of dimension m, f(•) the objective function vector of size n and g i (•) the i th constraint function of the vector of constraints g(•) of size q. Evolutionary algorithms such as NSGA-II (Non-dominated Sorting Genetic Algorithm II) [START_REF] Deb | A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II[END_REF] or SMPSO (Speed-constrained Multi-objective Particle Swarm Optimization) [START_REF] Nebro | Smpso: A new pso-based metaheuristic for multi-objective optimization[END_REF] are classically used to solve multi-objective optimization problems. However, due to the computationally intensive nature of the MultiDisciplinary Analysis involved in launch vehicle design, this type of algorithms is not suited for such problems.

In this work, an alternative technique is used based on Bayesian multi-objective optimization. It relies on surrogate models to build a statistical relationship between the input design variables and the quantities of interest (objective and constraint functions), to predict these latter using a limited number of data also called Design of Experiments (DoE).

The evaluation cost of the surrogate models is much cheaper compared to the exact functions. In this paper, Efficient Global Optimization (EGO) developed by Jones et al. [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] is considered. It relies on Kriging surrogate model [START_REF] Rasmussen | Gaussian processes in machine learning[END_REF] which is based on the Gaussian Process (GP) theory. The main advantage of Kriging compared to alternative surrogate modeling techniques is that in addition to the prediction, it provides uncertainty estimation of the model response. Based on these two information, an infill criterion (such as the Expected Improvement adapted to multi-objective problems) is used to iteratively add the most promising candidate to the existing DoE. This solution is then evaluated on the exact functions (MDA providing objective and constraint function values) and the GPs are updated and so on, until a stopping criterion is reached.

A. Gaussian Process (GP)

Let us initially consider a single objective unconstrained optimization problem. A surrogate model based on Gaussian Process can be built from a set of samples (DoE) of size N, Z N = z (1) , . . . , z (N ) the input data set (z ∈ R m ) and the corresponding objective function responses Y N = y (1) = f z (1) , . . . , y (N ) = f z (N ) . Then, it is possible to use this surrogate model to predict the exact function response f (•) at a new point without evaluating it. The advantage of this approach is relative to the computational cost reduction. Indeed, instead of evaluating the expensive black-box function f (•), the surrogate model (which is much cheaper to evaluate) is used. The interest of the GP surrogate model is that it provides a variance estimation of the prediction in addition to this latter, which makes it suitable in a global optimization framework (Figure 1).

A GP is used to describe a distribution over functions, it corresponds to a collection of infinite random variables, any finite number of which has a joint Gaussian distribution. A GP is characterized by its mean and covariance functions. In a GP regression, a GP prior is placed on the unobserved function f (•) (also referred as latent function) using a prior covariance function k Θ (z, z ) that depends on several hyper-parameters Θ and mean function µ(•). As the trend of the response is a priori unknown a constant mean function µ is assumed (then GP is named ordinary Kriging). Therefore, the GP may be written as f (z) ∼ N µ, k Θ (z, z ) and has a multivariate distribution on any finite subset of variables, in particular on the DoE Z N (noted f N ), f N |Z N ∼ N 1µ, K Θ N N where K Θ N N is the covariance matrix constructed from the parametrized covariance function k Θ (•, •) on Z N (in the following the dependence on Θ is dropped for notation simplicity). The choice of the covariance function determines our prior assumptions of the function to be modeled. A Gaussian noise variance is considered, such that the relationship between the latent function values f (Z N ) and the observed responses Y N is given by: p (y|f N ) = N y|f N , σ 2 I , with p(•) the density distribution of y given f N . This Gaussian noise is necessary when the objective function or the constraint responses are non deterministic. The marginal likelihood is obtained by integrating out the latent function f (•):

p (y|Z N , Θ) = N y| µ, K N N + σ 2 I (4)

Fig. 1 Example of Gaussian process prediction and associated confidence interval

In order to train the GP, it is possible to maximize the negative log marginal likelihood to find the optimal values of the hyper-parameters Θ, µ, σ. Gradient-based optimizers may be used to minimize the negative log marginal likelihood in order to determine the hyperparameter values of the trained GP.

After the training, the prediction at a new point z * is made by using the conditional properties of a multivariate normal distribution (Figure 1):

p (y * |z * , Z N , Y N , Θ) = N y * | ŷ * , ŝ * 2 (5)
with ŷ * , ŝ * 2 the mean prediction and the associated variance given by:

ŷ * = µ + k T z * K N N + σ 2 I -1 (y -1µ) (6) ŝ * 2 = k z * ,z * -k T z * K N N + σ 2 I -1 k z * + σ 2 (7)
where

k z * ,z * = k(z * , z * ) and k z * = k z (i) , z * i=1,..., N
An example of a 1D function and a GP built based on four observations (samples) is represented in Figure 1. The confidence interval is based on the associated GP variance (±3 ŝ * ). The variance is null at the observation and increases as the distance from an existing data sample increases (here σ = 0). GP is very interesting in the complex design because it provides the design with both the prediction of the model and its estimated uncertainty. This possibility to quantify the surrogate model prediction uncertainty is exploited in optimization to balance between global search exploration and exploitation. To solve an optimization problem with Bayesian optimization, an infill criterion such as the Expected Improvement [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] is used to iteratively add the most promising candidate to the existing DoE. This solution is then evaluated on the exact objective function and the GP is updated. This process is repeated until a stopping criterion is reached. Bayesian optimization has been extended to solve multi-objective problems as discussed in the next Section.

B. Multi-objective Infill Sampling Criteria

For multi-objective problems, different approaches have been proposed to select the sample candidates to be added to the current DoE such as aggregation-based techniques [START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with Gaussian process model[END_REF] or domination-based techniques [START_REF] Svenson | Multiobjective optimization of expensive-to-evaluate deterministic computer simulator models[END_REF][START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF]. In this study the approach proposed by Emmerich et al. [START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF] based on the Expected HyperVolume Improvement (EHVI) is used.

The notion of EHVI is an adaptation of the Expected Improvement (EI) to the multi-objective case. Consider an unconstrained multi-objective optimization problem characterized by n objectives, y = f(z) = [y 1 , . . . , y n ] stands for the objective function outputs. Considering a DoE of size N, Z N = z (1) , . . . , z (N ) the input data set and the corresponding objective function responses Y N = y (1) = f z (1) , . . . ,

y (N ) = f z (N ) . Let V = y ∈ R n | y L ≤ y ≤ y U be
H Y N = y ∈ V | ∃i ∈ {1, . . . , N }, f z (i) ≺ y (8)
H Y N is a subset of V whose points are dominated by the DoE responses Y N (the symbol ≺ stands for Pareto domination).

Considering a new candidate z (N +1) , since H Y N ⊆ H Y N +1 , the hypervolume improvement by adding this candidate is given by:

I N z (N +1) = |H Y N +1 | -|H Y N |.
EHVI Y N (z) = E |H Y N +1 | -|H Y N | = ∫ V\H Y N P Ŷ(z) ≺ p dp (9)
In the constrained case, a constrained infill criterion is considered which combines EHVI with the probability of feasibility [START_REF] Schonlau | Global versus local search in constrained optimization of computer models[END_REF] or the expected violation of the constraints [START_REF] Audet | A surrogate-model-based method for constrained optimization[END_REF]. Schonlau [START_REF] Schonlau | Global versus local search in constrained optimization of computer models[END_REF] proposed to integrate a measure of the feasibility of a candidate point into the unconstrained infill criterion (initially for single objective problems but extended to multi-objective problems here). The concept of probability of feasibility is defined as the probability that a candidate satisfies all the constraints (independent):

P f (z) = q i=1 Φ 0 -ĝi (z) σgi (z) (10)
with ĝi (•) the GP mean prediction of the constraint function g i (•), σ2 gi (•) the prediction variance and Φ(•) the cumulative distribution of the normal law. By multiplying this probability of feasibility to the Expected Hypervolume Improvement, the magnitude of the resulting quantity tends to zero in the design space areas where there is a low likelihood of the constraint feasibility and it tends to the EHVI value where there is a high likelihood of the constraint feasibility.

Fig. 3 Steps of Multi-objective EGO (MO-EGO)

At each iteration of the Multi-Objective EGO (MO-EGO), this infill criterion C(z) = E HV I(z) × P f (z) is maximized in order to identify the most promising candidate. This latter is evaluated on the exact objective and constraint functions and added to the DoE. All the surrogate models are updated and a new iteration can start if the convergence criterion is not satisfied. In general, the algorithm is stopped after a certain number of exact function evaluations depending on the available computational budget. The different steps of MO-EGO are summarized in Figure 3.

IV. Application to reusable launch vehicle design

In order to illustrate the proposed method to solve multi-mission launch vehicle design problem, an application to reusable launch vehicle design is presented in the following Section.

A. Description of vehicle concept, mission and design requirements

In function of different assumptions on addressable institutional and commercial satellite markets, various strategies for the reusability of the first stage of launch vehicles may be investigated. One option that may be explored [START_REF] Vila | Weighting options for the next generation of Ariane launchers[END_REF] is to design the first stage for both reusable and expendable configurations in order to address different missions. Following this strategy, two baseline missions are considered in this paper:

• SSO mission, with recovery, refurbishment and reuse of the launch vehicle first stage,

• GTO mission, in expendable mode with no recovery of the vehicle first stage. The proposed launch vehicle architecture should be able to perform these two missions. In order to identify the possible optimal architectures, the design of such a vehicle may be obtained through the solving of a multi-objective problem and determining the optimal Pareto front.

The considered launchers use lost cost LOx/CH4 engines in both stages with adaptation of the nozzle for the first and the second stages [START_REF] Iannetti | Prometheus, a LOx/LCH4 reusable rocket engine[END_REF]. For the SSO mission, all the additional subsystems that are added to the vehicle first stage to enable recovery are part of a reusability kit. This kit is removed for the expendable GTO mission and can be added to another first stage for a future reusable mission. Moreover, for the expendable mission, two solid rocket boosters are added to the central core. For the SSO mission, a gliding-back strategy is used. The first stage is equipped with several subsystems such as lifting surfaces (wings and control surfaces) which are added to the stage to enable a glided return to the launch site (Figures 6,7).

The sequences of the gliding-back strategy are presented in Figure 7. After the lift-off (step 1) and the separation of the first stage and second stage (step 2), the first stage is turned around and the main engines (one or several depending on the return strategy) are reignited (step 3) in order to inverse the horizontal velocity. This operation is performed usually in high altitude (above 50km) with limited drag. Then, a reentry into the atmosphere is performed (step 4) and The reusability kit for the gliding-back version of the 1 st stage is composed of: lifting surfaces (main wings and canards), a fairing, a skirt between the fairing and the upper tank, a vertical tail, and additional power and avionics.

Moreover, compared to classical expendable stage, the structures are reinforced to account for additional loads during the rocket engine boost to inverse the horizontal velocity, the reentry and the aerodynamic pull up maneuver, especially with respect to transverse loads which are not usual for such a type of vehicles.

Single-level, multi-objective MDO formulation

In order to design a reusable launch vehicle for the two missions described in the previous Section, a multi-objective MDO problem is set. The reusable and expendable architectures share the same first and second stages, the difference between the two is the absence of the reusability kit for the expendable mode, the presence of solid rocket boosters for the expendable strategy and the possibility to use all the propellants for the payload injection. The optimization objective for the reusable vehicle is to minimize its Gross Lift-Off Weight whereas for the expendable launcher, the optimization objective is to maximize the mass of the injected payload (Figure 9). Due to the fact that the two missions share the same first and second stages (except for the reusability kit) these two stages are sized with respect to the mission that leads to the maximal trajectory loads. Therefore a MDO coupling exists between the two missions and consequently design problems (Figure 8).

The propulsion, geometry and structural sizing of the main core are commune disciplines for the two configurations.

A specific module for the structural sizing of the reusability kit is present in the reusable configuration. A module called "maximal trajectory load" is in charge to compare the loads undertaken by the two vehicles during the trajectory and to take the maximum value in order to size the vehicle with respect to the most constrained one. In practice, in order to avoid loops between the MDA of the two missions, the coupling variables between the two configurations are removed and controlled at the system-level by the optimizer along with the design variables. The single-level, multi-objective

z min ≤ z ≤ z max [START_REF] Denis | ONERA's aerodynamic prediction code-MISSILE[END_REF] with z ∈ R 38 the design variable vector.

mathb f c r euse and c ex p are the control vector for the trajectory including pitch, angle of attack and yaw controls during the different phases of the trajectories (ascent and return for the reusable mission). S wing and λ wing are the surface area and the aspect ratio of the main wings for the reusable launcher, while β M prop1 is the allocation of the propellant mass for the first stage M pr op1 that is not used for the ascent trajectory and dedicated to the rocket engine boost to inverse the horizontal velocity during the return to launch site phase. Several inequality constraints are considered in this multi-objective problem, including the altitudes of apogee and perigee for the two missions (h apogee , h perigee ), the maximal acceptable loads (axial (NX), transverse (NZ), dynamic pressure (Pdyn) and flux). Moreover, a constraint representing the distance to the landing site d sitemax (both in altitude and range) is added in the optimization problem for the return trajectory. All the constraints are specified with respect to a maximal exceedance threshold. M prop1 , M pr op2 , N X d , N Z d and Pdyn d are common to the two vehicles. The loads are considered for the sizing of the stages (see Section IV.C.1 for more details).

Multi-level and multi-objective MDO formulation

Due to the design space dimension (z ∈ R 38 ), the non-convexity and the proportion of the feasible region over the entire design space for such multi-objective problems, it is complicated to solve with a single-level (see Section V.B for the application of the single-level MDO formulation). Indeed, the injection into orbit (altitudes of apogee and perigee, velocities, etc.) and the return trajectory constraints (final distance to site, admissible loads, etc.) are difficult to satisfy.

In order to efficiently solve this design problem, a multi-level MDO process is proposed (Figure 10). At the system-level, a Bayesian optimizer controls the decision variables that are shared between the two configurations, whereas, at the subsystem-level, two MDO problems are solved, one for the expendable configuration and one for the reusable launcher in order to facilitate the satisfaction of the constraints of the system-level. The two subsystem optimizers control the design variables specific to the considered configurations. The multi-level approach enables to distribute the problem complexity over different dedicated subproblem optimizations.

The proposed multi-level formulation is given by: System-level:

min f(z) = GLOW reuse (z * ), -M paylaod ex p (z * ) (27)

Fig. 10 Multi-level MDO formulation for multi-mission launch vehicle design

This formulation decomposes the global problem into three optimization problems that are easier to solve. At the system-level, the propellant masses for the two launchers are controlled as well as the sizing loads (axial, transverse and dynamic pressure) that must be shared by the two vehicles. The first and second stages are sized with respect to the most important loads of the two missions. Indeed, as the first stage will be used both in expendable and reusable configurations, it must be sized for these two types of trajectories. Therefore, by controlling these variables at the system-level, it ensures that the dry mass (e.g. tanks, thrust frame, skirt, intertanks, interstage, engines, subsystems) of the first and second stages are identical for the two missions (except for the reusability kit that is removed for the expendable mission). In addition, by controlling the sizing loads at the system-level, it enables to decouple the trajectory and the sizing disciplines and therefore to avoid a Fixed Point Iteration (FPI) between the disciplines to solve the MDA (see next Section). At the system-level, two constraints (J cons e x p , J cons r eu s e ) are handled to aggregate the constraints of the subsystem-levels to ensure that the found solutions at the subsystem-level are feasible with respect to the sizing constraints and the couplings between the trajectory and sizing. At the subsystem-level, each optimization problem controls the specific decision variables associated to the corresponding mission. δ M pr o p1 ∈ [0.9, 1.0] and δ M pr o p2 ∈ [0.9, 1.0] are auxiliary propellant mass variables to provide some flexibility in the optimization problem solving and enabling to not use the entire propellant masses provided by the system-level optimizer if not required. This multi-level MDO process enables to find candidate solutions satisfying the constraints, however it increases the computational cost at each iteration at the system-level, justifying the need for Bayesian optimization and the use of EGO derivation for multi-objective problems. Indeed, an iteration at the system-level involves two optimization problem solving at the subsystem-level. Even if these subsystem problems are easier to solve than the complete design problem, they still require important computational calculations to find the optimal subsystem design z corresponding to the system-level design variable vector value z * .

In order to estimate the objective function and the constraints for the two missions, a MultiDisciplinary Analysis is carried out with the discipline models that are described in the next Section.

C. MultiDisciplinary Analysis and disciplinary models

In order to assess the performance of the launch vehicle for the two different missions, a MultiDisciplinary Analysis (MDA) is used. The different disciplinary models (aerodynamics, propulsion, geometry, sizing, trajectory) are coupled through a Fixed Point Iteration method [START_REF] Ortega | Stability of difference equations and convergence of iterative processes[END_REF] (Figure 11).

Fig. 11 Illustrative MultiDisciplinary Analysis for launch vehicle design

In the proposed approach, in order to avoid FPI, the feedback coupling variables between the trajectory discipline and the sizing discipline are removed and controlled at the system-level (Figure 11). The feedback couplings correspond to the maximal loads undertaken during the trajectory.

Geometry and sizing

In this test case, the geometry and sizing discipline aims at estimating the mass budget of the launcher and its geometry. The 1 st and the 2 nd stages are LOx/CH4 stages. The dry mass of the stages is the sum of the masses of the structural masses (tanks, skirt, intertanks, thrust frame, interstages), turbopumps, combustion chamber, nozzle, pressurization system and avionics. The mass and geometry models are derived from the engineering models for the conceptual design of launch vehicle developed by Castellini [START_REF] Castellini | Multidisciplinary design optimization for expendable launch vehicles[END_REF]. Fig. 12 GLOW computations with all the involved components, extracted from [START_REF] Castellini | Multidisciplinary design optimization for expendable launch vehicles[END_REF] Several outputs from other disciplines are required as inputs to geometry and structural sizing to compute the GLOW:

GLOW

• Propulsion outputs: thrust, specific impulse, solid propellant grain geometry, chamber pressure, etc.

• Geometry outputs: geometry of all the different parts of the launch vehicle (tank volume, stage length, launcher surface, engine size, etc.)

• Trajectory outputs: maximal axial and transverse loads and maximal dynamic pressure endured by the launch vehicle during the flight.

For the reusable mission and especially the recovering subsystems (e.g. lifting and control surfaces) dedicated engineering models [START_REF] Rohrschneider | A comparison of modern and historic mass estimating relationships on a two stage to orbit launch vehicle[END_REF] are used to estimate their masses. An interdisciplinary coupling between the sizing discipline and the trajectory is required to model the dependencies between the dry mass and the loads undergone during the flight. Three types of loads are considered for the sizing of the stages: the axial (N X d ), transverse (N Z d ) loads and the dynamic pressure (Pdyn d ). These loads are provided by the trajectory simulation to the sizing discipline. These three variables stand for the feedback couplings in the MDA. In order to avoid the use of FPI, these three variables are controlled at the system-level as explained in the previous Section. Taking into account these loads is essential for the sizing of the launch vehicle stage. Indeed, the reuse of a launch vehicle with a gliding-back strategy leads to more important transverse load that conventional expendable launchers. These loads have to be taken into account in order to ensure that the vehicle will be able to perform the mission. The increase in the undertaken loads leads to an increase in the stage dry mass and the structural index (ratio of dry mass to propellant mass). In Figure 13, an illustration of the variation of the structural index with the maximal transversal load encountered during the flight is represented.

Fig. 13 Variation of the structural index for stage 1 as a function of the maximal transversal load for a given propellant mass

Aerodynamics

The aerodynamics discipline consists in computing the aerodynamics coefficients such as the drag and lift coefficients required to compute the aerodynamics loads during the launcher atmospheric flight. The calculations of the drag and lift coefficients are based on the ONERA code MISSILE [START_REF] Denis | ONERA's aerodynamic prediction code-MISSILE[END_REF] which relies on simplified aerodynamics theory and on an experimental data base to determine the aerodynamics forces and coefficients of complex launcher geometries.

Drag and lift coefficient tables as a function of the Mach number and the angle of attack are directly given to the trajectory discipline allowing to remove the feedback loop between the trajectory and the aerodynamics. This model is generally sufficient in the early design studies. In order to decrease the computational cost associated with the MISSILE execution, a Neural Network relying on a Multi-Layer Perceptron [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] is used to generate the aerodynamics coefficients from an original database. In Figure 14, a comparison between the Neural Network prediction and the exact MISSILE aerodynamics coefficients is presented for the lift to drag ratio of the reusable configuration.

Trajectory

In general, the trajectory discipline consists in solving an optimization problem. The objective of ascent trajectory optimization is to minimize the distance between the injection point and the given target. The trajectory model used the 20 following 3D dynamics equations, written in an Earth-centered, rotating Earth referential [START_REF] Castellini | Multidisciplinary design optimization for expendable launch vehicles[END_REF].

r = V sin(γ) (31)
V = T cos(θ -γ) -D m -g(r) sin(γ) + ω 2 E r cos(φ)(sin(γ) cos(φ) (32)
cos(γ) sin(φ) cos(ψ))

γ = [L + T sin(θ -γ)] cos(µ) mV + V r - g(r) V cos(γ) + 2ω E sin(ψ) cos(φ) (33)
+ ω 2 E r cos(φ)(cos(γ) cos(φ) + sin(γ) sin(φ) cos(ψ)) V λ = V cos(γ) sin(ψ) r cos(φ) (34)
φ = V cos(γ) cos(ψ) r (35)
ψ = [L + T sin(θ -γ)] sin(µ) mV cos(γ) + V cos(γ) sin(ψ) tan(φ) r + 2ω E (sin(φ) (36)
-cos(ψ) cos(φ) tan(γ)) + ω 2 E r sin(φ) cos(φ) sin(ψ) V cos(γ) m = -q (37)
Table 1 Trajectory variables and parameters. The ascent trajectory optimization consists of different flight phases: vertical take-off, pitch-over, gravity turn, bilinear tangent law, upper stage's coast and circularization burn for the SSO or injection into the GTO for the expendable mission. The vertical take-off phase enables to leave the launch pad before doing a pitch-over maneuver constituted depending on the target inclination and the current latitude (declination) [START_REF] Castellini | Multidisciplinary design optimization for expendable launch vehicles[END_REF]. During the gravity turn phase, the yaw is constrained to follow the target inclination law and the pitch is equal to the flight path angle. Once the exoatmospheric conditions have been reached, it is possible to control the pitch with a piece-wise linear interpolation of optimizable nodal value θ gt (t). For the second stage's flight, a bilinear tangent law [START_REF] Castellini | Multidisciplinary design optimization for expendable launch vehicles[END_REF] is assumed for the pitch angle and its parameters are decision variables. Finally, upper stage's coast and circularization burn are considered since the considered upper stage engine is reignitable. The engine is shut down once the current orbit apogee matches the target apogee. Then, the propellant mass required to circularize the orbit is determined based on the final conditions of velocity and inclination. During the ascent trajectory, several constraints are taken into account such as the maximal axial loads, the maximal dynamic pressure, the maximal flux, the maximal angle of attack during the atmospheric flight and the payload injection conditions on the perigee and apogee.

Symbol

All the previous described trajectory parameters are decision variables noted θ ex p and θ reuse for the two missions and are optimized along with launcher sizing variables.

For the return trajectory of the gliding-back strategy, the objective of the optimization problem is to minimize the distance between the actual landing point of the vehicle and the landing site under some maximal admissible load constraints (e.g. axial and transverse load factors, dynamic pressure, heat flux). The decision variables for the return phase are β M pr op1 (allocation of the propellant mass for the engine boost), a piece-wise linear interpolation of optimizable nodal value α(t) for the guidance during the return phase.

V. Multi-objective and Multidisciplinary Design Optimization: results and discussions

A. Multi-level, multi-mission MDO approach settings

The proposed multi-level and multi-mission MDO approach requires to optimize the subsystem problems which control the design variables specific to each mission. In order to optimize the subsystem problems, Covariance Matrix Adaptation -Evolution Strategy (CMA-ES) optimization algorithm [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF] is used due to the non-convex objective and constraint functions and the presence of numerous minima. Indeed, by using a classical SQP algorithm with different initializations, it converges to different local minima leading to a non robust convergence of the lower level (Figure 16). On the contrary, by exploring the design space, CMA-ES succeeds to find a better optimum than a gradient-based approach and is more robust with respect to the initialization (that is unknown and set in the middle of the design space). CMA-ES is parallelized in order to evaluate each candidate solution with a multiprocessing approach to increase the computational efficiency of the subsystem problem solving. A penalization approach is carried out to control the optimization problem constraints at the subproblem-level. To evaluate the efficiency of the proposed optimization strategy, a comparison with a multi-objective evolutionary algorithm [START_REF] Sierra | Improving PSO-based multi-objective optimization using crowding, mutation and epsilon-dominance[END_REF] named OMOPSO is carried out. Particle Swarm Optimization (PSO) [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF] is a bio-inspired meta-heuristic mimicking the social behavior of bird flocking or fish schooling. In PSO, the solutions are represented as particles inside a swarm, with a certain position and velocity. The position of a particle at a given generation defines the value of the decision variables at this generation, and it is calculated using a so-called velocity equation. This equation expresses the variation in the position between two generations. It is calculated as a balance between the best position of the particle and the global best position of all the particles. In the adaptation of a PSO to the multi-objective case, the main difference is the notion of the best individual. While in the single objective case only one best value may be obtained, in multi-objective problems, several equally good solutions may be obtained. A natural way to extend this notion is to consider every non-dominated solution over the generations as a leader. In the process of updating the leaders, the particles that are non dominated are kept in the set of leaders. In OMOPSO, the idea to choose for each particle a leader using a binary tournament selection based on the crowding distance. That is for each particle, two solutions from the set of leaders are randomly selected, then the solution with the greatest crowding distance is chosen. For more details on OMOPSO see [START_REF] Sierra | Improving PSO-based multi-objective optimization using crowding, mutation and epsilon-dominance[END_REF]. GPflow is used to model GP in the proposed multi-objective Bayesian Optimization approach [START_REF] Matthews | GPflow: A Gaussian process library using TensorFlow[END_REF]. From the ten repetitions, the proposed multi-level MDO approach enables to find a better Pareto front after 220 evaluations (20 samples for the initial DoE and 200 multi-objective EGO iterations) of the exact MDA compared to OMOPSO. Moreover, the proposed methodology is more robust to the initial DoE samples than OMOPSO. Indeed, over the ten repetitions, the multi-level MDO formulation ends at the same value of hypervolume (around 0.287 with a standard deviation of 0.00076) whereas OMOPSO presents a large variability in the final hypervolume (around 0.285 with a standard deviation of 0.0029) as illustrated in Figure 26. Moreover, the proposed multi-level MDO approach reaches the optimum hypervolume faster than OMOPSO on average, indeed, after 100 iterations, the mean value of hypervolume is 0.286 for EGO whereas it is only 0.281 for OMOPSO. Therefore, the proposed algorithm could be stopped earlier and save computational time. These results are presented in Table 2.

In order to further analyze the results, the initial Pareto front determined from the initial DoE of size 20 for the first repetition is represented in Figure 22. The proposed approach succeeded to add design points to improve the hypervolume (Figure 26) and the final Pareto front (Figure 23). The final Pareto front of EGO is better than the OMOPSO one. Moreover, for the Pareto front, no OMOPSO point dominates EGO points. From Figure 26, it can be seen that the hypervolume of EGO is always better than of OMOPSO even if the gap tends to decrease with the iteration there is still an important difference at the end of the optimization process. The other repetitions present similar trends.

Figure 27 illustrates the evolution of the EGO Pareto front through the iterations (at iteration 50, 100, 150 and 200). The EGO Pareto front spreads and moves forward to tend to an optimal Pareto front (which is unknown). The hypervolumes increases by adding appropriate design points which improve the Pareto front. In Figure 24, the evolution of the Pareto front for OMOPSO is depicted.

Compared to OMOPSO, multi-level MDO approach is able to provide a better Pareto front with a limited number of exact function evaluations. Moreover, using a multi-level approach allows to overcome the difficulty of constraint satisfaction in reusable launch vehicle design. The final Pareto front provides key elements to make an informed decision in terms of multi-mission launch vehicle design. Indeed, the proposed design process offers a quantification of the repercussions on the reusable GLOW for a SSO mission due to an increase of the payload in expendable mission to GTO orbit. From these graphs, it can be seen that in order to inject the maximal payload mass into orbit for the expendable mission In order to illustrate the difference between Pareto front candidates, two alternative solutions on the Pareto front are studied corresponding to the maximal payload mass and the minimal GLOW (Figure 28). configuration. These trajectories present strong similarities as they differ principally by the amount of propellants for the two stages of around 40t. However, differences can be seen during the re-entry in terms of mass needed to bring back the first stage at the landing site and in terms of maximal altitude before the re-entry into the atmosphere, resulting in differences in terms of mechanical and thermal solicitations (e.g. encountered loads, dynamic pressure, heat flux).

VI. Conclusions

In this paper, a multi-level multi-mission decoupled Multidisciplinary Design Optimization (MDO) formulation has been proposed in order to design a family of reusable launchers to address two specific missions: a medium range mass payload into Sun-Synchronous Orbit with return to the launch site of the first stage using a gliding-back architecture, and a medium range mass payload into Geostationary Transfer Orbit with an expendable architecture. The proposed approach relies on multi-objective Efficient Global Optimization (EGO) to efficiently solve the multi-mission problem and to limit the number of calls to the exact MultiDisciplinary Analysis thanks to the use of Gaussian Processes. The multi-mission aspect introduces couplings between the different missions in order to ensure the appropriate structural sizing of the launch vehicle. The multi-level MDO formulation allows to find feasible solutions with respect to trajectory and sizing constraints that are difficult to satisfy with using classical single-level approaches. The proposed strategy has been compared on a reusable launch vehicle design problem to a multi-level formulation using a multi-objective Particle Swarm Optimization algorithm (OMOPSO). The proposed approach converges faster to a better Pareto front than OMOPSO. This multi-mission performance study should be part of a broader trade-off between all the Pareto front candidates to include the refurbishment operations, the maintenance aspects and the associated cost of reusability in order to identify the most interesting candidate. In future works, a more complex multi-mission launch vehicle design problem involving higher fidelity codes for the different disciplines will be investigated using advanced Deep Gaussian Processes [START_REF] Damianou | Deep gaussian processes[END_REF][START_REF] Hebbal | Multi-objective optimization using Deep Gaussian Processes: Application to Aerospace Vehicle Design[END_REF].

  Brevault * , Mathieu Balesdent † ONERA, DTIS, Université Paris Saclay, F-91123 Palaiseau Cedex, France Ali Hebbal ‡ ONERA, DTIS, Université Paris Saclay, F-91123 Palaiseau Cedex, France, Université de Lille, CNRS/CRIStAL, Inria Lille Reusability of the first stage of launch vehicles may offer new perspectives to lower the cost of payload injection into orbit if sufficient reliability and efficient refurbishment can be achieved. One possible option that may be explored is to design the vehicle first stage for both reusable and expendable uses, in order to increase the flexibility and adaptability to different target missions. This paper proposes a multi-level Multidisciplinary Design Optimization (MDO) approach to design aerospace vehicles addressing multi-mission problems. The proposed approach is focused on the design of a family of launchers for different missions sharing commonalities using multi-objective MDO to account for the computational cost associated with the discipline simulations. The multi-mission problem addressed considers two missions: a reusable configuration for a Sun Synchronous Orbit with a medium payload range and recovery of the first stage using a gliding-back strategy; and an expendable configuration for a medium payload injected into a Geostationary Transfer Orbit. A dedicated MDO formulation introducing couplings between the missions is proposed in order to efficiently solve such coupled problem while limiting the number of calls to the exact MultiDisciplinary Analysis thanks to the use of Gaussian Processes and multi-objective Efficient Global Optimization. response vector g i (•) = i th constraint function Z N = design of experiment of size N k Θ (•, •) = covariance function parameterized by Θ N = normal distribution p(•) = probability distribution function E HV I = Expected HyperVolume Improvement GLOW = Gross Lift-Off Weight M pr op = propellant mass stage of launch vehicles may offer new perspectives to lower the cost of payload injection into orbit if sufficient reliability and low refurbishment costs can be achieved. In function of various hypotheses on addressable institutional and commercial markets, different strategies for the reusability of the launch vehicle first

  a finite hypervolume of the objective space where all the possible solutions lie, with y L = [min f 1 (z), . . . , min f n (z)] the ideal objective candidate and y U a chosen upper point (nadir point). The dominated hypervolume H Y N of the DoE responses Y N is defined as:

Figure 2 zFig. 2

 22 Fig. 2 Hypervolume improvement (hatched area) in a two objective case

  The characteristics of the SSO mission are: Target orbit SSO, 800km Launch site CSG, Kourou, French Guiana Payload mass 6t Vehicle type Two-Stage-To-Orbit (TSTO) reusable vehicle The characteristics of the GTO mission are: Target orbit GTO, 250km-36000km Launch site CSG, Kourou, French Guiana Payload mass To be determined Vehicle type Two-Stage-To-Orbit (TSTO) expendable vehicle

Fig. 4 Fig. 6

 46 Fig. 4 Illustration of the reusable launch vehicle architecture

Fig. 7

 7 Fig. 7 Steps of the gliding-back mode recovery

Fig. 8 Fig. 9

 89 Fig. 8 Single-level coupled MDO formulation for multi-mission launch vehicle design

  ,-turbopumpk,valvesk,pipingk, thrust,chamberk,igniterk,injectork,pre-combustorN, ,from,Propulsion Inputs,coming,from,Trajectory Inputs,coming,from,Sizing

Fig. 14

 14 Fig. 14 Comparison of the Neural Network prediction (line) and the exact MISSILE function evaluations (crosses)

V

  norm of the velocity vector (m.s -1 ) γ flight path angle (rad) φ latitude (rad) λ longitude (rad) ψ flight path heading (rad) µ bank angle (rad) θ pitch angle (rad) ω E angular velocity of the Earth (rad/s)

Fig. 15

 15 Fig. 15 Earth-centered, Earth-fixed reference frame

Figure 16

 16 Figure 16 represents a box-plot of 20 repetitions of expendable subsystem problem solving with different initializations of SQP for a given system-level design variable vector. The range of convergence for the objective function varies between 2.75t and 4.25t for the maximal payload mass which is clearly non robust to the initialization.

Fig. 16

 16 Fig. 16 Boxplot of the optimal objective function for the expendable subproblem optimization using a SQP algorithm -20 repetitions with random start

Figures 17 ,

 17 Figures 17,[START_REF] Rasmussen | Gaussian processes in machine learning[END_REF][START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with Gaussian process model[END_REF] represent a run of CMA-ES algorithm for the expendable optimization subproblem with the convergence plots of the objective function (Figure17), the design variables (Figure18) and the associated parallel plot (Figure19). The convergence of the design variable values and the penalized objective function is reached by CMA-ES

Fig. 19

 19 Fig. 19 Parallel plot (CMA-ES) -expendable subproblem

Fig. 22 InitialFigure 30

 2230 Fig. 22 Initial DoE samples, initial and final Pareto fronts -EGO multi-level formulation

Fig. 23 Fig. 24

 2324 Fig. 23 Initial and final DoE samples, final Pareto front -EGO multi-level formulation

Fig. 25 Fig. 26

 2526 Fig. 25 Comparison final Pareto fronts -EGO and OMOPSO, multi-level formulation

Fig. 27 Fig. 28 Fig. 29

 272829 Fig. 27 Evolution of Pareto front for EGO at iteration 50, 100, 150 and 200. Triangles are current DoE, line current Pareto front, circles final Pareto front

Figures 31 ,

 31 Figures 31,[START_REF] Hebbal | Multi-objective optimization using Deep Gaussian Processes: Application to Aerospace Vehicle Design[END_REF] and 33 present the ascent and re-entry to landing site trajectories for the expendable and reusable configurations corresponding to the minimal GLOW and maximal payload on the obtained Pareto frontier. These two solutions are the extreme points of the Pareto front. For these solutions, the altitude, the velocity and the mass are displayed as a function of time. The coasting phase (corresponding to the ballistic flight) is not represented in these figures. An illustration of the payload injection on GTO is represented in Figure29for the expendable maximal payload

Fig. 30 Fig. 32

 3032 Fig. 30 Parallel plots with respect to system-level design variables -multi-objective EGO (in upper graph, colors correspond to maximal payload mass objective; in lower graph, colors correspond to minimal GLOW objective)

  w.r.t. z = M pr op1 , M prop2 , S wing , λ wing , β M pr op1 , c r euse , c ex p , N X d , N Z d , Pdyn d

	s.t. | h apogee e x p (z) -h t apogee e x p |	≤ h a	(12)
	| h perigee e x p (z) -h t perigee e x p |	≤ h p	(13)
	NX max e x p (z)	≤ NX t max e x p	(14)
	Pdyn max e x p (z)	≤ Pdyn t max e x p	(15)
	Flux max e x p (z)	≤ Flux t max e x p	(16)
	α max e x p (z)	≤ α t max e x p	(17)
	| h apogee r eu s e (z) -h t apogee r eu s e |	≤ h a	(18)
	| h perigee r eu s e (z) -h t perigee r eu s e |	≤ h p	(19)
	NX max r eu s e (z)	≤ NX t max r eu s e	(20)
	NZ max r eu s e (z)	≤ NZ t max r eu s e	(21)
	Pdyn max r eu s e (z)	≤ Pdyn t max r eu s e	(22)
	Flux max r eu s e (z)	≤ Flux t max r eu s e	(23)
	d sitemax r eu s e (z)	≤ d t sitemax r eu s e	(24)
	α max r eu s e (z)	≤ α t max r eu s e

  )h t apogee e x p | ≤ h a | h perigee e x p (z, z * )h t perigee e x p | ≤ h p NX max e x p (z, z * ) ≤ NX t )h t apogee e x p | ≤ h a | h apogee r eu s e (z, z * )h t apogee r eu s e | ≤ h a | h perigee r eu s e (z, z * )h t perigee r eu s e | ≤ h p NX max r eu s e (z, z * ) ≤ NX t Flux max r eu s e (z, z

		Multi]objectiveLoptimizer
		zN=[Mprop1ZLMprop2ZLNXdZNZdZPdynd]
	ReusableLconfLoptimizer	ExpendableLconfLoptimizer
	z=[δMprop1ZδMprop2ZSwingZ LLLLLλwingZLβMprop1ZLθreuse]L		z=[δMprop1ZδMprop2ZLθexp]L
	Propulsion
			Propulsion
	GeometryL[
	sizing		GeometryL[
			sizing
	Aerodynamics
			Aerodynamics
	Ascent
	trajectory		Ascent
		Return		trajectory
	ReusableLconfiguration	trajectory	ExpendableLconfiguration
		GLOWreuse
		MpayloadLexp
		Constraints
	s.t.	J cons e x p (z * ) ≤ ex p		(28)
		J cons r eu s e (z * ) ≤ reuse		(29)
	z * min ≤ z max e x p Pdyn max e x p (z, z * ) ≤ Pdyn t max e x p Flux max e x p (z, z * ) ≤ Flux t max e x p α max e x p (z) ≤ α t max e x p z min ≤ z ≤ z max	min	GLOW reuse (z, z max r eu s e NZ max r eu s e (z, z * ) ≤ NZ t max r eu s e Pdyn max r eu s e (z, z * ) ≤ Pdyn t max r eu s e

w.r.t. z * = M prop1 , M prop2 , N X d , N Z d , Pdyn d * ≤ z * max (30) Subsystem-level: min -M payload ex p (z, z * ) w.r.t. z = δ M pr o p1 , δ M pr o p2 , c ex p s.t. | h apogee e x p (z, z * * ) w.r.t. z = δ M pr o p1 , δ M pr o p2 , S wing , λ wing , β M prop1 , c reuse s.t. | h apogee e x p (z, z * * ) ≤ Flux t max r eu s e d sitemax r eu s e (z, z * ) ≤ d t sitemax r eu s e α max r eu s e (z) ≤ α t max r eu s e z min ≤ z ≤ z max

Acknowledgments

This work is funded by the ONERA (Office National d'Etudes et de Recherches Aérospatiales -The French Aerospace Lab) project HERACLES (Hypersonic Efficient and Reusable Aerospace Concepts for Launcher Evolution Strategies, 2017-2020). The PhD thesis of A. Hebbal is funded by ONERA and the University of Lille. Support of CNES is gratefully acknowledged.

B. Result analysis and discussions

To illustrate the need for a multi-level MDO formulation, the single-level decoupled approach presented in Section IV.B.1 using an OMOPSO algorithm [START_REF] Sierra | Improving PSO-based multi-objective optimization using crowding, mutation and epsilon-dominance[END_REF].

A population of 20 individuals are considered and OMOPSO runs for 500 iterations corresponding to 10000 evaluations of the subsystem coupled models (for the reusable and expendable configurations). Among all the evaluations, only around 10% of solutions are feasible with respect to the system-level constraints (Figure 20). Moreover, the Pareto front found by OMOPSO with the single-level formulation is far from the Pareto front obtained with the multi-level formulation using EGO (Figure 21).

The proposed multi-level MDO approach starts with an initial DoE size of 20 generated by Latin Hypercube Sampling (in dimension 5) and 200 iterations are run before stopping the algorithm. The proposed multi-level MDO approach is compared to a reference OMOPSO optimization algorithm [START_REF] Sierra | Improving PSO-based multi-objective optimization using crowding, mutation and epsilon-dominance[END_REF] in terms of quality of the obtained Pareto front (Figures 24 and26). The OMOPSO has a population of 20 individuals and runs for 11 iterations before stopping (therefore OMOPSO and EGO have the same number of exact MDA evaluations). To evaluate the robustness to the initial DoE (or initial population of OMOPSO), ten repetitions of the multi-level MDO approach and OMOPSO optimization are carried out. On a standard desktop machine (using Intel Xeon CPU E5-2650 with 16 cores for parallelization), one


Le Flécher
email: sentenac@laas.frajsd@cin.ufpe.br

Visual predictive control of robotic arms with overlapping workspace

Keywords: Agricultural robotics, multi-arms, Non-linear Model Predictive Control, Image Based Visual Servoing

This paper deals with multi-arms fruits picking in orchards. More specifically, the goal is to control the arms to approach the fruits position. To achieve this task a VPC strategy has been designed to take into account the dynamic of the environment as well as the various constraints inherent to the mechanical system, visual servoing manipulation and shared workspace. Our solution has been evaluated in simulation using on PR2 arms model. Different models of visual features prediction have been tested and the entire VPC strategy has been run on various cases. The obtained results show the interest and the efficiency of this strategy to perform a fruit picking task.

Introduction

To meet the increasing food demands from a growing world population, agriculture will need to double its production by 2050. Alternatively, it will need to reduce its environmental impact to fight against climate change and avoid to harm soils. Among the solutions presented in [START_REF] Foley | Solutions for a cultivated planet[END_REF][START_REF] Grift | A review of automation and robotics for the bio-industry[END_REF], robotics has been designated as one of the most promising strategies.

Agriculture robotics may operate in various environments such as open fields, orchards, green houses, ... In this work, one focuses on robots performing weeding, spraying, trimming, or harvesting in orchards. The autonomous achievement of any of these tasks requires the same capacities for the robotic system: (i) autonomously navigate in an orchard; (ii) detect and localize a set of natural landmarks of interest; (iii) approach the landmarks with some robotic manipulators; (iv) manipulate/interact with the landmarks of interest without damaging them. After addressing the autonomous navigation problem in [START_REF] Durand-Petiteville | Design of a sensor-based controller performing u-turn to navigate in orchards[END_REF][START_REF] Flecher | Implementation on a harvesting robot of a sensor-based controller performing a u-turn[END_REF], the focus is now on the landmark approach problem when considering a couple of robotic manipulators embedded on a mobile platform. This topic has already been addressed during the design of systems harvesting cucumbers (Van Henten et al., a https://orcid.org/0000-0002-2683-8675 2003) or tomatoes [START_REF] Zhao | Dualarm Robot Design and Testing for Harvesting Tomato in Greenhouse[END_REF]. Despite the promising results, the harvesting efficiency needs to be increased to meet the viability demand from farmers (Zhao et al., 2016a). To do so, one of the most adopted approach consists of embedding several robotics arms on the mobile platform, as it is done in (Zhao et al., 2016a) for a melon harvester. However, it is pointed in [START_REF] Vougioukas | A Study of Fruit Reachability in Orchard Trees by Linear-Only Motion[END_REF]) that fruits are not uniformly distributed on commercial orchards trees. Indeed, to maximize the fruits accessibility, growers use specific trees with high fruits density, thus facilitating the manual harvesting. For this reason, during the design of a harvesting system, it seems then interesting to consider the use of several arms whose workspace overlap. Thus, even if the arms do not collaborate to manipulate one single object, they share a common workspace. It is then possible to use several arms to pick fruits in the same high density area, which decreases the overall harvesting time.

Regarding the control of the robotics arms, it seems interesting to focus on reactive controllers. Thus, when considering the harvesting problem, a tree is a highly dynamic environment because of the flexibility of the branches and the weight of the fruits: each time a fruit is harvested, the local distribution of fruits and branches is modified. It is then proposed to use an image based visual servoing (IBVS) scheme for its reactive capabilities [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF]. This kind of controller will track the fruit of interest and use its coordinates in the image to control the arm. However, IBVS schemes working in the image space, they do not provide any guarantee regarding the trajectory of the camera in the workspace. This represents a major issue when considering arms sharing a common workspace. A well adopted solution to take into account constraints on position, velocity, or workspace is the model predictive control scheme. Indeed, when minimizing the cost function, the solver can deal with numerous constraints [START_REF] Diehl | Fast motions in biomechanics and robotics: optimization and feedback control[END_REF]. Thus, to obtain a reactive controller dealing with constraints, it is proposed to design an IBVS-based NMPC (nonlinear model predictive control) controller, also named Visual Predictive Control (VPC). For this particular case, the cost function is defined as the error between the image coordinates of the visual features representing the landmark of interest at the current pose and the desired one.

Similar control schemes have been already proposed. For example, in [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF] several prediction models are used to simulate the control of a 6 DOF camera subject to visibility constraints. In [START_REF] Sauvee | Image Based Visual Servoing through Nonlinear Model Predictive Control[END_REF], the designed controller deals with torques, joints and image constraints and is evaluated in simulations. In [START_REF] Hajiloo | Robust On-Line Model Predictive Control for a Constrained Image Based Visual Servoing[END_REF] an IBVS-based MPC controller using a tensor product model is developed. It also takes into account actuators limitations and visibility constraints.

In this paper, it is proposed to extend the VPC case to the control of two arms sharing a common workspace. Thus, in addition to drive the two endeffectors to their desired pose using the image space, the controllers have to guarantee the non-collision between both arms.

A non-linear predictive control strategy to perform a reactive control based on IBVS method is simulated on a robot equipped of multiple arms with eye-inhand cameras. Each arm needs to achieve separate tasks sharing the same workspace. These tasks are subject to constraints such as the restricted field of view, the kinematic (joints limits in position and velocity) but also the self collision constraints.

The paper is organized as follows. In section II the robot model to estimate the location of visual features is defined. In section III, the control strategy is explained, as well as the cost function described and the diverse constraints designed. Section IV is dedicated to the simulation results obtained while section V summarizes and develops the expectations for the future.

Modeling

In this section, one first focuses on the modeling of the system, i.e., the robotic arms, the cameras and the targets. Then, two models to predict the visual feature coordinates of a point are presented. One ends with a discussion regarding the viability of these models for the visual predictive control problem.

System modeling

The robotic system considered in this work is composed of n a identical arms, each of them containing m q revolute joints. A camera is attached to the end effector of each arm. To model, the system, let us first define a global frame F b = (O b , x b , y b , z b ) attached to the system base. Moreover, n a frames

F c i = (O c i , x c i , y c i , z c i ), with i ∈ [1, ..., n a ]
, are used to represent the pose of the camera attached to the i th arm.

One denotes q i j and qi j , with j ∈ [1, ..., m q ], respectively the angular position and velocity of the j th joint of the i th arm. Thus, for the i th arm, one obtains the configuration vector Q i = [q i1 , q i2 , ..., q im q ] T and the control input vector Qi = [ qi1 , qi2 , ..., qim q ] T . Finally, we define for the whole system

Q = [Q T 1 , ..., Q T n a ] T and Q = [ QT 1 , ..., QT n a ]
T . The cameras mounted on the end effectors are modeled using the pinhole model. Thus, the projection matrices H I i /c i mapping the 3D coordinates (x, y, z) of a point in the camera frame F c i to its 2D projection (X,Y ) on the image plan F I i is defined as follows:

   X Y z 1    =    f /z 0 0 0 0 f /z 0 0 0 0 1 0 0 0 0 1       x y z 1    (1)
In this work, the arms are controlled to make each camera reach a pose defined by the coordinates of point visual features in the image space. To do so, one uses n a landmarks made of four points. When considering the i th camera/target couple, the coordinates of the projection of each point on the image is denoted

S il = [X il ,Y il ] T , with i ∈ [1, ..., n a ] and l ∈ [1, 2, 3, 4].
Thus, the visual feature vector for each camera is defined as S i = [S i1 , S i2 , S i3 , S i4 ] T , whereas the system one is given by S = [S T 1 , ..., S T n a ] T . In the same way, a vector for the visual feature desired coordinates of each camera is defined as

S * i = [S * i1 , S * i2 , S * i3 , S * i4 ] T .

Prediction models

A visual predictive control scheme requires to be able to predict the values of the visual features for a given sequence of control inputs. In [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF], the authors present two ways to perform the prediction step for a flying camera: the local and global approaches. This section is devoted to the presentation of both methods for the case of a 6 degrees of freedom camera embedded on a robotic arm.

The global model

The global approach is based on the robot geometric model and the camera projection matrix. The idea consists of computing the 3D coordinates of the features in the next camera frame, then to project them onto the camera image plan. To do so, let us consider a discrete time system with a sampling period T s , where t k+1 = t k + T s , and the notation .(k) = .(t k ).

The configuration Q i (k + 1) of the i th arm, with i ∈ [1, ..., n a ], at the next iteration obtained after applying the control inputs Qi (k) is given by:

Q i (k + 1) = Q i (k) + Qi (k)T s (2)
Knowing the arm configuration at instants t k and t k+1 , it is then possible to compute the homogeneous transformation matrices between the camera frame F c i and the base one

F b . ōF b (.) = H b/c i (.) ōF c i (.) (3)
where ō are the homogeneous coordinates of a point feature expressed in the base or camera frame. Finally, coupling (3) with (1), one obtains a prediction of the coordinates of a point feature in the image space.

S il (k + 1) =

H I i /c i (k + 1)H -1 b/c i (k + 1)H b/c i (k)H -1 I i /c i (k)S il (k) (4)

The local model

The local approach relies on the differential equation mapping the derivative of the visual feature coordinates Ṡil to the robot control input vector Qi . To obtain this equation, one uses the robot's Jacobian J i , mapping the camera kinematic screw T c i and Qi (5), and the interaction matrix L il , mapping Ṡil to T c i (6).

T c i = J i Qi (5) Ṡil = L il T c i (6
) where

L il = -1 z l 0 X l z l X l Y l -(1 + X 2 l ) Y l 0 -1 z l Y l z l 1 +Y 2 l -X l Y l X l (7)
By combining equations ( 5) and ( 6), one obtains the following local model for each arm/camera couple:

Ṡil = L il J i Qi (8)

Usability of the models

When considering using one of the models on an experimental system, several aspects have to be taken into account. First, both models require a measure/estimate of the visual feature depth z (see equations ( 1) and ( 7)). However, when using a monocular camera, this value can only be retrieved by estimation.

Another widely adopted option, consists of using a constant value, such as the depth value at the initial or desired pose [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF].

The second issue is related to the integration of the local model to obtain the prediction equation. A first option consists of using an advanced numeric scheme to accurately compute the integration, whereas a second one relies on a linear approximation of the model, such as:

S il (k + 1) = S il (k) + L il J i Qi (k)T s (9)
Thus, depending on the user choices regarding these issues, the prediction feature coordinates will be more or less accurate. In section 4, different configurations will be evaluated for our specific application in order to select the most appropriate and realistic solution.

Controller design

In this section the control of the arms motion to reach their own targets and compounding necessary constraints are presented. To achieve these goals, NMPC control has been implemented to find the optimal inputs vector which minimizes a cost function subjects to constraints. To do so, we first present the cost function to be minimized in order to reach the desired targets. Then, the expression of each constraint will be defined.

Visual Predictive Control

The proposed VPC scheme couples NMPC (Nonlinear Model Predictive Control) with IBVS (Image-Based Visual Servoing). Thus, similarly to NMPC, it consists of computing an optimal control sequence Q(.) that minimizes a cost function J N P over a prediction horizon of N P steps while taking into account a set of user-defined constraints C( Q(.)). Moreover, similarly to IBVS, the task to achieve is defined as an error in the image space. Thus, the cost function to minimize is the sum of the quadratic error between the visual feature coordinates Ŝ(.) predicted over the horizon N p and the desired ones S * . It then possible to write the VPC problem as:

Q(.) = min Q(.) J N P (S(k), Q(.)) (10)
with

J N P (S(k), Q(.)) = k+N p -1 ∑ p=k+1 [ Ŝ(p) -S * ] T [ Ŝ(p) -S * ] (11) subject to Ŝ(.) = (4) or (8) (12a) Ŝ(k) = S(k) (12b) C( Q(.)) ≤ 0 (12c)
As mentioned in the previous section, the visual feature coordinates can be predicted either using the global or local models (equation ( 12a)). To compute the prediction, one uses the last measured visual features (equation ( 12b)). Finally, equation (12c) defines the set of constraints to be taken into account when minimizing the cost function ( 11). Thus, solving equation ( 10) leads to the optimal sequence of control inputs Q(.). In this work, and as it is usually done, only the first element Q(1) of the sequence is applied to the system. At the next iteration, the minimization problem is restarted, and a new sequence of optimal control inputs is computed. This loop is repeated until the task is achieved, e.i., the cost function ( 11) is vanished.

Constraints

As seen in the previous section, it is possible to add a set of constraints in the minimization problem. In this section, one presents the constraints related to the geometry of the robot, the velocity of the joints, and the field of view of the camera. One finishes with defining a set of constraints that guarantee the non-collision of the arms despite sharing a common workspace.

System constraints

The joint velocity lower and upper bounds are respectively denoted Qmin and Qmax . They are two N Q long vector, with N Q = n a m q N p . They allow to define the lower and upper limit of each joint of each arm over the prediction horizon. Consequently, the constraint on Q(.) is described in the domain κ = [ Qmin , Qmax ], and the velocity constraint is then written as:

Q -Qmax Qmin -Q ≤ 0 (13)
Mechanical limits follow the same principle. The lower and upper bounds of the angular values are denoted Q min and Q max . They are two N Q long vector, with N Q = n a m q N p . These constraints are then given by:

Q -Q max Q min -Q ≤ 0 (14)
Finally, the visual constraints are defined by the image dimension. Ones denotes the bounds by S min and S max . They are two N S long vector, with N S = 8n a N p . Thus, all the predicted visual features coordinates need to be kept into this space, and the constraints can then be written as:

S -S max S min -S ≤ 0 (15)

Overlapping workspace constraint

Due to the overlapping workspace, it is mandatory to guarantee the non-collision between the arms. To do so, one computes the shortest distance d i, j|i , j between the j th joint of the i th arm and the j th joint of the i th arm. The closest distance is evaluated as the norm of the segment perpendicular to both assessed links as it is illustrated in figure 1. Finally, in order to avoid collision, one defines the following constraint:

d min -d ≤ 0 ( 16
)
where d = [d 1,1|2,1 , ..., d n a -1,m q |n a ,m q ] and d min is a predefined distance.

Figure 1: Computation of the shortest distance d i, j|i , j between links L i j and L i j .

Terminal constraint

To make the system converge, it is mandatory to guarantee that the predictive control problem is feasible [START_REF] Lars | Nonlinear model predictive control[END_REF]. To check the feasability, one adds a final constraint in the optimization problem. It is defined as the error between the prediction of the visual feature coordinates Ŝ(k + N p -1) obtained at the end of the prediction horizon, and the desired ones S * . If the solver cannot compute a sequence of control inputs Q that respects this constraint, then the problem is not feasible, and there is no guarantee regarding the system convergence.

Ŝ(k + N p -1) -S * ≤ 0 (17)
All the constraints are now defined and can be integrated into the non-linear function C( Q(.)):

C( Q) ≤ 0 ⇒            Q -Qmax Qmin - Q Q -Q max Q min -Q S -S max S min -S Dmin -D Ŝ(k + N p -1) -S *            ≤ 0 (18)
Thus, the optimal input must belong to the domain defined by C( Q(.)).

Simulation

In this section, one presents the results obtained by simulating the control of two arms with MATLAB software. A 7 DOF model has been selected, based on the arms embedded on the PR2 robot. In these simulations the task consists of positioning the cameras attached to each end-effector with respect to two given targets. Each desired position is defined in the image space by S i = [x 1 , y 1 , x 2 , y 2 , x 3 , y 3 , x 4 , y 4 ], which are the coordinates of four points representing the corresponding target. Moreover, the prediction horizon is set up at N p = 5. Regarding the constraints, the angular velocity limits for the first prediction step are given by κ = [-0.2 0.2] for each joint, whereas it is extended to κ = [-5 5] for the rest of the prediction step. This allows to decrease the computation while guaranteeing a feasible control input for the first step. The size of the image is S min = [-4 -4] and S max = [4 4]. The angular constraints on each joint are setup as described in the PR2 documentation [START_REF] Garage | Pr2 user manual[END_REF]. Finally the minimum collision avoidance distance is fixed to D min = 0.1. The initial coordinate of the visual features are computed based on the initial configuration of the robot given by q init = [0 π 2 0 0 0 0 0 0 π 2 0 0 0 0 0] T for every simulation. Finally, one fixes T s = 50ms.

Comparison between different prediction methods

The first set of simulation assesses the accuracy of different predictive models presented in section 2.2.

To do so, one drives the arms from q init to q * = [-π 8 π 2 0 0 0 0 0 0 π 2 0 0 0 0 0] T . During this motion, the coordinates of the visual features predicted over the prediction horizon, are compared with a ground truth. In this case, it consists of the prediction obtained with the global model computed with the real depth of each visual feature.

Figure 2 shows the evolution over the simulation of the quadratic error between the predicted coordinates of one visual feature and the ground truth. Six models are compared by combining the global model (solid lines), local Euler approximation (dashed lines) and local Runge Kutta approximation (dotted lines) with the depth values real z (black) and desired z * (blue). It can be seen that the global models and the local Runge Kutta are significantly more accurate than the local Euler one. Indeed, using this latter prediction scheme introduces a large error, and might not allow the system to converge toward the desired configuration. In a second figure 3, the previous errors were summed over the simulation to offer a better view of model performances. When considering the use of the real depth, the global model is the most accurate. However, this result is obtained in a simulation, and its use in a real context would require an estimator of the depth. Thus, one now considers one more realistic case where the depth is approximated with the constant value z * . In this case, the local Runge Kutta is the most accurate model. Thus, for the rest of the simulations the local Runge Kutta model with z * is used to evaluate the performances in a more realistic scenario. π 8 0 0 0] T from the initial one q init while dealing with the previously presented constraints. This visual predictive control was performed using the local model with Runge Kutta approximation and the desired depth z * . The trajectories in the images of visual features are displayed in figure 4 and show the achievement of the task for both cameras. Indeed, the proposed controller was able to drive both arms to make the cameras converge towards the desired features from the initial ones (respectively represented by the red and green crosses in figures 4(a) and 4(b)). To achieve the task in the image space, the arms had to reach the poses presented in figure 5. During the realization of this task, the image limits were never reached. Similarly there was no risk of collision during the servoing. Moreover, it can be seen in figures 6(a) and 6(b) that the constraints on the joint configurations were taken into account. Similarly, in figures 6(c) and 6(d), one can see that the constraints on the joint velocities were taken into account to achieve the task. Indeed, the applied commands stay within the given range κ = [-0.2 0.2] all along the simulation.

Collision avoidance

In this last simulation, the system has to reach the desired configuration q * = [ π 16 3π 8 0 π 4 0 π 8 0π 16 π 2 π 4 π 8 0 0 0] T . This set is run using the Runge Kutta local model with the desired depth values. One more time the desired visual features are reached from the initial poses (figure 7). However, as it can be seen in figure 8, the arms had to cross each other to reach the desired poses. Thanks to the shared environment constraint included in the control problem, the collision was avoided. Indeed, figure 9 shows that the distance between the different links was never smaller than the fixed value D min = 0.1. In addition, the constraints on the joints angular values are respected. In figure 10, it can be seen that the right arm joints never reach their angular limits represented by red dashed lines (the joints represented in figures 10(e) and 10(g) do not have limits on PR2 robot). Similarly, the control inputs of the right arm present saturations due to the angular velocity constraints, represented by red dashed lines in figure 11, but never goes beyond.

Conclusion

In this paper, we have addressed the problem of fruit picking in orchards performed by a multi-arms robot. The goal is to bring the different end-effectors close to their designated fruits while taking into account the dynamic environment as well as the constraints inherent to the system mechanics, the visual manipulation and the shared workspace. The proposed solution relies on a VPC strategy to fulfill these requirements. First, several prediction models were presented. Next, the visual servoing task and the constraints were described including them in a NMPC strategy. Finally, simulations based on PR2 manip-(a) q 11 joint angular values (b) q 12 joint angular values (c) q 13 joint angular values (d) q 14 joint angular values (e) q 15 joint angular values (f) q 16 joint angular values (g) q 17 joint angular values Figure 10: Configuration of the right arm over the simulation -Blue solid: Angular joint value q 1 j -Red dashed: Maximum angles ulators were performed to assess the performance of each predictive model and to validate the completion of the servoing. The obtained results show the efficiency of this approach to perform various tasks in a shared workspace. However, multiple steps still need to be accomplished. First, this method must be integrated on a real robot. This step will require to optimize the computational time to run it in real time. In addition, a constraint avoiding the target occlusions in the image must be integrated. Finally, this manipulation controller will be coupled with our previous works about the autonomous navigation in orchards to perform a complete fruits picking task.

Figure 2 :

 2 Figure 2: Quadratic errors of prediction over the simulation -Solid: Global model -Dashed: Local Euler -Dotted: Local Runge Kutta -Black: Real z -Blue: Desired z

Figure 4 :

 4 Figure 4: Visual features location using local Runge-Kutta model and real z -Blue dotted: Trajectories -Red crosses: Final locations -Green crosses: Initial locations

Figure 5 :

 5 Figure 5: Final poses of the arms -Red dots represent the two targets -Green crosses are the initial poses of the end effectors

Figure 7 :

 7 Figure 7: Visual features location using local Runge-Kutta model and real z -Blue dotted: Trajectories -Red crosses: Final locations -Green crosses: Initial locations

  Figure 11: Right arm joints velocities -Blue solid: Joint velocity value -Red dashed: Maximum velocities