# Prioritized optimization by Nash games : towards an adaptive multi-objective strategy

Jean-Antoine Désidéri, Régis Duvigneau

► To cite this version:

Jean-Antoine Désidéri, Régis Duvigneau. Prioritized optimization by Nash games: towards an adaptive multi-objective strategy. ESAIM: Proceedings and Surveys, 2021, 71, pp.54-63. 10.1051/proc/202171106. hal-03430912

HAL Id: hal-03430912

https://inria.hal.science/hal-03430912

Submitted on 16 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PRIORITIZED OPTIMIZATION BY NASH GAMES : TOWARDS AN ADAPTIVE MULTI-OBJECTIVE STRATEGY
## APPLICATION TO AIRCRAFT SIZING FOR OPTIMAL FLIGHT PERFORMANCE

JEAN-ANTOINE DÉSIDÉRI<fnref n="1" /> AND RÉGIS DUVIGNEAU<fnref n="1" />

**Abstract.** This work is part of the development of a two-phase multi-objective differentiable optimization method. The first phase is classical: it corresponds to the optimization of a set of primary cost functions, subject to nonlinear equality constraints, and it yields at least one known Pareto-optimal solution $\mathbf{x}_A^*$. This study focuses on the second phase, which is introduced to permit to reduce another set of cost functions, considered as secondary, by the determination of a continuum of Nash equilibria, $\{\bar{\mathbf{x}}_\varepsilon\} (\varepsilon \ge 0)$, in a way such that: firstly, $\bar{\mathbf{x}}_0 = \mathbf{x}_A^*$ (compatibility), and secondly, for $\varepsilon$ sufficiently small, the Pareto-optimality condition of the primary cost functions remains $O(\varepsilon^2)$, whereas the secondary cost functions are linearly decreasing functions of $\varepsilon$. The theoretical results are recalled and the method is applied numerically to a Super-Sonic Business Jet (SSBJ) sizing problem to optimize the flight performance.

**Résumé.** Ce travail s'inscrit dans le développement d'une méthode d'optimisation différentiable multiobjectif priorisée qui s'applique en deux phases. La première phase est classique: elle correspond à l'optimisation d'un ensemble de fonctions coûts, considérées comme principales, sous contraintes d'égalité non linéaires, et fournit au moins une solution Pareto-optimale $\mathbf{x}_A^*$ connue. On se concentre dans cette étude sur la deuxième phase, qui est introduite pour permettre la réduction d'autres fonctions coûts, dites secondaires, par la détermination d'un continuum d'équilibres de Nash, $\{\bar{\mathbf{x}}_\varepsilon\} (\varepsilon \ge 0)$, de telle sorte que: d'une part $\bar{\mathbf{x}}_0 = \mathbf{x}_A^*$ (compatibilité), et d'autre part, pour $\varepsilon$ suffisamment petit, la condition de Pareto-optimalité des fonctions coûts principales est maintenue à $O(\varepsilon^2)$ près, alors que les fonctions coûts secondaires décroissent linéairement en $\varepsilon$. Les résultats théoriques sont rappelés et la méthode est appliquée numériquement à un problème de dimensionnement d'un avion d'affaires supersonique (SSBJ) pour en optimiser la performance de vol.

## INTRODUCTION

In the engineering environment, the process of numerical optimization is complex due to several coupled elements that we outline in the example of the shape optimization of a civil aircraft wing or configuration. The process is

* multi-disciplinary, since aerodynamic forces, structural and thermal loads, acoustics, in particular, should be assessed and, or optimized;

* multi-objective, since, e.g. in the aerodynamic design alone, it requires the control of cost and constraint functions related to forces (drag and lift), and moments (pitching moment in particular);

<footnote n="1">¹ Université Côte d'Azur, Inria, CNRS, LJAD, 2004 Route des Lucioles, BP 93, 06902 Sophia Antipolis Cedex (France); e-mail: jean-antoine.desideri@inria.fr regis.duvigneau@inria.fr</footnote>
* multi-point, since the design must validated over the entire flight envelope, from the subsonic regime of take-off and landing to the transonic cruise regime, and also for certain extreme conditions in view of certification.

Additionally most of the functions of interest are not directly expressed in terms of the optimization variables, such as geometrical parameters, but are instead functionals of state variables through the partial-differential equations of fluid dynamics, elasticity and of all other physical disciplines to be accounted for, via a complex modeler-mesher-solver coupled sequence [7]. As a result, these functions are evaluated at high cost, and the calculation of the gradient is usually not trivial.

