Jean-Antoine Désidéri

Julien Wintz

Nathalie Bartoli

Christophe David

Sébastien Defoort

Combining Pareto Optimality with Nash Games in Multi-Objective Prioritized Optimization of an Aircraft Flight Performance

Keywords: differentiable multi-objective optimization, criterion prioritization, Nash game optimisation multiobjectif différentiable, hiérachisation des critères, jeu de Nash

This report reflects a cooperation between the Inria Project Team Acumes and the Information Processing and Systems Department (DTIS) of Onera Toulouse. We present a case study of prioritized multi-objective optimization of the flight performance of an Airbus-A320-type aircraft in terms of take-off fuel mass (f 1 ), operational empty weight (f 2 ) and ascent-to-cruisealtitude duration (f 3 ). The optimization is subject to functional constraints: interval bound on static margin and upper bound on wing span. Six design variables are considered defining geometry, mean aerodynamic chord and take-off thrust potential. Designs are evaluated by means of the FAST-OAD open-source software developed by ONERA and ISAE-SUPAERO, and the prioritized optimization is conducted by the Inria Nash-MGDA software. From a point of Pareto-optimality of the sole pair (f 1 , f 2 ), considered preponderant, the climb duration f 3 is reduced drastically by a one-shot procedure establishing a continuum of Nash equilibria. The application of this procedure is far less computationally demanding than is the numerical establishment of the Pareto-front associated with the three cost functions together, while the Pareto-optimality of the preponderant cost functions (f 1 , f 2 ) is only marginally degraded.

Introduction

In multi-objective optimization, the notion of Pareto-optimality is central [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. In order to recall this concept, we consider the following multi-objective continuous optimization problem min x∈Ω ad {f 1 (x), f 2 (x), . . . , f M (x)} subject to: c k (x) = 0 (k = 1, . . . , K) [START_REF] Allmendinger | What if we increase the number of objectives? theoretical and empirical implications for many-objectve combinatorial optimization[END_REF] in which Ω ad is the admissible domain, usually a subset of R n , possibly bound by inequality or interval constraints on the components of the optimization variable x, and {f j (x)} (j = 1, . . . , M ) and {c k } (k = 1, . . . , K) are the cost and (functional) constraint functions. We assume that the possible bound or interval constraints on the optimization variables are accounted for by a special parameterization as in Section 4. Hence we refer here to nonlinear constraints of other types, that we designate as "functional constraints", in the sense that they are calculated together with the cost functions by some nonlinear procedures, subject or not to the solution of a system of partial-differential equations (PDEs). Note that bound or interval (functional) constraints can always be formulated as equality constraints by the introduction of slack variables, thus increasing the dimension n, and in fact, this technique is used and made explicit in the present application (see Section 4). Hence the above formulation is fairly general.

For the above multi-objective problem, a design-point x ∈ Ω ad is said to "dominate in efficiency" the design-point y ∈ Ω ad , if and ony if (iff) its performance is at least equal to the performance of y with respect to (w.r.t.) all cost functions, ∀j = 1, . . . , M :

f j (x) ≤ f j (y) (2)
and strictly superior w.r.t one cost function at least (strict inequality for one j at least). Clearly, this notion defines a relationship of partial order. The "Pareto set" is the subset of the admissible domain made of the non-dominated design-points, and the "Pareto front" its image in function space. Strictly speaking, the Pareto set and front constitute the mathematical solution to the problem [START_REF] Allmendinger | What if we increase the number of objectives? theoretical and empirical implications for many-objectve combinatorial optimization[END_REF], and determining these sets is considered by many as the ultimate designer's graal. However this viewpoint may be nuanced in the context of complex multi-disciplinary optimization often named concurrent engineering. For purpose of discussion, let us consider the somewhat idealized, nevertheless instructive situation in which a given population of design-points is evaluated several times w.r.t. an ordered list of cost functions {f 1 , f 2 , . . . } assumed to be continuous and convex over a closed and convex domain. Given the integer k (k ≥ 2), if S k denotes the Pareto set w.r.t. the first k functions, {f 1 , f 2 , . . . , f k } with the understanding that S 1 = arg min f 1 . Then we generally have the inclusions

S 1 ⊆ S 2 ⊆ • • • ⊆ S k . . . (3)
(see Appendix A for a justification). Hence each new evaluation with a larger number of cost functions can only have the effect of enlarging the Pareto set. One could imagine some absurd situation in which at the end, no design point from an original finite population could be excluded by the Pareto-optimality criterion, and the entire given population would constitute the Pareto set itself 1 . Without going to this unrealistic extreme, we conclude that the Pareto set becomes less and less informative as the number of cost functions increases, and this is counter-productive in the design of a complex engineering system. An alternate conclusion is that if the number of cost functions is increased, in a numerical context, the population of candidate design-points to be tested should be enriched substantially, thus adding another element that increases the computational demand to determine the Pareto set.

Let us return to [START_REF] Allmendinger | What if we increase the number of objectives? theoretical and empirical implications for many-objectve combinatorial optimization[END_REF] in the case of a typical engineering system for which the runtime of the computational procedure that evaluates cost and constraint functions is not negligible. Then, this factor may critically condition the type of potential methods to be used to determine the Pareto 1 This would be the case for the finite population {x i } (i = 1, . . . , N ), w.r.t. the cost functions f j (x) = xx j 2 (j = 1, . . . , N ).

RR n°9490

set. Even in case of a moderate runtime, say of the order of the minute, the approximate determination of the Pareto set may be computationally demanding enough to exclude the practitioner's interactive action in the course of the optimization process. Anticipating Section 7, in the optimization of an aircraft flight-performance considered in this report, the average runtime for one evaluation on a laptop of the sole cost functions, without their gradients, has been measured to be about 30 s with our software, and for the numerical determination of the Pareto front associated with two cost functions by the classical NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] to more than 22 hrs for a population of 50 individuals. Hence developing computationally efficient optimization strategies in such a context remains of great importance, particularly in the perspective of the treatment of systems governed by partial-differential equations (PDEs), most often out of reach today.

Nevertheless, if the Pareto set has been identified numerically, even laboriously, the decision by the designer of which non-dominated design-point to elect is still to be made. In a commonly-used basic decision scheme, the choice is made by evaluating subsequently an additional criterion at the Pareto-optimal design-points, and choosing one of the best points w.r.t. this new criterion. In such a scheme, the decision is made in two steps. The introduction of an auxiliary criterion, in effect, relies on a hierarchical split made between the cost functions for which computing the Pareto front is consented, designated as "primary cost functions" hereafter, and the new cost function (or functions) used to finalize the decision, designated as "secondary cost function(s)". The prioritized multi-objective optimization method promoted here generalizes this split in a more elaborate and systematic way and has been demonstrated to be far more efficient than the basic decision scheme in several applications to engineering.

Our general method includes two optimization phases:

1. First Phase: numerical identification of the Pareto front associated with the constrained optimization of the "primary cost functions", and the election of a starting point x A .

2. Second Phase: one-shot computation of a continuum of Nash equilibria originating from x A permitting to reduce the "secondary cost functions" while best preserving the Paretostationarity of the primary cost functions.

The construction of the algorithm permitting to achieve this is briefly recalled in Section 2 for sake of completeness, while more detailed descriptions can be found in several publications: the standard method, associated software, and first case studies were initially introduced in [START_REF] Désidéri | Direct and adpative approaches to multi-objective optimization[END_REF]7]; the optimization test-case of the Super-Sonic Business Jet (SSBJ) was also presented in [START_REF]Prioritized optimization by Nash games: towards an adaptive multi-objective strategy[END_REF]; [START_REF] Désidéri | Prioritized multi-objective optimization of a sandwich panel[END_REF] relates to the sizing of an aluminum sandwich panel with honeycomb core to maximize its strength under bending loads or a blast and involved four objective functions given by analytical models. Finally, the algorithm was again defined in details in [START_REF]Adaptation by Nash games in gradient-based multi-objective/multi-disciplinary optimization[END_REF] along with convergence proofs, and some summary of accomplished case studies. The latter publication is the most complete reference. The case study considered presently is also one of multi-objective optimization of the flight performance of an aircraft. However here the aircraft is transonic, of A320-type, and the goals of the optimization are defined differently. The evaluation of the cost functions was accomplished by using the open-source software FAST-OAD developed by ONERA and ISAE-SUPAERO (see Section 4).

