# Prioritized multi-objective optimization of an aircraft flight performance based on Nash games from preponderant Pareto-optimal points

Jean-Antoine Désidéri*

Julien Wintz*

Mickaël Binois*

Nathalie Bartoli†

Christophe David†

Sébastien Defoort†

## Abstract

A hierarchical multi-objective optimization process is applied to the performance of an Airbus-A320-type aircraft. Six design variables are used to calibrate wing geometry, mean aerodynamic chord and take-off thrust potential. Functional constraints are enforced (interval bound on static margin and upper bound on wing span). The system is subject to the classical Breguet laws of Flight Mechanics and to the numerical integration of the ascent phase. It is thus an Ordinary-Differential-Equation-constrained system and it is evaluated by the open-source software FAST-OAD. Three cost functionals are minimized: fuel mass at take-off, operational empty weight, and ascent duration. The first two, considered of preponderant importance are prioritized, and the third only secondary. The method permits, at moderate numerical cost, to identify quasi-Pareto-optimal by Nash games originating from Pareto-optimal solutions of the sole prioritized cost functions.

## 1 Introduction

In concurrent engineering, the numerical identification of the Pareto front by a robust method such as NSGA-II [4] is often a highly demanding numerical task, the cost of which may increase rapidly with the number of cost functions accounted for. However in many real design situations, a subset of the cost functions, also known as objective functions in the literature, has preponderant importance, and this incites to a hierarchical, or prioritized treatment. We have proposed to compute both the Pareto set and Pareto front for these primary cost functions only, and to account subsequently for the secondary cost functions by Nash games [9, 7]. The efficacy and cost efficiency of the prioritized approach are presently assessed in a case of multi-objective optimization test-case aimed to be a representative preliminary design phase exercise.

<footnote>*Université Côte d'Azur, Inria, CNRS, LJAD, Sophia Antipolis, France jean-antoine.desideri@inria.fr</footnote>

<footnote>†ONERA / DTIS, Université de Toulouse, France nathalie.bartoli@onera.fr</footnote>
The precise numerical construction is recalled in Section 2 where the convergence result is stated. The proposed algorithm is described in Section 3 from an implementation point of view. In Section 4, the methodology is applied to an aircraft performance optimization test-case involving three objective functions, The numerical results are commented in detail including a qualitative comparison of numerical cost efficiency depending on whether NSGA-II or Bayesian optimization is used to compute Pareto fronts. Finally, conclusions and perspectives are drawn in Section 5.

# 2 Theoretical background

## 2.1 Mathematical problem definition and objective

One considers a multi-objective differentiable optimization problem involving a set of $M$ cost functions $\{f_j(\mathbf{x})\} (j = 1, \dots, M; \mathbf{x} \in \Omega_a \subseteq \mathbb{R}^N)$, the first $m$ of which $(1 \le m < M)$ are declared preponderant and named primary, and the remaining ones $(j = m + 1, \dots, M)$, secondary; all are subject to $K$ equality constraints, $c_k(\mathbf{x}) = 0 (k = 1, \dots, K)$. Objective and constraint functions are all smooth, say $C^2(\Omega_a)$.

A first phase of optimization (omitted here) has led to the knowledge of a point $\mathbf{x}_A^*$ in the open interior of the admissible domain $\Omega_a$. This point is assumed to be Pareto-optimal for the constrained minimization of the sole primary cost functions.

The second phase of optimization is aimed to construct a continuum of Nash equilibria $\{\bar{\mathbf{x}}_\varepsilon\}$ parameterized by $\varepsilon$ ($\varepsilon \ge 0$) in way such that: firstly, $\bar{\mathbf{x}}_0 = \mathbf{x}_A^*$ (compatibility), and secondly, as $\varepsilon$ increases but remains sufficiently small, the Pareto-optimality condition of the primary cost functions remains $O(\varepsilon^2)$, whereas the secondary cost functions decrease linearly, or faster.

## 2.2 Preliminary calculations

The preliminary calculations have been defined in full in [8][7]. They are made of three main steps:

* Definition of a convex "Primary Steering Function" $f_A^+(\mathbf{x})$