Consequently, many objective functions are potentially significant in the design optimization process and some may reveal, at the stage of an a posteriori validation, more critical than anticipated. Hence, if a preliminary phase of optimization accounted for the most essential design objectives, a second phase may turn out to be necessary to correct the design in view of reducing additional cost functions. The second phase should be conducted in a way such that the optimality of the first design should be the least possible degraded. Our approach is precisely meant to introduce this type of adaptivity in the multi-objective optimization process.

In Section 1, the problematics is stated in mathematical terms, a theoretical result formulated, and a numerical method proposed. In Section 2 the method is applied to a Super-Sonic Business Jet (SSBJ) sizing problem devised for optimal flight performance.

# 1. THEORETICAL BACKGROUND

## 1.1. Mathematical problem definition

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$ ($1 \le m < M$) of which are thought of preponderant, or primary, and the remaining ones ($j = m + 1, \dots, M$), secondary, and all are subject to $K$ equality constraints, $c_k(\mathbf{x}) = 0 (k = 1, \dots, K)$. Objective and cost 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$ 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 $\{\bar{\mathbf{x}}_\varepsilon\}$ parameterized by $\varepsilon$ ($\varepsilon \ge 0$) in way such that: firstly, $\bar{x}_0 = \mathbf{x}_A^*$ (compatibility), and secondly, for $\varepsilon$ sufficiently small, the Pareto-optimality condition of the primary cost functions remains $O(\varepsilon^2)$, whereas the secondary cost functions are decreasing functions of $\varepsilon$ (linearly, or faster).

## 1.2. Preliminary calculations

The preliminary calculations have been defined in full in [4]. 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 * on any symbol denotes throughout an evaluation at $\mathbf{x} = \mathbf{x}_A^*$, $\nabla$ is the gradient operator w.r.t. $\mathbf{x}$, $\{\alpha_j^*\} (j = 1, \dots, m)$ are the known coefficients of a convex combination, $\mathbf{P}$ is the projection operator
onto the linear subspace orthogonal to all constraint gradients, {$\nabla c_k^*$} ($k = 1, \dots, K$). 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 \quad (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 (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 [4, 5]:

$$\mathbf{H}'_A = \mathbf{P}\mathbf{H}_A^{+,*}\mathbf{P} = \Omega\mathbf{H}\Omega^t \quad (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 eigenvectors, column-vectors of matrix $\mathbf{\Omega}$, are conventionally arranged in a special way: 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 eigenmodes 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}) \quad (4)$$

where $\mathbf{u} \in \mathbb{R}^{n-p}$, $\mathbf{v} \in \mathbb{R}^p$, and $p$ is an integer free parameter to be chosen such that $1 \le p < n-K$. 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} \quad (5)$$

and the following scaled logarithmic gradients are associated with the secondary cost functions:

$$\mathbf{g}_j = \mathbf{S}^{-\frac{1}{2}} \frac{\nabla_{\mathbf{v}} f_j^*}{f_j^*} = \mathbf{S}^{-\frac{1}{2}} \frac{\mathbf{\Omega}_{\mathbf{v}}^t \nabla f_j^*}{f_j^*} \quad (j = m + 1, \dots, M) \quad (6)$$

where the scaling matrix $\mathbf{S} = (\nabla_{\mathbf{v}\mathbf{v}}^2 f_A^+)^*$ is the lower $p \times p$ diagonal block of matrix $\mathcal{H}$ of (3), positive-definite by convexity fix. Hence these gradients reflect the sensitivities of the secondary cost functions when $\mathbf{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) [4], one identifies the minimum-norm element $\mathbf{\omega}_B^*$ in the convex hull of the gradients $\mathbf{g}_j$:

$$\mathbf{\omega}_B^* = \sum_{j=m+1}^{M} \alpha_j^* \mathbf{g}_j \quad (7)$$

and the coefficients {$\alpha_j^*$} of this convex combination are used to define the secondary steering function:

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

and this completes the preliminaries.
## 1.3. Nash game formulation and theoretical results

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

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

For every fixed value of $\varepsilon$, two virtual players, Player A and Player B, respectively controlling the sub-vectors $\mathbf{u}$ and $\mathbf{v}$ in (4), are involved in a Nash game defined by the following 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}$,

whereas:

- 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 [4]:

* The Nash equilibrium point $\bar{\mathbf{x}}_\varepsilon$ exists for all sufficiently small $\varepsilon$, and

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