In the base algorithm proposed in https://mgda.inria.fr, the user provides procedures for computing the cost and constraint functions. Here the algorithm has been split to include a preliminary "Design of Experiment" permitting to execute a number of evaluations, those necessary to construct the initial meta-models, in a batch mode, and limiting the number of calls to the evaluation procedure to one per computed Nash equilibrium point. These modifications are detailed in Section 5. Lastly, the results of the numerical campaign are described in Section 7 where the efficiency gains are analyzed in details, before some concluding remarks are drawn in Section 8.

Base algorithm

From here and throughout, a hierarchy between the cost functions is introduced according to the running index j: Inria 1. Primary cost functions: {f j (x)} (j = 1, . . . , m) 2. Secondary cost functions: {f j (x)} (j = m + 1, . . . , M ) Objective and cost functions are all assumed to be smooth, say C 2 (Ω ad ). A first phase of optimization (omitted here) has led to the election of a design-point x A in the open interior of the admissible domain Ω ad that is Pareto-optimal for the minimization of the sole primary cost functions under the constraints.

The second phase of optimization is aimed to construct a continuum of Nash equilibria {x ε } parameterized by ε (0 ≤ ε ≤ 1) in way such that: firstly, x 0 = x A (compatibility), and secondly, for ε sufficiently small, the Pareto-optimality condition of the primary cost functions remains O(ε 2 ), whereas the secondary cost functions decrease linearly with ε (or faster).

Preliminary calculations

The preliminary calculations have been defined in full in [START_REF] Désidéri | Direct and adpative approaches to multi-objective optimization[END_REF][START_REF]Adaptation by Nash games in gradient-based multi-objective/multi-disciplinary optimization[END_REF] and only summarized here. They are made of three main steps:

• Definition of a convex "Primary Steering Function" f + A (x) • Definition of a "Territory Splitting" • Definition of a "Secondary Steering Function" f B (x).
The cost functions are all assumed to be strictly positive. This may require the preliminary application of an adequate exponential transform, which preserves regularity and sense of variation of the function. This permits us to consider gradients in the logarithmic form, and thus examine the variations of these functions in terms of proportions from their starting values.

The Pareto-stationarity condition satisfied at x = x A is expressed as follows:

m j=1 α j P∇f j f j = 0 (4)
where the supercript on any symbol denotes throughout an evaluation at x = x A , ∇ is the gradient operator w.r.t. x, {α j } (j = 1, . . . , m) are the known coefficients of a convex combination (to be determined), P is the projection operator onto the linear subspace orthogonal to all constraint gradients, {∇c k } (k = 1, . . . , K). Then the primary steering function is defined as follows:

f + A (x) = m j=1 α j f j (x) f j + c 2 x -x A 2 ( 5
)
where c is a nonnegative constant adjusted to enforce the local convexity of f + A (x) and the associated Lagrangian. Then, by virtue of (4), the primary steering function is stationary and convex at x A , and therefore achieves a local minimum.

The territory splitting is based on the orthogonal decomposition of the following n × n reduced Hessian matrix [START_REF] Désidéri | Computational Intelligence in Aerospace Sciences[END_REF][START_REF] Désidéri | Direct and adpative approaches to multi-objective optimization[END_REF]:

H A = PH +, A P = ΩHΩ t (6)
where

H +, A = (∇ 2 f + A ) = H A + cI n and H A = ∇ 2 f A . Matrix H A is positive semi-definite, of rank n -K exactly since matrix H +,
A is strictly-positive definite by convexity fix. The eigenvectors, column-vectors of matrix Ω, are conventionally arranged in a special way: the first K are associated with the directions of the constraint gradients (null space of P), and the remaining ones are ordered by decreasing eigenvalue. In this way, the tail eigenmodes are associated with the smaller sensitivities of the primary steering function. The splitting is associated with the following change of variables:

x = x A + Ω u v := X(u, v) (7)
where u ∈ R n-p , v ∈ R p , and p is an integer free parameter to be chosen such that and 1 ≤ p < n -K. In this setting, the following notations are used:

Ω = (Ω u Ω v ) , Ω u = Ω uu Ω vu , Ω v = Ω uv Ω vv (8)
and the following scaled logarithmic gradients are associated with the secondary cost functions:

g j = S -1 2 ∇ v f j f j = S -1 2 Ω t v ∇f j f j (j = m + 1, . . . , M ) (9)
where the scaling matrix S = (∇ 2 vv f + A ) is the lower p × p diagonal block of matrix H of ( 6), positive-definite by convexity fix. Hence these gradients reflect the sensitivities of the secondary cost functions when v varies, and a scaling by elements related to the primary steering function as well. Following the strategy employed in the Multiple Gradient Descent Algorithm (MGDA) [START_REF] Désidéri | Direct and adpative approaches to multi-objective optimization[END_REF], one identifies the minimum-norm element ω B in the convex hull of the gradients g j :

ω B = M j=m+1 α j g j (10)
and the coefficients {α j } of this convex combination are used to define the secondary steering function:

f B (x) = M j=m+1 α j f j (x) f j (11)
and this completes the preliminaries.

Nash game formulation

One defines a continuation parameter ε (0 ≤ ε ≤ 1) and the auxiliary cost function:

f AB (x) = (1 -ε)f + A (x) + εf B (x) (12)
For every fixed value of ε, two virtual players, P layer A and P layer B, respectively controling the sub-vectors u and v in (7), are involved in a Nash game defined by the following strategies:

-Strategy A: P layer A attempts to minimize f + A X(u, v) by the strategy u, subject to the equality constraint c X(u, v) = 0, and by accounting for P layer B's fixed strategy v, whereas:

-Strategy B: P layer B attempts to minimize f AB (X(u, v)) by the strategy v, subject to no constraints, but by accounting for P layer A's fixed strategy u.

Theoretical results

The following results have been established [START_REF] Désidéri | Direct and adpative approaches to multi-objective optimization[END_REF][START_REF]Adaptation by Nash games in gradient-based multi-objective/multi-disciplinary optimization[END_REF]:

• The Nash equilibrium point x ε exists for all sufficiently small ε, and

x 0 = x A (13)
an essential property referred to as compatibility or consistency.

• If φ A (ε) = f + A (x ε ), φ A (0) = 0 and f + A = 1 + O(ε 2 ) (14)
which reflects the fact that, as ε varies, the Pareto-stationarity of the primary cost functions is degraded by O(ε 2 ) only.

• If φ B (ε) = f B (x ε ), φ B (0) = -σ B , where σ B = ω B 2 ≥ 0.
Hence, the secondary steering function f B (x) decreases linearly as ε increases, and by virtue of properties of the MGDA construction, all the secondary cost functions do as well.

Thus the objectives of our construction are achieved.