* Definition of a "Territory Splitting"

* Definition of a "Secondary Steering Function" $f_B(\mathbf{x})$.

The cost functions are all assumed to be strictly positive. This may require the preliminary application of an adequate exponential transform. This permits us to consider gradients in the logarithmic form, and examine the variations of these functions in terms of proportions.
The Pareto-stationarity condition satisfied at $\mathbf{x} = \mathbf{x}_A^*$ is expressed as follows:

$$ \sum_{j=1}^{m} \alpha_j^* \frac{\mathbf{P} \nabla f_j^*}{f_j^*} = 0 \qquad (1) $$

where the superscript * is used throughout to indicate an evaluation at $\mathbf{x} = \mathbf{x}_A^*$; $\nabla$ is the gradient operator w.r.t. $\mathbf{x}, \mathbf{P}$ is the projection operator onto the linear subspace orthogonal to all constraint gradients $\{\nabla c_k^*\} (k = 1, \dots, K)$ (identified by application of the QR algorithm), and $\{\alpha_j^*\} (j = 1, \dots, m)$ are the known coefficients of a convex combination. Then the primary steering function is defined as follows:

$$ f_A^+(\mathbf{x}) = \sum_{j=1}^{m} \alpha_j^* \frac{f_j(\mathbf{x})}{f_j^*} + \frac{c}{2} \|\mathbf{x} - \mathbf{x}_A^*\|^2 \qquad (2) $$

where $c$ is a nonnegative constant adjusted to enforce the local convexity of $f_A^+(\mathbf{x})$ and the associated Lagrangian. Then, by virtue of Eq. (1), the primary steering function is stationary and convex at $\mathbf{x}_A^*$, and therefore achieves a local minimum.

The territory splitting is based on the orthogonal decomposition of the following $N \times N$ reduced Hessian matrix [8, 10]:

$$ \mathbf{H}'_A = \mathbf{P} \mathbf{H}^{+,*}_A \mathbf{P} = \mathbf{\Omega} \mathcal{H} \mathbf{\Omega}^t \qquad (3) $$

where $\mathbf{H}_A^{+,*} = (\nabla^2 f_A^*)^* = \mathbf{H}_A^* + c \mathbf{I}_N$ and $\mathbf{H}_A^* = \nabla^2 f_A^*$. Matrix $\mathbf{H}_A'$ is positive semi-definite, of rank $N-K$ exactly since matrix $\mathbf{H}_A^{+,*}$ is strictly-positive definite by convexity fix. The eigenvalues in the diagonal matrix $\mathcal{H}$ are real and strictly positive. The eigenvectors, column-vectors of the orthogonal matrix $\mathbf{\Omega}$, are conventionally arranged in a special order: the first $K$ are associated with the directions of the constraint gradients (null space of $\mathbf{P}$), and the remaining ones are ordered by decreasing eigenvalue. In this way, the tail eigen-modes are associated with the smaller sensitivities of the primary steering function. The splitting is associated with the following change of variables:

$$ \mathbf{x} = \mathbf{x}_A^* + \mathbf{\Omega} \begin{pmatrix} \mathbf{u} \\ \mathbf{v} \end{pmatrix} := \mathbf{X}(\mathbf{u}, \mathbf{v}) \qquad (4) $$

where $\mathbf{u} \in \mathbb{R}^{N-p}, \mathbf{v} \in \mathbb{R}^p$, and $p$ is an adjustable integer satisfying $1 \le p < N-K$. A large $p$ gives the secondary cost functions a great weight. In this setting, the following notations are used:

$$ \mathbf{\Omega} = (\mathbf{\Omega}_{\mathbf{u}} \ \mathbf{\Omega}_{\mathbf{v}}), \quad \mathbf{\Omega}_{\mathbf{u}} = \begin{pmatrix} \mathbf{\Omega}_{\mathbf{u}\mathbf{u}} \\ \mathbf{\Omega}_{\mathbf{v}\mathbf{u}} \end{pmatrix}, \quad \mathbf{\Omega}_{\mathbf{v}} = \begin{pmatrix} \mathbf{\Omega}_{\mathbf{u}\mathbf{v}} \\ \mathbf{\Omega}_{\mathbf{v}\mathbf{v}} \end{pmatrix} \qquad (5) $$
and scaled logarithmic gradients of the secondary cost functions are calculated:

$$g_j = S^{-\frac{1}{2}} \frac{\nabla_v f_j^*}{f_j^*} = S^{-\frac{1}{2}} \frac{\Omega_v^t \nabla f_j^*}{f_j^*} \quad (j = m+1, \dots, M) \quad (6)$$

where the scaling matrix $S = (\nabla \nabla_{vv}^t f_A^*)^*$ is the lower $p \times p$ diagonal block of matrix $\mathcal{H}$ of Eq. (3), hence it is guaranteed to be positive-definite by convexity fix. 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) [8][7], one identifies the minimum-norm element $\omega_B^*$ in the convex hull of the gradients $g_j$:

$$\omega_B^* = \sum_{j=m+1}^{M} \alpha_j^* g_j \quad (7)$$

The coefficients $\{\alpha_j^*\}$ of this convex combination are used to define the secondary steering function $f_B(\mathbf{x})$, thus completing the preliminaries:

$$f_B(\mathbf{x}) = \sum_{j=m+1}^{M} \alpha_j^* \frac{f_j(\mathbf{x})}{f_j^*} \quad (8)$$

## 2.3 Nash game formulation and theoretical results

One defines a continuation parameter $\epsilon$ ($0 \le \epsilon \le 1$) and the auxiliary cost function:

$$f_{AB}(\mathbf{x}) = f_A^+(\mathbf{x}) + \epsilon (f_B(\mathbf{x}) - f_A^+(\mathbf{x})) \quad (9)$$

For every fixed value of $\epsilon$, two virtual players, Player A and Player B, respectively controlling the sub-vectors $\mathbf{u}$ and $\mathbf{v}$ in Eq. (4), are involved in a Nash game defined by the following concurrent strategies:

- **Strategy A:** Player A attempts to minimize $f_A^+(\mathbf{X}(\mathbf{u}, \mathbf{v}))$ by the strategy $\mathbf{u}$, subject to the equality constraint $c(\mathbf{X}(\mathbf{u}, \mathbf{v})) = 0$, and by accounting for Player B's fixed strategy $\mathbf{v}$, concurrently with:

- **Strategy B:** Player B attempts to minimize $f_{AB}(\mathbf{X}(\mathbf{u}, \mathbf{v}))$ by the strategy $\mathbf{v}$, subject to no constraints, but by accounting for Player A's fixed strategy $\mathbf{u}$.

The following results have been established [8][7]:
• The Nash equilibrium point $\bar{x}_{\epsilon}$ exists for all sufficiently small $\epsilon$, and

$$ \bar{x}_0 = \mathbf{x}_A^* \qquad (10) $$

an essential property of the formulation referred to as compatibility or consistency.

• If $\phi_A(\epsilon) = f_A^+(\bar{x}_\epsilon)$, $\phi'_A(0) = 0$ and

$$ f_A^+ = 1 + O(\epsilon^2) \qquad (11) $$

which reflects the fact that, as $\epsilon$ varies, the Pareto-optimality of the primary cost functions is degraded by $O(\epsilon^2)$ only.

• If $\phi_B(\epsilon) = f_B(\bar{x}_\epsilon)$, $\phi'_B(0) = -\sigma_B$, where $\sigma_B = \|\omega_B^*\|^2 \ge 0$. The secondary steering function $f_B(\mathbf{x})$ decreases linearly as $\epsilon$ increases, and by virtue of properties of the MGDA construction, all the secondary cost functions do as well [8][7].

Figure 1a illustrates this theoretical result, depicting in particular the continuum of Nash equilibria tangent at $\mathbf{x}_A^*$ to the primary Pareto front.

<figure><img src="image_1.png" /><figcaption>Figure 1: Left: sketch of primary Pareto front and a representative continuum of Nash equilibria tangent to it (a); right: computed primary Pareto front and three continua of Nash equilibria denoted by L, C or R (b).</figcaption></figure>
### 3 Numerical implementation, constraints and coordination loop