an essential property referred to as compatibility or consistency.

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

$$f_A^+ = 1 + O(\varepsilon^2) \quad (11)$$

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

* If $\phi_B(\varepsilon) = f_B(\bar{\mathbf{x}}_\varepsilon)$, $\phi'_B(0) = -\sigma_B$, where $\sigma_B = \|\omega_B^*\|^2 \ge 0$. Hence, the secondary steering function $f_B(\mathbf{x})$ decreases linearly as $\varepsilon$ 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.

## 1.4. Numerical implementation

The present method can be applied formally to simple cases defined analytically. In the more general case in which the formal derivation is not at hand, we propose to resort to models of substitution, and formulate the Nash game in terms of these meta-models. Precisely, quadratic models can be constructed once for all for the two steering functions based on local differential elements at $\mathbf{x} = \mathbf{x}_A^*$ (typically approximated by central differencing), and function values over a database to account on second derivatives more globally. For constraints, quadratic meta-models can also be constructed, but it is recommended to upgrade these models, as the determination of the continuum of Nash equilibria proceeds, because constraints participate critically to the definition of Pareto-optimality, and this notion cannot be altered legitimately. Such techniques have been implemented over the software platform http://mgda.inria.fr according to procedures fully defined in [6].

# 2. APPLICATION TO A SSBJ SIZING PROBLEM FOR OPTIMAL FLIGHT PERFORMANCE

The test-case is inspired from Multi-Disciplinary Optimization (MDO) problematics raised by Dassault Aviation within the ANR (French Research Agency) "Optimisation Multidisciplinaire" Project (2006-2009) documented in [2, 3]. Dassault Aviation had provided partners with a software permitting to exercise numerical procedures for multidisciplinary or multi-objective problems in flight-mechanics design [7]. This software has been used again here for purpose of illustration of our method for prioritized multi-objective optimization.

The software is made of several procedures that calculate certain performance quantities according to classical laws of flight-mechanics, the Breguet laws in particular, as functions of 15 design variables {$X_i$} ($i = 1, \dots, 15$),
all subject to prescribed interval constraints $X_{i,min} \le X_i \le X_{i,max}$. For the design of a SuperSonic Business Jet (SSBJ) these variables and their bounds are described in Table 1.

<table><caption>TABLE 1. Physical design variables in the flight-mechanics test-case and their specified bounds</caption><thead><tr><th>symbol<br>(X<sub>i</sub>)</th><th>significance<br>(unit)</th><th>lower bound<br>X<sub>i,min</sub></th><th>upper bound<br>X<sub>i,max</sub></th></tr></thead><tbody><tr><td>z</td><td>cruise altitude (m)</td><td>8000</td><td>18500</td></tr><tr><td>xmach</td><td>cruise Mach number</td><td>1.6</td><td>2.0</td></tr><tr><td>S</td><td>wing reference surface (m<sup>2</sup>)</td><td>100</td><td>200</td></tr><tr><td>phi0w</td><td>wing leading-edge sweep angle (°)</td><td>40</td><td>70</td></tr><tr><td>phi100w</td><td>wing trailing-edge sweep angle (°)</td><td>-10</td><td>20</td></tr><tr><td>xlw</td><td>wing taper ratio</td><td>0.05</td><td>0.50</td></tr><tr><td>t_cw</td><td>wing relative thickness</td><td>0.04</td><td>0.08</td></tr><tr><td>phi0t</td><td>vertical-tail leading-edge sweep angle (°)</td><td>40</td><td>70</td></tr><tr><td>phi100t</td><td>vertical-tail trailing-edge sweep angle (°)</td><td>0</td><td>10</td></tr><tr><td>xlt</td><td>vertical-tail taper ratio</td><td>0.05</td><td>0.50</td></tr><tr><td>t_ct</td><td>vertical-tail relative thickness</td><td>0.05</td><td>0.08</td></tr><tr><td>dfus</td><td>fuselage diameter (m)</td><td>2.0</td><td>2.5</td></tr><tr><td>wfuel</td><td>fuel mass (kg)</td><td>15,000</td><td>40,000</td></tr><tr><td>alpha</td><td>landing maximum angle of attack (°)</td><td>10</td><td>15</td></tr><tr><td>xfac</td><td>mlw/tow, landing to take-off mass ratio</td><td>0.85</td><td>0.95</td></tr></tbody></table>

In order to respect the prescribed interval constraints, we have used, for each design variable $X_i$, a transformation of the following type