These theoretical results are illustrated in Figure 1 in the case of two primary convex cost functions (f 1, f 2) whose associated Pareto front is represented by a convex red curve ("Primary Pareto Front, PPF), and an arbitrary number of secondary cost functions. The continuum of Nash equilibria is represented by the green curve, tangent to the PPF at the chosen starting point x A . In practice, the algorithm is applied to local quadratic metamodels as described in [7]. Some implementation variations have been introduced to gain efficiency. These are defined in Section 5.

m j=1 αjP∇fj/f j = O(ε 2 ) fj/f j = 1 -σBε + O(ε 2 ) (σB ≥ 0) (ε ≥ 0, small enough) x A
3 FAST-OAD simulation software and coupling with the MGDA software

The FAST-OAD software

FAST-OAD (Future Aircraft Sizing Tool -Overall Aircraft Design) is a framework for performing rapid Overall Aircraft Design. The software was developed jointly by ISAE-SUPAERO 2 and ONERA 3 for the preliminary design, analysis and optimization of aircraft. It is a dimensioning code that integrates models describing, with appropriate fidelity levels, the different physics involved in aircraft design: aerodynamics, propulsion, mass, geometry, performance, flight qualities, certification requirements, etc. FAST-OAD makes these different disciplines interact through multidisciplinary analysis and optimization techniques to design and evaluate an aircraft that meets the specifications defined by the user.

The software relies on the OpenMDAO [START_REF] Gray | OpenMDAO: An Open-Source Framework for Multidisciplinary Design, Analysis, and Optimization[END_REF] framework. It allows easy switching between models for a same discipline, and also adding/removing/developing models to match the user's needs. Details can be found in [4] and in the official documentation: https://fast-oad.readthedocs. io/en/v1.3.5/ . See https://github.com/fast-aircraft-design/FAST-OAD for installation of the Python software. This question is revisited in Section 3.3.

The mgda.inria.fr software platform

The MGDA platform (https://mgda.inria.fr) is intended to provide the user with information or software for multi-objective differentiable optimization. The platform is made of four components of increasing complexity:

Basic MGDA tool for finding a search direction. This component allows for computing a descent direction common to an arbitrary set of cost functions whose gradients are provided in situations other than Pareto stationarity.

Directions for solving a constrained problem. This component refers to the research report Inria "Quasi-Riemannian Multiple Gradient Descent Algorithm for constrained multiobjective differential optimization" [START_REF] Désidéri | Quasi-Riemannian multiple gradient descent algorithm for constrained multiobjective differential optimization[END_REF] for solving constrained problems by using the basic MGDA tool.

Tool for solving a prioritized optimization problem, or Nash-MGDA. This component permits to solve a multi-objective optimization problem in which the cost functions are given in two sets:

• a primary set of cost functions subject to constraints for which a Pareto optimal point is provided by the user (after using the previous tool or any other multiobjective method, possibly an evolutionary algorithm)

• a secondary set of cost functions to be reduced while maintaining quasi Pareto optimality of the first set.

In this report, we use a slightly modified version of the base code permitting to externalise the initial calls to FAST-OAD necessary to build up the local metamodels.

Stochastic MGDA This component provides references to SMGDA, an extension of MGDA applicable to certain stochastic formulations.

The platform was developed by the SED Inria Department (Service for software Experimentation and Development), in particular by T. Kloczko, N. Niclausse, M. Vesin and J.Wintz. The build-up of this platform has required a focus on two software entities, namely MGDA itself, and the server providing its online counterpart experience, hereafter referred to as mgda-www.

Technical considerations on MGDA

The MGDA base code has been developed using fortran 77 (for the core) up to 2003 (see C++ additions below) using the blas and lapack libraries as its sole dependencies. Its uses cmake as a build system and is either packaged as a rpm package, suitable for the red-hat family of linux distributions, or as a self-contained binary archived, suitable for most linux distributions. It is also packaged for macOS and windows platforms.

What follows describes the original Fortran implementation of MGDA along with some C++ additions later helpful to use this codebase in the context of the online platform.

The main code and subroutines are written in Fortran 77, except for one procedure. The execution follows the preparation steps indicated in details in [7] permitting to define the initial local metamodels for constraints and cost functions, including an orthogonal decomposition of the working space on which relies the territory splitting between the virtual players of the Nash games. A procedure that permits to identify the minimum-norm element in the convex hull of given gradients was written in C to permit recursivity.

C++ additions. In order to handle the variety of use cases exposed by the online platform, some features have been added to the original code base, namely to:

• handle a heterogeneous set of options • add some parsing facilities • implement a pilot to handle them all and interface the original fortran procedures

• intertwine these procedures with the pilot using ISO C BINDING.

Fortran's ISO C Binding module provides named constants, types and procedures that are useful in a mixed-language (C and Fortran) program.

Technical considerations on mgda-www

The online user interface over the software stack described before is implemented as a web application leveraging the nodejs asynchronous event-driven JavaScript runtime, for both its frontend part, and its backend part.

Frontend. Presented as a conventional website, each component follow the same structure:

• Synopsis: a generic description of the problem at hand • Run: a set of forms to retrieve inputs and options from the user, actually run the code and present its output

• Example: more in-depth example runs to help the power user, or collaborator with an access on the binaries for its own academic use, to implement or interface its own use case with the mgda procedures.

Backend. The backend is the component that connects the frontend of mgda-www together with the actual code. However, in order to make these calls to the (server side) mgda binaries consistent with the the various step of the wizards implemented by the frontend, the platform implements and leverages a RESTful Application Programming Interface (API). Put simply, a REST API, allows a client (a web browser) to send messages to a server in order for it to perform specific actions and return their result asynchronously to be rendered with the web page.

The RESTful implementation in mgda-www is very basic. It consists in the following routes:

• upload: This allows the user to upload its own dataset, to be fed either to the MGDA or Nash-MGDA procedures.

• compute: This checks and validates the presence of all required files before spawning and executing the relevant procedure call at hand, finally retrieving its output and presenting them either in the frontend or as a download option.

• download: This allows the user to download the output file of the previous remote procedure call. In the standard usage of the platform, the user provides fortran-compatible procedures that permit the separate computation of (i) the primary cost functions, (ii) the secondary cost functions, and (iii) the constraints functions that are called during the computation of the continuum of Nash equilibria to upgrade the metamodels. However here, all three sets are computed together by one call to FAST-OAD, and this is done only once per equilibrium point. (See Section 6.)

Coupling the two platforms

This coupling has been realized in compliance with the instructions of the document [3] by using the associated library call_py_fort.

A FAST-OAD run wraps up an OpenMDAO problem, which is, in a nutshell, the assembly of components that each have input and output variables. Of course, the outputs of some component can be the inputs of some other ones, so that the whole system can be solved. FAST-OAD allows to define the problem to solve (or to optimize) through a configuration file that makes easy to add/remove/replace any component. By doing that, the input data of the problem can be very different from one problem to the other, but FAST-OAD comes with facilities to build the needed input data files. A FAST-OAD problem can be fully run from a command line interface or from the Python API.

Usage of Python API, including pre-processing and post-processing utilities are currently provided through Python notebooks.

A typical run of FAST-OAD uses following types of user files:

Configuration file (.yml): This file defines the OpenMDAO problem by defining:

• what components will be in the problem

• the files for input and output data

• the problem settings

• the definition of the optimization problem if needed (not used here).

Inria

Input/output data files (.xml): These files contain the information of the variables involved in the system model:

• The input file contains the global inputs values required to run all the model. The user is free to modify the values of the variables in order that these new values are considered during a model run.

• The output file contains all the variables (inputs + outputs) values obtained after a model run.

Mission definition file (.yml): This file defines the mission profile that is used for computing aircraft performances using a mass point formulation. It is used only if a flight mission is simulated and none of the bundled mission profiles is used, which is the current case.

To the end of coupling FAST-OAD together with MGDA in order to address the problem at hand, these files have been tailored as follows:

• The configuration file defines the overall aircraft design problem in a classical way: the aircraft is sized though an iterating process devised to achieve the targeted mission while respecting the specified values for input parameters. Simple, low-fidelity models are used for sake of computation speed.

• The input file (.xml) defines the values for all input parameters that are not modified in our process (e.g. passenger count, mission range, various models settings, etc).

• The mission definition file (.yml) defines the mission profile that is used for computing fuel consumption and climb duration. The targeted mission is a climb-cruise-descent sequence, followed by a diversion, a holding phase, and a final descent.

• FAST-OAD launcher (.py): This Python file allows MGDA to launch FAST-OAD. It provides a function that is called by the Fortran program, and launches a FAST-OAD computation.

The function uses a dictionary as argument. This dictionary is populated by the Fortran program (using the set_state function from call_py_fort) prior to calling the Python function. The Fortran program retrieves the FAST-OAD results using the set_state function from call_py_fort.

Optimization test-case

We consider the optimization of six variables that are essential to the design of an aircraft of the Airbus A320-type, to maximize its flight performance measured by three cost functions two of which prioritized under several functional constraints calculated by the simulation software FAST-OAD.

Cost and constraint functions

Primary cost functions

• Fuel mass at take-off (by multiples of 10 4 kg): f 1 ;

• Operating empty weight, OEW (by multiples of 10 4 kg): f 2 .

Secondary cost function

• Ascent time (by multiples of 10 3 s): f 3 .

Functional constraints

• 0.05 ≤ static margin ≤ 0.10

• wing span ≤ 36 (m). These physical variables are all bound by interval constraints, namely:

Physical and dimensionless optimization variables

a i ≤ φ i ≤ b i (i = 1, . . . , 6) (15)
for user-specified {a i } and {b i }. These constraints are handled via the parameterization chosen to define a simple one-to-one relationship between φ i and x i in a way such that:

• x i varies in the interval [-1, 1];

• φ i is a monotone-increasing function of x i from φ i (-1) = a i to φ i (1) = b i .
The interval constraint on φ i is equivalent to maintaining the function

ψ i = φ i - a i + b i 2 b i -a i 2 (16)
in the interval [-1, 1]. This can be accomplished simply via a sine transform:

ψ i = ψ i (x i ) = sin πx i 2 , x i ∈ [-1, 1], (17)
so we let:

φ i (x) = a i + b i 2 + b i -a i 2 sin πx i 2 , x i ∈ [-1, 1]. ( 18
)

Slack variables and numerical formulation of constraints

Both constraints are formulated as equality constraints, say c 1 (x, x 7 ) = 0 for the interval constraint on static margin c 2 (x, x 8 ) = 0 for the inequality constraint on wing span where, for the moment, x ∈ R 6 is the vector of dimensionless optimization variables, whereas (x 7 , x 8 ) are slack variables to be defined. Subsequently, (x 7 , x 8 ) are integrated into the definition of the vector x thereafter handled as an 8-component vector in the numerical process.

The particular form given to (c 1 , c 2 ) should reflect two objectives:

1. enforcing automatically the inequality constraints on state margin and wing span; Concerning the above second requirement, convexity is required to facilitate the back steering of these functions to 0 via a Newton-type iteration on the slack variables when a small violation of the constraint is observed. The functions are also expected to be convex in x, or require only a small convexity fix. Monotonicity is not essential but convenient. It ensures that each constrained physical variable (static margin or wing span) varies in the same sense as its associated slack variable (resp. x 7 or x 8 ).

Interval constraint on static margin The second derivatives of the static margin w.r.t. the optimization variables have been approximated by central differences at the starting point x A . All six axial derivatives were found negative: where the subscript 0 refers to the starting point. Thus, the following dimensionless function that varies in the opposite sense is defined

∂ 2 static : margin ∂x 2 i 0 = (-0.
g 1 = g 1 (x) = 0.10 -static : margin 0.05 = 2 -20 × static : margin
and the constraint is formulated as follows:

0 ≤ g 1 ≤ 1.
and handled simply via a sine function:

c 1 (x, x 7 ) = g 1 (x) -sin γ 1 = 0, γ 1 = 1 -x 7 K 1 ,
in which to ensure the positivity of the sine and its convexity in x 7 at least inside the macro lattice (that is here,

x 7 ∈ [-1, 1]
), the argument γ 1 should be maintained in [0, π]. Actually, the more stringent condition γ 1 ∈ [0, π/2] will now be enforced by the choice of the constant K 1 . For this, let

θ 1 = sin -1 (g 1 ) 0 ∈ [0, π/2]
and assume that θ 1 ≤ π/4; in the case θ 1 > π/4, slightly different formulas would be used with identical efficiency. At x A , the condition c 1 = 0 must be enforced, and this requires that

1 -(x 7 ) 0 K 1 = θ 1 + 2kπ, or π -θ 1 + 2kπ. Let us try 1 -(x 7 ) 0 = K 1 θ 1 .
At the upper limit of the macro lattice, x 7 = 1, and γ 1 = 0. At the lower limit, x 7 = -1, and

γ 1 = 2/K 1 = π/2 provided K 1 = 4/π; then (x 7 ) 0 = 1 -K 1 θ 1 = 1 - 4θ 1 π ∈ [0, 1].
At intermediate points, γ 1 ∈ [0, π/2] as required. Finally:

γ 1 = (1 -x 7 ) π 4 .
Inequality constraint on wing span The axial second derivatives of wing span have been approximated by central differences at the starting point, and all six were found to be negative: This suggests to define c 2 in a way such that it varies in the opposite sense. One lets

∂ 2 wing : span ∂x 2 i 0 = -2.
g 2 = 1 - wing : span 36
and the constraint is formulated as follows:

0 ≤ g 2 ≤ 1.
Mimicking the development of the previous subsection, one lets:

c 2 (x, x 8 ) = g 2 (x) -sin γ 2 , γ 2 = (1 -x 8 ) π 4 ,
and

θ 2 = sin -1 (g 2 ) 0 , (x 8 ) 0 = 1 -K 2 θ 2 = 1 - 4θ 2 π ∈ [0, 1],
in which the assumption θ 2 ≤ π/4 was made.

Design of Experiment

A Design of Experiment is carried out in the preparation phase to calculate a database used to elaborate initial metamodels for the constraint and steering functions. Lattices of points in R 8 are defined, and function values at these points are calculated. The resulting database serves to the construction of the initial quadratic metamodels. More specifically, the data are used to approximate the gradients and the diagonal elements of Hessians by central differencing, and to approximate the other second derivatives by least squares. The lattices are first constructed in R 6 by discretizing the 6 physical variables and five "physical functions" (the 3 cost functions, static margin and wing span) are evaluated at these points by calls to FAST-OAD. The corresponding database is then extended to dimension 8 by discretizing also the slack variables. The latter do not affect the physical functions and only affect the constraint functions (c 1 , c 2 ). As a result, for each 8-component vector x, we dispose of the values of the 3 cost functions (f 1 , f 2 , f 3 ) and 2 constraint functions (c 1 , c 2 ).

The starting values of the physical variables {(φ i ) 0 } are provided by the user, yielding the {(ψ i ) 0 }, from which the components of x A are calculated:

(x i ) 0 = 2 π sin -1 (ψ i ) 0 . ( 19
)
For each component, the distance to the limits is calculated:

d i = min((x i ) 0 + 1, 1 -(x i ) 0 ) (20)
as well as d = min i d i ( 1 -distance to the boundary). A somewhat smaller distance h cut ≤ d is chosen.

Micro lattice in R 6 . This lattice is chosen to approximate at x A the first partial derivatives of various functions and their Hessian diagonal elements. The approximation is made by central differences with stepsize h = 10 -2 . This lattice is thus centered at x A and contains 16 additional points.

Macro and medium-size lattice in R 6 . These lattices are used to compute the non-diagonal elements of Hessians by least-squares, and thus provide a more global information. The macro lattice is centered at 0 ∈ R 6 . For each pair of coordinates (x i , x j ) (15 in total), 8 other points correspond to a uniform discretization of [-1 2 , 1 2 ] 2 , and 8 other points to a uniform discretization of [-1, 1] 2 . This lattice contains 15 × 16 + 1 = 241 points. The medium-size lattice has the same structure but is centered at x A (restricted to R 6 ) and its largest size is h cut instead of 1. This lattice also contains 241 points. Note that by discretizing x i at the points 0, ± 1 2 , ±1, the corresponding sine function values are 0, ±1/ √ 2, ±1 and these correspond to a Chebyshev distribution of points which is the recommended interpolation support in the so-called "best approximation".

In total, micro, medium-size and macro lattices contain 13 + 2 × 241 = 495 points.

Lattices extended to R 8 . In their extensions to R 8 the numbers of datapoints are:

-17 points for the micro lattice -28 × 16 + 1 = 449 points for the macro and medium-size lattices.

This makes a total of 17 + 2 × 449 = 915 data-points, out of which only 495 are actually evaluated by calls to the simulation procedure, and the other ones extrapolated from the parameterization since the slack variables do not affect the physical functions, but only the somewhat artificial constraint functions (c 1 , c 2 ) through explicit formulas.

Weights and least-squares matrix. The least-squares procedure to approximate off-diagonal elements of Hessians has been described for the base algorithm in [7]. Two revisions were brought to this procedure impacting the least-squares matrix and the inversion process. Firstly, since the database of points is not completely centered at x A , the least-squares matrix is no longer scalar and has been recomputed. Secondly, the functions values have been weighted to reduce the influence of the points of the macro lattice that are very remote from x A . The weights were defined to be 1 inside the medium-size lattice and to have the decay of the error-function of the 1 -distance from it. The adjustment of the parameters was not really optimized.

The matrix is calculated and factorized once, and then used several times, firstly for the constraint metamodels, and subsequently for the two steering functions metamodels.

Coordination loop

The coordination loop has been described for the base algorithm in [7]. The discussion is here limited to the handling of the slack variables that has been revised to maintain the constraint functions rigorously to 0 throughout the determination of the continuum of Nash equilibria.

At the start of the procedure to compute a new Nash equilibrium point, the FAST-OAD software is called to evaluate the physical functions, as well as the cost and constraint dimensionless functions. Then the slack variables are updated to drive the cost functions (c 1 , c 2 ) to 0, and this has the effect of upgrading the constraint metamodels. Then begins the coordination between the constrained minimization of the primary steering function f + A by using the subvector u (for fixed v) as the strategy of the virtual player A competing with the unconstrained minimization of the functionf AB by using the subvector v (for fixed u). The convergence of this coordination results in the definition of a new Nash equilibrium point. The procedure is continued until the coordination process fails to find a limit after a prescribed number of iterations, or because the defect of convexity makes the coordination itself inapplicable. Then the plotting procedure is launched using the available computed data.

Design of Experiment

The continuum of Nash equilibria was calculated in three cases, indicated by the symbols L, C and R, on Figure 4 corresponding to three different starting points. To prepare each calculation, a database of 495 points about the starting point was computed each time providing values of (f 1 , f 2 , f 3 ) and (c1, c 2 ) at 915 points according to the process detailed in Section 5. This calculation took about 4hrs on a laptop in the Case C; somewhat less in Case L, and more in Case R.

7.3 Second phase of optimization: Nash equilibria

7.3.1 Trace of continua in the (f 1 , f 2 ) plane
In all three cases, we observe on Figure 4 that the trace of the continuum (in green) is tangent to the Pareto front at the starting point, as it should by theory. After a segment very close to the front, the trace experiences an abrupt change in curvature. The performance achieved by these Inria continua must now be assessed. The numerical process fails at some ε inferior but close to is the limit of convexity ε max calculated by the numerical method. Not only the process has permitted to reduce the climb duration from the starting value, but it has also produced design-points that are more efficient in terms of this variable than all the design-points of the (f 1 , f 2 ) Pareto front.

Three-dimensional visualizations

These visualizations depict the continua of Nash equilibria as curves in the (f 1 , f 2 , f 3 ) space (in green) in Cases L (upper one), C (center) and R (lower one).

In Figure 6, the continua are compared with the values of the triplet (f 1 , f 2 , f 3 ) attained at the points of the Pareto set w.r.t. (f 1 , f 2 ) alone. In Figure 7, the continua are compared with the Pareto front associated with (f 1 , f 2 , f 3 ).

As already mentioned, the Pareto sets were numerically approximated by applying the classical NSGA-II algorithm [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] using the pymoo toolbox [START_REF] Blank | Pymoo: Multi-objective optimization in Python[END_REF]. In the first case, 50 individuals were used, 50 generations computed, corresponding to 2500 calls to FAST-OAD, and the run time on a lap-top was 22 hrs. In the second case, 100 individuals were used, 100 generations computed, corresponding to 10000 calls to FAST-OAD, and the run time on a lap-top was 77 hrs. In the Introduction, it was observed that when the cost functions are convex, the continuous Pareto set can only be enlarged by considering three cost functions instead of two. Here, discrete Pareto sets are considered for which a shift of the Pareto optimal points is possible; additionally, convexity is not guaranteed. Nevertheless the enrichment in red symbols from the first to the second figure is indicative of the enlargement.

The Nash equilibria are observed on Figure 7 to be above but close to the Pareto front associated with the whole set of cost functions, without prioritization. The paths appear directed towards local sinks of the secondary cost function f 3 in the (f 1 , f 2 , f 3 ) Pareto front visualized by the interpolated grid.

Cost Functions along continuum

The variation with ε of the three cost functions f 1 (ε)/f 1 , f 2 (ε)/f 2 and f 3 (ε)/f 3 and the two steering functions f + A (ε) and f B (ε) is represented in Figures 8,9 and 10 in Cases L, C and R respectively. In all three cases, the method is effective up to the end of the computed path which is only a part of the continuum, including the portion after the abrupt change in curvature. These plots illustrate the usual "funnel-shaped" pattern, systematically observed in other cases [START_REF] Désidéri | Direct and adpative approaches to multi-objective optimization[END_REF]7,[START_REF]Prioritized optimization by Nash games: towards an adaptive multi-objective strategy[END_REF][START_REF] Désidéri | Prioritized multi-objective optimization of a sandwich panel[END_REF][START_REF]Adaptation by Nash games in gradient-based multi-objective/multi-disciplinary optimization[END_REF]: the primary steering function admits an initial horizontal tangent and increases moderately with ε, and so do the primary cost functions (f 1 , f 2 ), which remain nearly Pareto-stationary; the secondary steering function admits a tangent with a strictly negative slope, and the secondary cost function f 3 remains very close to it, and decreases importantly. In each case, the designer could favor a point prior to the last computed Nash equilibrium point. Nevertheless, as an indication, Table 1 provides the percentage variation of the three cost functions realized at this point.

Case

f 1 f 2 f 3 L +0.6 +1.1 -25.5 C
+3.9 +4.6 -62.9 R -4.8 +1.4 -65.4

Table 1: Percentage variation of the three cost functions (fuel mass, OEW and climb duration) from starting point to lastly-computed point

Constraints along continuum

As mentioned in Section 6, the slack variables are upgraded after each Nash equilibrium point without alteration of the physical functions and this permits to maintain the constraint functions (c 1 , c 2 ) to exactly 0. As a result, for the computed Nash equilibrium points the interval constraint on static margin and bound on wing span are satisfied throughout. This is evident from The interval constraint on the static margin and the bound constraint on the wing span are satisfied throughout.

Dimensionless and physical optimization variables along the continuum

These variables are observed to remain smooth all along, except for the slack variables at the approach of the convexity break point. We have considered a simplified multi-objective optimization of the flight performance of an aircraft, of Airbus A320 type. The problem involved:

• six optimization variables affecting wing geometry, mean aerodynamic chord, and maximum take-off thrust;

• two cost functions (f 1 , f 2 ) (fuel mass at take-off, Operating Empty Weight) considered primary;

• one cost function f 3 (climb duration), considered secondary;

• an interval constraint on the static margin;

• a bound constraint on the wing span.

The notion of prioritization relied on the split made between primary and secondary cost functions. The performance was calculated by a call to the open-source code FAST-OAD, and the optimization conducted using a customized version of the Nash-MGDA code of the platform mgda.inria.fr . In a first optimization phase, the Pareto set and front associated with the sole primary cost functions were identified by using the classical NSGA-II algorithm. Three design-points were chosen in the Pareto set, two extreme (L and R), and one central (C). For each design-point, a Design of Experiment (DoE) was conducted externally to support the construction of local meta-models.

The second optimization phase was then initiated from either L, C or R. In the customized version of the Nash-MGDA code that relies on the pre-calculated data from the DoE, only one call to the simulation procedure, here FAST-OAD, was executed per computed Nash equilibrium point. Additionally, an improved treatment of the slack variables permitted to satisfy the functional constraints rigorously throughout the entire process.

In all three cases (L, C and R), the method was found effective to reduce substantially the secondary cost function f 3 . In the central case C, the reduction could be of more than 60% while f 1 was diminished, an improvement, and f 2 was increased of less than 5%, a small degradation; both effects partially compensating each other resulted in an even smaller degradation of the Pareto-stationarity condition.

The focus in this report was on the demonstration of the potential of the prioritized optimization approach. Nevertheless, these preliminary numerical experiments permit to draw some sensible conclusions. In particular, it was observed that the mostly influential variables are different depending on the location on the Pareto front of the point chosen to initiate the calculation of the Nash equilibria:

• Case L (large OEW, small fuel mass): Horizontal Tail Aspect Ratio (to be reduced); Wing Taper ratio and Maximum Take-Off Thrust (to be increased);

• Case C (lighter sum at take-off OEW+mass fuel): Wing Taper Ratio (to be reduced); Maximum Take-Off Thrust (to be increased);

• Case R (small OEW, large fuel mass): Wing Aspect Ratio and Taper Ratio (to be increased); Horizontal Tail Aspect Ratio (to be reduced).

Concerning cost efficiency, the following estimated runtimes provide only approximate indications since the computations were not all performed on the same computer.

• One call to FAST-OAD: from 15 s to 45 s in general, depending on the point of evaluation which affects the inner iteration Note that the present optimization test-case was indeed multi-objective, but not multi-disciplinary or multi-point. Thus a unique call to the performance evaluation procedure (FAST-OAD) was sufficient to compute all three cost functions. Nevertheless, the computation of the Pareto front associated with the three cost functions was substantially more costly than with two of them. In multi-disciplinary applications that we target in the long run, the penalty would be even more severe. That future context could thus reveal even more favorable to the mixed approach Paretofront/Nash-games in contexts where the prioritization does make sense. The computations relied on the gradient-free simulation procedure FAST-OAD, thus demonstrating a form of robustness of the prioritized approach. However applications in which the functional gradient is at hand could likely be handled even more efficiently by virtue of a more robust construction of the local metamodels.

Hence, in summary, when a prioritization of the cost functions in primary and secondary cost functions makes sense, it is computationally significantly more economical and effective to limit the computation of the Pareto front to the primary cost functions and to compute rapidly one or several continua of Nash equilibria that reduce the secondary cost functions, while best preserving the Pareto-optimality of the essential, primary cost functions. The Nash equilibria have been found close to the Pareto front associated with the whole set of cost functions, without prioritization, thus providing excellent designs.

In the future, we plan to realize a more complete and realistic experiment involving more than one secondary cost functions.

Inria

A Nested Pareto sets and fronts

Very intuitively, as the number of cost functions is increased, one expects the Pareto-optimality criterion to discriminate less, and the Pareto set and front to be enlarged without loss. This question is of some relevance since it is natural in case of concurrent engineering to introduce additional criteria gradually. In [START_REF] Allmendinger | What if we increase the number of objectives? theoretical and empirical implications for many-objectve combinatorial optimization[END_REF], this issue has been considered theoretically and empirically in the context of combinatorial optimization. However, we were not able to find a proper "textbook reference" for this nesting property, likely to be classical, that would apply more appropriately to the context of continuous optimization. For sake of completeness, and in particular to support (3), a justification is provided here for the basic result according which, in a continuous and convex setting, Pareto sets and fronts are generally nested as the number of cost functions is increased.

Consider a set of k + 1 (k ≥ 2) continuous and convex cost functions {f 1 (x), . . . , f k+1 (x)} defined over the closed and convex domain Ω ad ⊆ R n . Let {S k , F k } be the Pareto set and front associated with the first k cost functions, and {S k+1 , F k+1 } the Pareto set and front associated with the whole set of cost functions.

Let Fk be an arbitrary element of F k , image of at least one point x in the kth-dimensional function space, that is Fk = ( f1 , f2 , . . . , fk ) = (f 1 (x), f 2 (x), . . . , f k (x)).

(21)

Let Fk+1 = ( Fk , f k+1 (x)) be the image of x in the k +1st-dimensional function space. Then, either Fk+1 ∈ F k+1 , or there exists x ∈ Ω ad that dominates x in efficiency w.r.t. {f 1 , f 2 , . . . , f k+1 }, i.e.:

f 1 (x ) ≤ f 1 (x) = f1 , . . . f k (x ) ≤ f k (x) = fk , f k+1 (x ) ≤ f k+1 (x) = fk+1 , (22)
and at least one of the above inequalities is strict. The latter cannot be one of the first k since otherwise x would also dominate x in efficiency in S k . Therefore, it must be the last one:

f 1 (x ) = f 1 (x) = f1 , . . . f k (x ) = f k (x) = fk , f k+1 (x ) < f k+1 (x) = fk+1 . ( 23
)
Since all the functions are continuous and the domain is closed, f k+1 (x) subject to the constraints defined by the above equalities admits a minimum. Without loss of generality, we may assume that x realizes this minimum. Then we claim that x ∈ S k+1 . Indeed otherwise there would exist an element x ∈ Ω ad that would dominate x in efficiency in S k+1 , i.e.:

f 1 (x ) ≤ f 1 (x ) = f 1 (x) = f1 , . . . f k (x ) ≤ f k (x ) = f k (x) = fk , and f k+1 (x ) ≤ f k+1 (x ) < f k+1 (x) = fk+1 , (24)
and at least one of the above large inequalities should be strict. The latter could not be one of the first k since otherwise x would dominate x in efficiency in S k . Hence, the following should hold

f 1 (x ) = f 1 (x ) = f 1 (x) = f1 , . . . f k (x ) = f k (x ) = f k (x) = fk , f k+1 (x ) < f k+1 (x ). ( 25
)
But this is in contradiction with the fact that x , among all the elements for which (23) holds, minimizes f k+1 (x ). To reject this contradiction, we conclude that indeed x belongs to S k+1 , and its image in the k + 1st-dimensional function space

F k+1 = (f 1 (x ), . . . , f k+1 (x )) (26)
belongs to F k+1 . But F k+1 = ( Fk , f k+1 (x)). Hence in the k + 1st-dimensional function space, Fk+1 admits the projection F k+1 in the direction of the f k+1 axis that belongs to the front F k+1 . In conclusion, given the arbitrary element Fk ∈ F k , and its extension Fk+1 = ( Fk , f k+1 (x)) to the k + 1st-dimensional function space, either Fk+1 belongs to the front F k+1 , or it admits a projection in the direction the f k+1 axis onto the front F k+1 .

We presently examine the inclusion of Pareto sets. For this, we further make the hypothesis that the existence of a line segment of the domain over which {f 1 (x), f 2 (x), . . . , f k (x)} are all constant is excluded. Then:

1. Uniqueness: Fk ∈ F k admits no antecedent in S k other than x.

However, none of these inequalities can be strict, since otherwise x α would dominate x in efficiency in S k . Hence these inequalities can be replaced by equalities, and we conclude that the first k cost functions are constant over the line segment [x, y]. But this eventuality has been excluded by assumption. Therefore the existence of y = x is rejected: x is unique. Now, let x be an arbitrary element of S k and make the tentative hypothesis that x / ∈ S k+1 . Then, there exists an admissible element y that dominates x in efficiency w.r.t. {f 1 , f 2 , . . . , f k+1 }, that is:

f 1 (y) ≤ f 1 (x), f 2 (y) ≤ f 2 (x), . . . , f k+1 (y) ≤ f k+1 (x), (28)
and at least one of these inequalities is strict. That inequality cannot be one of the first k since otherwise y would dominate x in efficiency also w.r.t. {f 1 , f 2 , . . . , f k } which is excluded since x ∈ S k . Hence in (28):

f 1 (y) = f 1 (x), f 2 (y) = f 2 (x), . . . , f k (y) = f k (x), (29)
and f k+1 (y) < f k+1 (x). But, we have shown that a point of the Pareto front of the functions {f 1 , f 2 , . . . , f k } is the image of a unique element of S k , and since y = x, we arrive at a contradiction. The contradiction is waived by rejecting the tentative hypothesis (x / ∈ S k+1 ), and we conclude that instead x ∈ S k+1 . Since x is arbitrary in S k , this conclusion implies that S k ⊆ S k+1 and the remaining follows.

Remark 1

Excluding situations in which the first k cost functions would be constant over some line segment is not very restrictive. Nevertheless, the following example was fabricated to exhibit a situation in which the exclusion does not hold:

x = (x, y, z) ∈ Ω ad = R 3 , f 1 (x) = φ(x, y), f 2 (x) = ψ(x, y), f 3 (x) = z 2 , {φ, ψ} convex. (30) Then in R 3 , S 2 is an infinite cylinder invariant by translation parallel to the z axis, while S 3 is limited to its trace in the plane z = 0. Of course, S 2 is larger than S 3 only if we consider Ω ad to be R 3 , which we would not do in a study limited to the pair {f 1 , f 2 }.

  f j = 0 (m=2 on plot) 2nd phase: Adaptation to 2ndary cost functions (j > m): continuum of Nash equilibria

Figure 1 :

 1 Figure 1: Primary Pareto Front (in red) and Continuum of Nash Equilibria (in green) in the case of m = 2 primary cost functions (f 1 , f 2 ) and an arbitrary number of secondary cost functions

Figure 2 :

 2 Figure 2: Index page for the first component of the platform

Figure 3 :

 3 Figure 3: Forms to allow interacting with the server side binaries

Six optimization variables subject to interval bound constraints are considered: 1 .

 1 distance from nose to 25% of wing Mean Aerodynamic Chord (MAC), φ 1 = W M AC25 ∈ [16, 18], to be controlled by x 1 , 2. wing aspect ratio, φ 2 = W AR ∈[START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF][START_REF] Désidéri | Computational Intelligence in Aerospace Sciences[END_REF], to be controlled by x 2 , 3. horizontal tail aspect ratio, φ 3 = HT AR ∈ [1.[START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF][START_REF] Désidéri | Quasi-Riemannian multiple gradient descent algorithm for constrained multiobjective differential optimization[END_REF], to be controlled by x 3 , 4. wing taper ratio, φ 4 = W T R ∈ [0, 1], to be controlled by x 4 , 5. wing sweep angle at 25%MAC, φ 5 = W SW EEP 25 ∈ [20, 30], to be controlled by x 5 , 6. maximum take-off thrust, φ 6 = M T OT ∈ [100000, 140000], to be controlled by x 6 .

Inria 2 .

 2 defining c 1 (resp. c 2 ) as a smooth, convex and monotone-increasing function of x 7 (resp. x 8 ) at least in the interval [-1, 1].

7. 1 RFigure 4 :

 14 Figure 4: Pareto front associated with the primary cost functions (f 1 , f 2 ) (red symbols); traces of three continua of Nash equilibria to reduce f 3 (green)

7. 3 . 2

 32 Climb duration f 3 vs fuel mass f 1 along the continuum

Figure 5 :

 5 Figure 5: Climb duration f 3 vs fuel mass f 1 along the continuum -Case C

2 f 1 f 3 Figure 6 : 3 Figure 7 :

 213637 Figure 6: Values of (f 1 , f 2 , f 3 ) attained by the points of discrete Pareto set associated with (f 1 , f 2 ) alone (red symbols), and the three continua of Nash equilibria initiated from it in Cases L (upper one), C (center) and R (lower one) (green curves)

  Case 1 -continuous cruise altitude -start at point L

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Cost functions along continuum as functions of ε -Case L

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: Dimensionless and physical optimization variables along the continuum -Case L Mostly influential: Horizontal Tail Aspect Ratio (to be reduced); Wing Taper ratio and Maximum Take-Off Thrust (to be increased)

••

  1st Phase: Two-objective (f 1 , f 2 ) Pareto front: about 21.5 hrs by NSGA-II (2500 calls to FAST-OAD v. 1.3.4) Design of Experiment about x A (495 points): about 4 hrs (about 30 s per call), in Case C • 2nd Phase: Continuum alone (22 equilibrium points): about 7 mn runtime, in Case C • Three-objective (f 1, f 2, f 3) Pareto front: about 77 hrs by NSGA-II.

2 .

 2 Inclusions: S k ⊆ S k+1 , and the image of S k in the k + 1st-dimensional function space is included in F k+1 .Proof: Suppose that another point y of the Pareto set S k realizes the same first k cost function values as x, and let x α = (1 -α)x + αy (0 ≤ α ≤ 1) be an element of the line segment [x, y] which lies entirely in the assumed convex admissible domain Ω ad . Then, by convexity of the cost functions: ∀j ≤ k, f j (x α ) ≤ (1 -α)f j (x) + αf j (y) = fj .

  100969, -32.5971, -32.6567, -57.5602, -71.5263, -74.8816)

  064 × 10 -2 , -15.9242, -1311.66, -1311.92, -1424.07, -1522.41 .

Table 2

 2 which provides the first five and last five sets of values in Case C.

	Nash	φ 1	φ 2	φ 3	Static	Wing
	equil.	(fuel)	(OEW)	(duration) Margin	Span
	point (×10 4 kg) (×10 4 kg)	(×10 3 s)	(×10 -2 )	(m)
	0	2.0047	3.8497	2.6175	5.7333	27.6631
	1	1.9975	3.8541	2.5542	5.7322	27.8288
	2	1.9933	3.8570	2.5147	5.7337	27.9297
	3	1.9904	3.8591	2.4830	5.7347	27.9971
	4	1.9885	3.8607	2.4570	5.7356	28.0447
	5	1.9871	3.8622	2.4308	5.7365	28.0792
	. . .	. . .	. . .	. . .	. . .	. . .
	18	2.0090	3.9126	1.5560	5.8797	28.0927
	19	2.0197	3.9300	1.4082	5.8301	28.0798
	20	2.0367	3.9567	1.2423	5.8767	28.0647
	21	2.0618	3.9940	1.0757	5.9505	28.0492
	22	2.0953	4.0420	0.9233	6.0008	28.0247
	∆ max	+3.9%	+4.6%	-62.9%

Table 2 :

 2 Physical cost functions and constraints along the continuum in Case C

RR n°9490

https://www.isae-supaero.fr/en/

https://www.onera.fr/en RR n°9490

Inria

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399