The present method can be applied directly to simple cases defined analytically. The software provided by the platform http://mgda.inria.fr can be used for this purpose according to procedures detailed in [6]. In the more general case in which the formal derivation is not entirely at hand, we proceed by applying the method to local metamodels some of which are upgraded in the course of the computational process. A customized variant of the software is then employed. We consider

- a split of the optimization variables $\mathbf{x} = (\mathbf{x}_{\phi}, \mathbf{x}_{s}) \in \mathbb{R}^{N}$ in which each component of $\mathbf{x}_{\phi} \in \mathbb{R}^{n}$ is in a one-to-one correspondence with one of the $n$ natural physical variables of the problem, and those of $\mathbf{x}_{s} \in \mathbb{R}^{s}$ are slack variables ($N = n + s$);

- “physical functions” of $\mathbf{x}_{\phi}$ including cost functions and auxiliary functions in which the “physical constraints” of the problem are expressed; an evaluation procedure to compute these physical functions in terms of $\mathbf{x}_{\phi}$ is at disposal; the number of calls to this procedure is our measure of computational effort;

- “constraint functions” $\{c_k(\mathbf{x})\} (k = 1, \dots, K)$ devised to enforce the physical constraints by the somewhat artificial condition $c_k(\mathbf{x}) = 0, \forall k$; these functions are locally represented by metamodels; once a candidate Nash equilibrium point ($\mathbf{x}_{\phi}, \mathbf{x}_{s}$) has been determined, the satisfaction of the physical constraints by $\mathbf{x}_{\phi}$ is tested; if the test is positive, in preparation of the determination of the next equilibrium point, the slack variables $\mathbf{x}_{s}$ are adjusted to reset the pointwise values of these metamodels to 0, upgrade the metamodels, and avoid accumulation of even minute errors on the $\{c_k(\mathbf{x})\}$; otherwise the process interrupts.

A Design of Experiment (DoE) is firstly conducted over a dataset of points in physical space ($\mathbb{R}^n$) combining a micro lattice, a medium-size lattice and a macro lattice. The micro lattice serves to approximate first derivatives of functions of interest at $\mathbf{x} = \mathbf{x}_A^*$, and the diagonal elements of their Hessians by central differencing. It is made of $2n$ points at $\pm h$ (for a small $h$) about $\mathbf{x}_A^*$ over each axis. The macro lattice and medium-size lattice have similar structure. They are made of 16 points per pair of coordinates $\{x_i, x_j\}$. The macro lattice is centered at the origin and extends over the entire admissible set. The medium-size lattice is centered at $\mathbf{x}_A^*$ and its sides are adjusted to lie strictly within the domain. Altogether the macro lattice and the medium-size lattice provide a more global domain representation in view of approximating the off-diagonal elements of Hessians by least squares. In total, the evaluation of the entire dataset requires $16n(n-1) + 2n + 1$ additional calls to the simulation software. The dataset is then extended to
dimension $N$ by discretizing $\mathbf{x}_s$ only and by evaluating the $\{c_k(\mathbf{x})\}$ functions correspondingly; this evaluation is fast since the physical part has already been computed, and the remaining part is explicit. Quadratic metamodels are constructed once for all for the two steering functions. Similar quadratic metamodels are also constructed at first for the constraint functions, but these are upgraded throughout the computational process of Nash equilibria. We emphasize that this upgrade is essential to preserve the Pareto-stationarity condition. The Hessian of $f_A^+$ and first derivatives of constraint functions are approximated by central differencing. This serves to define the orthogonal decomposition in Eq. (4).