$$X_i = \frac{X_{i, \max} + X_{i, \min}}{2} + \frac{X_{i, \max} - X_{i, \min}}{2} \sin x_i \quad (12)$$

and used the dimensionless vector $\mathbf{x} = \{x_i\}$ as the optimization variable.

The software permits to compute in particular

* mass at take-off, $M$ (to be minimized),

* range, $R$ (to be maximized),

* approach speed, $V$ (to be minimized),

* take-off distance, $D$ (to be minimized, or maintained below a specified upper bound).

In the present test-case, the take-off distance $D$ was subject to the following inequality constraint:

$$D \le 1828 \text{ m.} \quad (13)$$

This led us to introduce an additional variable, $x_{16}$, as a slack variable, and to replace the above inequality by the following equality constraint:

$$c_1(\mathbf{x}) = \frac{D}{1828} - 1 + x_{16}^2 = 0 \quad (14)$$

thus increasing by one the dimension of the vector $\mathbf{x}$.
Three cost functions were defined in relation to mass, range and approach speed:

$$f_1(\mathbf{x}) = \exp\left[\gamma\left(\frac{M}{M^*} - 1\right)\right] \qquad (15)$$

$$f_2(\mathbf{x}) = \exp\left[\gamma\left(\frac{R^*}{R} - 1\right)\right] \qquad (16)$$

$$f_3(\mathbf{x}) = \exp\left[\gamma\left(\frac{V}{V^*} - 1\right)\right] \qquad (17)$$

where starred quantities correspond to initial values at $\mathbf{x} = \mathbf{x}_A^*$ and are case-dependent as well as the positive numerical parameter $\gamma$. In this setting, the parameter $\gamma$ appears as a scaling factor in the expression of the initial gradients of the cost functions, since it permits to render their magnitudes close to 1. The form given to (15) and (17) reflects the objective to minimize mass and approach speed, and to (16) to maximize range. The three functions $f_1$, $f_2$, $f_3$ are uniformly positive and equal to 1 initially:

$$f_1^* = f_2^* = f_3^* = 1. \qquad (18)$$

## 2.1. Problematics, Pareto front and prioritized optimization.

The cost functions $f_1$ and $f_2$ representing mass and range have been prioritized, and $f_3$ representing approach speed considered as a secondary cost function, all three functions being subject to the constraint $c_1(\mathbf{x}) = 0$ related to take-off distance.

In a first experiment, not involving the slack-variable $x_{16}$, the mass-range Pareto front associated with the pair $(f_1, f_2)$ subject to the equality constraint $c_1(\mathbf{x}) = 0$ handled by penalization was determined approximately by applying the *Pareto Archived Evolution Strategy* (PAES) [1]. This task revealed more laborious than we had anticipated due to the small scales affecting the design variables. After some 100,000 evaluations of the functions, eliminating all the points violating the constraint (13) and all the dominated points, resulted in the approximate Pareto front represented by the symbols in black of Figure 1.

Five points A, B, C, D and E associated with mass-range Pareto-optimal solutions subject to the distance constraint were selected from the right-most point of the approximate discrete Pareto front (largest mass and range), to the left-most point (smallest mass and range). Each of these points was used to conduct an independent numerical experiment, in which the continuum of Nash equilibria originating from it has been calculated with the objective to diminish the approach speed. In this way, we obtained five different paths initiated from the front and represented on Figure 1 in different colors.

We refer to [6] for the precise description of this experiment. Here, the discussion is limited to the numerical results.

Let us first examine the most standard case corresponding to the green path initiated at B. In this experiment, the parameter $\gamma$ was set to 10 to better visualize the variations of the functions. Evidently, it is confirmed that the path is tangent to the Pareto front and constitutes a satisfactory approximation of it over a significant segment.

A great part of the path remains relatively close to the front. For example, the point B' at which $\varepsilon/\varepsilon_{max} = 0.8$, where $\varepsilon_{max}$ is the theoretical upper bound for convexity defined in [6] is indicated by the symbol ●. Thus, the portion of the path beyond this point only represents the 20% tail of the continuum.

The variation with $\varepsilon$ of the three cost functions, and two steering functions $f_A^+$ and $\tilde{f}_B$ (meta-model of $f_B$) is given on Figure 2 on which the values at point B' are indicated by the symbol ●.

Several observations can be made from this plot:

* The general pattern of the function plot is horn-shaped; here bounded from above by the increasing primary steering function and from below by the decreasing secondary steering function.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Discrete mass-range Pareto front and five continua of Nash equilibria (abscissa: mass in kg; ordinate: -range in m)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Variation with $\varepsilon$ of the cost functions $f_1, f_2, f_3$ and the two steering functions $f_A^+$ and $\tilde{f}_B$ in Case B. The values at point B' are indicated by the symbol •.</figcaption></figure>

* The primary steering function $f_A^+$ has indeed a null initial derivative. The two primary cost functions vary inversely: $f_1$ diminishes, and $f_2$ increases. For example, at B':

$$f_1 \approx 0.866, f_2 \approx 1.307 \Rightarrow M/M^* \approx 0.986, R^*/R \approx 1.027. \quad (19)$$

Hence mass was reduced of 1.4%, a gain, and range diminished 2.7%, a loss. The function $f_A^+$ has a more rapid increase due to the convexity-fix term.

* The initial derivative of the secondary steering function $\tilde{f}_B$ is well given by the theoretical prediction indicated by a dashed tangent line of slope $-\sigma_B$. The function $f_3$, calculated a posteriori, follows accurately its meta-model $\tilde{f}_B$. At B':

$$f_3 \approx 0.571 \Rightarrow V/V^* \approx 0.944. \quad (20)$$

Hence the approach speed was reduced of some 5.6%. This is confirmed by the plot of V as a function of mass along the continuum given by Figure 3, over which the value at B' is indicated by the symbol •.

Correspondingly, the variations with $\varepsilon$ of the optimization variables $\{x_i\}$ ($i = 1, \dots, 16$) are given by Figure 4 over which the values at B' are indicated by the symbols •. Only the first four of these variables have a definite variation, the other ones remain nearly constant. Also note that the angular point in the plot of $x_3(\varepsilon)$ does not reflect a derivative discontinuity, but a convention adopted in the graphics representation [6].

In Figure 5, we present the variation with $\varepsilon$ of the constraint $c_1$ related to the take-off distance. Over a broad portion of the interval, the constraint is satisfied with very high accuracy thanks to the local constraint meta-model, upgraded prior to the calculation of each new Nash equilibrium point. For example, at point B', the value is $c_1 \approx -3.7 \times 10^{-5}$, while TOL = $10^{-4}$. Beyond this point, the constraint value ceases to be negligible rapidly. These values are firstly negative which implies that an unnecessarily severe strict inequality is satisfied, that is, some significant distance away from the Pareto front, the approach speed is overly reduced at the cost of an unduly mass increase and range reduction. Further away, when approaching $\varepsilon = \varepsilon_{max}$, an instability is triggered. As a result, the tail of the continuum is meaningless. Perhaps a more sophisticated constraint meta-model could permit to prolongate the continuum further. This point is currently being examined.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Approach speed as a function of mass along the continuum. The value at point B' is indicated by the symbol •.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Variation with $\epsilon$ of the optimization variables $\{x_i\}$ ($i = 1, \dots, 16$). The values at B' are indicated by the symbols •.</figcaption></figure>

We finish this description of Case B by some information relative to the number $\lambda$ of outer iterations necessary to coordinate the two subvectors **u** and **v** at a given $\epsilon$ for the specified tolerance $TOL = 10^{-4}$, and the number $\mu$ of Newton's method inner iterations to solve for **u** for fixed **v** with accuracy tolerance of $TOL/100 = 10^{-6}$. Note that over the hundred $\epsilon$ steps, 94 were successful to determine the Nash equilibrium. In all these cases, $\mu$ never exceeded 3, while the termination value of $\lambda$ is given by Figure 6 in terms of $\epsilon$. Over more than 60% of the continuum, $\lambda$ remains at most equal to 2. Then when approaching the limit of convexity it increases, and finally the process diverges, or more rigorously speaking, is interrupted. These results reflect some arbitrariness related to the setting of numerical parameters; they would change with a different setting of the accuracy tolerances, or with the specification of a different maximum allowable number of iterations in the inner or outer loop.

We now discuss the results related to the continuum originating from point A, at the extreme right of Figure 1 corresponding to the case of largest mass and range. These results are similar to those of Case B except that the portion of the continuum close to the Pareto front extends somewhat less. Again, a lighter aircraft could permit a smaller approach speed, but at the cost of a smaller range.

Thirdly, the continuum initiated at point C remains very close to the discrete Pareto front over a large segment, over which it is an excellent approximation to it.

Fourthly, the continuum initiated at point D is very small and obtained with difficulty using the very large value of 1000 for the parameter $\gamma$. In this area, there seems to occur an inversion of trend, more visible from the continuum initiated at point E, using $\gamma = 40$, and a very low upper bound on the condition number ($\kappa = 1.5$) in which results a large convexity-fix constant. The continuum there makes a loop, indicating a form of instability.

Evidently, the results in Cases D and E do not bring improved practical designs.

In conclusion of this case study, we emphasize the following points:

* the test-case is of actual technical interest;

* it involves 15 design variables subject to interval constraints that were handled by a change of variables, and one functional constraint handled by the introduction of a slack variable; many of these variables however remained almost constant along the continuum;

* the scales in the design variables were found very small and made the establishment of the Pareto front by the evolutionary strategy laborious.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Variation with $\varepsilon$ of the landing distance constraint along the continuum. The value at B' is indicated by the symbol •.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Number $\lambda$ of outer iterations necessary to coordinate the sub-vectors **u** and **v** at a given $\varepsilon$ as a function of $\varepsilon$. (The process was interrupted whenever the convergence test still was not met after the 15th iteration.)</figcaption></figure>

Nevertheless, we were able to identify a continuum of Nash equilibria originating from a large portion of the Pareto front. This continuum complies with the theoretical findings, in particular w.r.t. tangency to the Pareto front. Along the continuum, the cost functions exhibit a horn-shaped pattern indicating a primary steering function that increases with the continuation parameter $\varepsilon$ and a secondary steering function that decreases, initially linearly, in compliance with the theoretical results. The constraint is very well satisfied over a large part of the interval $[0, \varepsilon_{max}]$, where $\varepsilon_{max}$ is the known theoretical limit of convexity of the meta-model-based problem. Some extreme cases required an unusual adjustment of the parameter $\kappa$ (which controls the added convexity term), and $\gamma$ which controls the magnitudes of initial gradients.

In these experiments, we have determined a priori an approximate but entire Pareto front by PAES. This was for purpose of numerical illustration and verification of the theoretical relationship between the front and the continuum. However, only a few Pareto-optimal points were actually necessary to initiate the continuums, and the first phase of the experiment could have been limited to the determination of these few points. Also note that this determination could have been conducted by other means, perhaps more efficiently, and the Pareto front roughly sketched by the continuums themselves alone. This would have been an alternative to the evolutionary strategy.

Lastly, to establish one such continuum typically required about 45 sec of run time on a standard laptop for program assembly and compilation, execution (mostly in the 15 phases of preparation of the Nash games), and the production of plots using the gnuplot software. In this test-case that involved 16 variables, the database contained 1920 points about each $x_A^*$. Consequently, most of the computational work was devoted to the establishment of the meta-models and the evaluation of the primary agglomerated function $f_A$, and relatively less to the calculation of the Nash equilibria per se. For example, in Case B, 94 equilibrium points were calculated. The numerical experiment required a total of 2048 calls to the evaluation procedures of the primary and secondary cost functions, and 5055 calls to the procedure of constraint evaluation.

The authors wish to acknowledge the very fruitful technical discussions with M. Ravachol about the supersonic business jet problematics, and are indebted to Dassault Aviation that provided the flight mechanics software.
REFERENCES

[1] J. Knowles and D. Corne. The Pareto Archived Evolution Strategy: A New Baseline Algorithm for Pareto Multiobjective Optimisation. In *Proceedings of the 1999 Congress on Evolutionary Computation, CEC 1999*, January 1999.

[2] R. Filomeno Coelho and P. Breitkopf, editor. *Optimisation multidisciplinaire en mécanique 1 - démarche de conception, stratégies collaboratives et concourantes, multiniveaux de modèles et de paramètres*. Mécanique et Ingénierie des Matériaux. Hermes Lavoisier, Paris, 2009. www.hermes-science.com www.lavoisier.fr.

[3] R. Filomeno Coelho and P. Breitkopf, editor. *Optimisation multidisciplinaire en mécanique 2 - réduction de modèles, robustesse, fiabilité, réalisations logicielles*. Mécanique et Ingénierie des Matériaux. Hermes Lavoisier, Paris, 2009.

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

[5] J.-A. Désidéri and 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.

[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] M. Ravachol. *Multidisciplinary Design Optimization in Computational Mechanics*, chapter Multilevel Multidisciplinary Optimization in Airplane Design. Wiley-ISTE, May 2010.