At start ($\mathbf{x} = \mathbf{x}_A^*$), the physical constraints are satisfied exactly since the assumed Pareto optimality of $\mathbf{x}_A^*$ is subject to this condition, and the slack variables have been initialized accordingly. After identifying a number of Nash equilibrium points (outer loop), the metamodels for the constraint functions $\{c_k(\mathbf{x})\}$ are locally reset to 0; then, a situation similar to initial applies at $\mathbf{x} = \bar{\mathbf{x}}_\epsilon$ for a certain $\epsilon > 0$. The search for the next equilibrium point then proceeds with a coordination iteration for the pair $(\mathbf{u}, \mathbf{v})$. This inner iteration begins with the constrained minimization of the primary steering function $f_A^+$ by means of the subvector $\mathbf{u}$ (for fixed $\mathbf{v}$) in compliance with the strategy of virtual Player A (this minimization is carried out by a Newton-type algorithm); then follows the unconstrained minimization of the function $f_{AB}$ by means of the subvector $\mathbf{v}$ (for fixed $\mathbf{u}$) in compliance with the strategy of virtual Player B; this is realized by inversion of the linear system of optimality. The coordination is judged converged when the variation in $\|\mathbf{v}\|$ has achieved a value less than a small tolerance. Then the pair $(\bar{\mathbf{u}}, \bar{\mathbf{v}})$ defines the candidate equilibrium point $(\mathbf{x}_\phi, \mathbf{x}_s)$. The test on physical constraints is made. If it is positive, the slack variables $\mathbf{x}_s$ are corrected, usually by minute variations, to enforce the rigorous satisfaction of the condition $c_k = 0, \forall k$. These small variations affect the slack subvector $\mathbf{x}_s$ only, and have no impact on the physical functions. This completes the definition of the new Nash equilibrium point. The constraint metamodels are then upgraded and the outer loop is continued. The process terminates when the constraints are violated, or when the determination of $\bar{\mathbf{u}}$ fails in the approach of the convexity limit, or because the coordination process has not converged after a prescribed number of iterations.

# 4 Application to an aircraft optimization exercise

## 4.1 Test case definition

We consider the optimization of six design variables impacting wing and tail geometry, engine thrust and longitudinal stability of an Airbus A320-type aircraft to reduce two primary cost functions, and one secondary, subject to functional constraints. The test case is meant to be representative of a process taking place in a preliminary aircraft design phase.
Six design optimization variables in physical dimensions (n = 6):

* $\phi_1$: distance from nose to 25% of wing Mean Aerodynamic Chord (MAC) (denoted WMAC25 on figures); to be in $[a_1, b_1] = [16, 18]$;

* $\phi_2$: wing aspect ratio (denoted WAR); to be in $[a_2, b_2] = [5, 11]$;

* $\phi_3$: horizontal tail aspect ratio (denoted HTAR); to be in $[a_3, b_3] = [1.5, 6]$;

* $\phi_4$: wing taper ratio (denoted WTR); to be in $[a_4, b_4] = [0, 1]$;

* $\phi_5$: wing sweep angle at 25% of mean aerodynamic chord (denoted WSWEEP25); to be in $[a_5, b_5] = [20, 30]$;

* $\phi_6$: maximum take-off thrust (denoted MTOT); to be in $[a_6, b_6] = [100000, 140000]$.

The above physical variables are the components of the vector $\mathbf{x}_{\phi} = \{\phi_i\}$ and are parameterized via dimensionless variables $\{x_i\}$ according to

$$ \phi_i(\mathbf{x}) = \frac{1}{2}(a_i + b_i) + \frac{1}{2}(b_i - a_i) \sin\left(\frac{\pi x_i}{2}\right), \quad x_i \in [-1, 1] \quad (i = 1, \dots, 6) \quad (12) $$

Two primary cost functions ($m = 2$) are chosen:

* fuel mass, $f_1$ (by multiples of $10^4$ kg);

* operating empty weight (OEW), $f_2$ (by multiples of $10^4$ kg).

Physical constraints are considered:

* geometrical bound: wing span $\le 36$ (m);

* functional constraint for longitudinal stability: $0.05 \le$ static margin $\le 0.10$.

One secondary cost function ($M = 3$) is added:

* duration of ascent to cruise altitude, $f_3$ (by multiples of $10^3$ s).

A first slack variable $x_7$ to enforce the interval constraint on the static margin is used. Defining the following function observed to be locally convex

$$ g_1 = g_1(\mathbf{x}) = \frac{0.10 - \text{static margin}}{0.05} \quad (13) $$

the constraint writes $0 \le g_1 \le 1$ and is enforced by letting

$$ c_1 = g_1 - \sin \gamma_1 = 0, \quad \gamma_1 = (1 - x_7) \frac{\pi}{4}, \quad x_7 \in [-1, 1]. \quad (14) $$

A second slack variable $x_8$ ($s = 2$) is introduced to enforce the bound constraint on wing span.
The function $g_2$ observed to be locally convex is defined:

$$g_2 = 1 - \frac{\text{wing span}}{36} \qquad (15)$$

So the constraint writes $0 \le g_2 \le 1$ and is enforced by letting

$$c_2 = g_2 - \sin \gamma_2 = 0, \quad \gamma_2 = (1 - x_8) \frac{\pi}{4}, \quad x_8 \in [-1, 1]. \qquad (16)$$

In this setting, the dimensionless admissible domain is $[-1, 1]^8$.

## 4.2 Evaluation and optimization software

The designs are evaluated by means of the FAST-OAD open-source software developed by ONERA and ISAE-SUPAERO [3] which combines the classical Breguet laws with a time integration of the ascent phase, and performs a full resizing (dimensions, masses) of the aircraft for each set of design variables through a multidisciplinary coupled approach. Hence, this evaluation is constrained by an ordinary-differential equation. The prioritized optimization is conducted by a customized version of the INRIA Nash-MGDA software [5]. From a point of Pareto-optimality of the sole pair $(f_1, f_2)$, considered preponderant, the climb duration $f_3$ is to be reduced along a continuum of Nash equilibria.

## 4.3 Numerical results using NSGA-II

The first phase of optimization consisted in establishing the “Primary Pareto Front” (PPF) associated with the minimization of $(f_1, f_2)$ alone under constraints. This was done firstly by the genetic algorithm NSGA-II [4] using a population of 50 individuals. A well-converged front was obtained using an indicative budget of 2500 calls to the FAST-OAD software, at a computation cost of about 22 hrs on a PC. Then, three sets of Nash equilibria were computed, each one originating from a point of the discrete PPF and appearing on Figure 1b as a green path denoted L, C or R. The three paths are tangent to the PPF and remain close to it initially over a variable length. On Figure 2, the paths are represented in the $(f_1, f_2, f_3)$ space. These results have been described in details in [11], and we now focus our discussion mostly on the central case C. In this case, the DoE involved some 500 preliminary evaluations and this required about 4.5 hrs of computing time. Figure 3a illustrates the "funnel-shaped" convergence plot typical pattern:

* above the $\varepsilon$-axis, the primary steering function $f_A^+$ (black solid line), and below it, still above the axis, the primary cost functions $(f_1, f_2)$ (red and orange solid lines), all three
<figure><img src="image_2.png" /><figcaption>Figure 2: Primary Pareto front and three continua of Nash equilibria represented in ($f_1, f_2, f_3$).</figcaption></figure>

tangent to the axis by virtue of the initial Pareto-stationarity; both cost functions are only
very moderately degraded (less than 5% at endpoint);

* below the $\varepsilon$-axis, the secondary steering function $f_B$ (black hatched line) and the secondary cost function (blue solid line) have the same initial negative slope, confirming the linear decay; $f_3$ follows the metamodel $f_B$ very closely; at endpoint $f_3$ has been importantly reduced (more than 60%).

Figure 3b indicates how does the optimization modify the physical variables with $\varepsilon$, in Case C: the most influential variables have been wing taper ratio (highly reduced), wing aspect ratio (slightly reduced), and take-off maximum thrust (increased). The computation of the continuum of 22 equilibrium points in Case C required about 7 mns on a laptop and has permitted the user interaction. Along the path, the functional constraints have been enforced rigorously [11].

Cases L and R corresponded to very different, and extreme, initial designs and the optimization resulted in very different modifications. Compared to the central case (C), the variations were briefly the following [11]:

- Case L: the initial design corresponded to a heavier aircraft with less embarked fuel. The optimization resulted in an increase of the maximum take-off thrust and small geometrical changes (mostly a reduction of the horizontal tail aspect ratio).

- Case R: the initial design corresponded to a lighter aircraft with more embarked fuel. The optimization resulted mostly in more sensible geometrical changes (increase of horizontal
<figure><img src="image_3.png" /><figcaption>Figure 3: Left: funnel-shaped variation plot of the cost functions along the continuum in Case C (a); right: variation of the design variables in physical dimensions (b).</figcaption></figure>

tail ratio and wing taper ratio, reduction of wing aspect ratio).

Lastly, the Pareto front associated with the triplet $(f_1, f_2, f_3)$ was calculated with a larger population of 100, and a larger, somewhat arbitrary, budget of 10,000 evaluations, requiring some 77 hrs of computing time. With the addition of the cost function $f_3$, the Pareto front (blue dots on Figure 4) is enhanced. The graphic software gnuplot permits to visualize an interpolated surface (light violet grid). This surface seems to contain the $(f_1, f_2, f_3)$ triplets corresponding to the 3D-extension of the PPF over about two thirds of it, but a gap is visible in Case R and near it. This is the result of the better convergence of the three-objective Pareto front. Nevertheless, the three continua of Nash equilibriums are tangent to the extended PPF, and all three tend to reduce $f_3$. The interpolated surface is clearly not convex, showing several local minima and a maximum inside the domain. This reveals that the cost functions are not all three convex, and this is particularly so near the peak of $f_3$ close to the starting point of Case R. Nevertheless, the green paths all three tend to reduce $f_3$ and remain close to the Pareto front. This verifies that the discrepancy with Pareto-optimality is very small ($O(\varepsilon^2)$ in theory). In Case C, the path points in a valley towards the minimum calculated value of $f_3$.

## 4.4 Numerical experiments using Bayesian optimization

We consider Bayesian optimization (BO) [13] as a replacement for NSGA-II to establish the PPF with less evaluations. BO relies on Gaussian process (GP) metamodels to predict the behavior of the objectives and guide the search towards promising regions of the parameter space. To handle multiple objectives, multi-output GPs are built in combination with sequential infill criteria aiming at finding Pareto optimal points, such as the hypervolume improvement [12] or
<figure><img src="image_4.png" /><figcaption>Figure 4: Two-Objective Primary Pareto Front by NSGA-II, three continua of Nash equilibriums initiated from it, three-objective Pareto front and interpolated grid</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: Three-objective Pareto fronts obtained by NSGA-II (10000 evaluations) and Bayesian Optimization (DoE of 20 evaluations + 100 sequential iterations).</figcaption></figure>
the probability of improvement (PI [15]). The use of the SEGOMOE toolbox [1] with its recent extension to multi-objective problems [14] is illustrated in Figure 5. An initial DoE of 20 points is used to build first metamodels and then 100 iterations are performed using PI based criterion, taking around 1h to complete. Compared to NSGA-II, BO required much less evaluations to obtain a three-objective Pareto front. The front is denser, but also less spread, especially at the extremities. This reflects the increasing difficulty of covering the entire Pareto front as more objectives are added, while searching for specific solutions instead is also advocated e.g., by [2]. With prioritized objectives, finding the PPF with BO is more cost-efficient than with NSGA-II. Further reductions of the number of evaluations in the prioritized approach could be obtained by considering GP metamodels rather than quadratic ones.

# 5 Conclusion and perspectives

The efficacy and cost-efficiency of a prioritized approach have been explored in the context of concurrent engineering, in which several cost functions are to be minimized. The prioritized approach is justified particularly in the following situations:

1. when the cost functions are of different importance, which encourages to give greater importance to preponderant cost functions;

2. when the computational cost of using a specific software to identify Pareto fronts numerically increases rapidly with the number of cost functions accounted for.

Our prioritized approach is a two-step procedure consisting in first computing the Pareto front for the preponderant cost functions only, and to construct from it a continuum of Nash equilibria proved by asymptotic theory to be quasi-Pareto optimal for the whole set. The method has been applied to a test-case of aircraft flight performance optimization involving two primary cost functions and one secondary. Most of our numerical experiments were carried out using the classic NSGA-II for identifying Pareto fronts. Indeed, solutions very close to the Pareto front of the whole set of cost functions were identified by the prioritized algorithm (efficacy), far more economically than by computing the complete Pareto front (cost-efficiency). Our experiments with BO involved a relatively small initial database, and Pareto fronts could be computed much more economically, but resulting in a sparser Pareto front. The prioritized approach resulted in a much denser continuum focused in areas of interest. This comparison would deserve more extensive analysis using larger databases in particular. Additionally, the prioritized approach is promising for engineering cases involving more than three cost functions. Lastly, in Partial Differential Equation-constrained optimization, our ultimate goal, cost efficiency should also be evaluated in the context of parallel computing.
References

[1] N. Bartoli, T. Lefebvre, S. Dubreuil, R. Olivanti, R. Priem, N. Bons, J. R. Martins, and J. Morlier. Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design. *Aerospace Science and Technology*, 90:85-102, 2019.

[2] M. Binois, V. Picheny, P. Taillandier, and A. Habbal. The Kalai-Smorodinsky solution for many-objective Bayesian optimization. *The Journal of Machine Learning Research*, 21(1):5978-6019, 2020.

[3] C. David, S. Delbecq, S. Defoort, P. Schmollgruber, E. Benard, and V. Pommier-Budinger. From FAST to FAST-OAD: An open source framework for rapid overall aircraft design. *IOP Conference Series: Materials Science and Engineering*, 1024(1):012062, jan 2021.

[4] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan. A fast and elitist multiobjective genetic algorithm: NSGA-II. *IEEE Transactions on Evolutionary Computation*, 6(2):182-197, 2002.

[5] J.-A. Désidéri. MGDA Software Platform for Multi-Objective Differentiable Optimization. http://mgda.inria.fr.

[6] J.-A. Désidéri. Platform for prioritized multi-objective optimization by metamodel-assisted Nash games. Research Report 9290, Inria, September 2019. https://hal.inria.fr/hal-02285197.

[7] J.-A. Désidéri. Adaptation by Nash games in gradient-based multi-objective/multi-disciplinary optimization. In *JANO13, Mathematical Control and Numerical Applications, Khouribga, Morocco*, volume 372. Springer Proceedings in Mathematics and Statistics, 2021. https://hal.inria.fr/hal-03430972.

[8] J.-A. Désidéri and R. Duvigneau. Direct and adaptive approaches to multi-objective optimization. Research Report 9291, Inria, September 2019. https://hal.inria.fr/hal-02285899.

[9] J.-A. Désidéri and R. Duvigneau. Prioritized optimization by Nash games: Towards an adaptive multi-objective strategy. In *ESAIM: Proceedings and Surveys*, volume 71, pages 54-63. EDP Sciences, 2021. https://hal.inria.fr/hal-03430912.

[10] J.-A. Désidéri, R. Duvigneau, and A. Habbal. *Computational Intelligence in Aerospace Sciences, V. M. Becerra and M. Vassile Eds.*, volume 244 of *Progress in Astronautics and
*Aeronautics, T. C. Lieuwen Ed.-in-Chief*, chapter Multi-objective design optimization using Nash games. American Institute for Aeronautics and Astronautics Inc., Reston, Virginia, 2014.

[11] J.-A. Désidéri, J. Wintz, N. Bartoli, C. David, and S. Defoort. Combining Pareto optimality with Nash games in multi-objective prioritized optimization of an aircraft flight performance. Research Report RR-9490, Inria - Sophia Antipolis ; Acumes, Oct. 2022.

[12] M. T. Emmerich, K. C. Giannakoglou, and B. Naujoks. Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels. *IEEE Transactions on Evolutionary Computation*, 10(4):421–439, 2006.

[13] R. Garnett. Bayesian optimization. *Cambridge University Press*, 2023.

[14] R. Grapin, Y. Diouane, J. Morlier, N. Bartoli, T. Lefebvre, P. Saves, and J. Bussemaker. Regularized infill criteria for multi-objective Bayesian optimization with application to aircraft design. In *AIAA AVIATION 2022*, 2022.

[15] D. R. Jones. A taxonomy of global optimization methods based on response surfaces. *Journal of global optimization*, 21:345–383, 2001.