Gradient-based Pareto front approximation applied to turbomachinery shape optimization

Multi-objective optimization has been rising in popularity, especially within an industrial environment, where several cost functions often need to be considered during the design phase. Traditional gradient-free approaches, such as evolutionary algorithms, can be employed to compute a front of equally suitable solutions a designer can choose from, with a high computational cost though, particularly in high-dimensional design spaces. In this paper, a gradient-based algorithm is developed for efficiently tracing the Pareto front in bi-objective aerodynamic shape optimization problems, where an adjoint method is used for the computation of the objective functions’ gradients with respect to the design variables. After obtaining a starting point on the front, a prediction–correction approach is employed to compute new Pareto points. Satisfying the Karush–Kuhn–Tucker conditions provides a prediction for the next point, which is, then, corrected by solving a minimum distance problem. The prediction step, though, requires the costly computation of the Hessian matrix. This is avoided here using the BFGS (Nocedal and Wright in Numerical optimization, Springer, New York, 2006) technique. The proposed method is first demonstrated in a less expensive lift/drag optimization of an isolated airfoil and then applied to the bi-objective optimization of a 3D compressor stator. The extension of the proposed method to cases with more than two objectives is straightforward, on condition that an algorithm is found to coordinate the way the Pareto front is swept.


Introduction
Optimization in aerodynamics and engineering, in general, may involve several contradictory objectives. These multiobjective optimization (MOO) problems can be solved using stochastic approaches, such as evolutionary algorithms (EAs), which are able to capture the front of non-dominated solutions with a single run. EAs treat the objective function evaluation tool as a "black box" and, thus, no additional software development is required. Provided that they are run for a sufficiently long time, EAs are able to find the global optima, by computing convex, non-convex and discontinuous Pareto fronts [1,2]. However, the main disadvantage associated with EAs is that their computational cost scales with the number of design variables [3]. As a remedy to this problem, metamodel-assisted EAs (i.e. EAs regularly implementing low-cost surrogate evaluation models within the evolution) and/or dimensionality reduction techniques using, for instance, principal component analysis [4] can be applied.
On the other hand, gradient-based approaches deal with MOO problems by performing a series of single-objective optimization (SOO) runs. The most popular gradient-based strategy is the weighted-sum approach, in which a convex combination of the individual objectives is minimized. Differently weighting the objectives leads to different nondominated solutions and, thus, the process is repeated until the Pareto front is adequately populated. However, a uniform spread of weights rarely results in a good distribution of points along the Pareto front. If the Pareto front is nonconvex, none of the combinations of weights may compute internal front points.
To approximate the Pareto front with a number of welldistributed optimal points, several approaches have been proposed. The Normal Boundary Intersection (NBI) method, where the Pareto front is sought by projecting points lying on the convex hull of local minima of each individual objective, is presented in [5]. A uniformly spaced grid of points on the convex hull results in a good distribution of Pareto optimal solutions. Although this method is applicable to high objective space dimensions, it is not guaranteed that the whole Pareto front is found, whereas special treatment is required at the non-convex parts of the front.
These issues of the NBI can be addressed by the adaptive weighted-sum method [6], where, instead of prescribing a-priori weights to populate the Pareto front, weights are dynamically modified. A small number of SOO problems are initially solved for given combinations of weights and then the unexplored Pareto regions are populated by specifying additional inequality constraints. The extra computational cost as well as the dependency on some user-defined parameters are the main weaknesses of this method. In addition, extension to problems with more than two objectives is not straightforward. An improvement to the NBI method for high objective space dimensions is proposed in [7], successively populating the inner parts of bounded Pareto fronts. For such cases, the boundary of the Pareto front can be computed using objective reduction techniques [8].
A different scalarization technique using steepest descent is described in [9], where convex and non-convex well-distributed Pareto fronts are generated. This is also achieved in [10] by augmenting the Karush-Kuhn-Tucker (KKT) conditions with additional objective-space-based distance constraints and, then, using Newton's method to solve the coupled system of equations, which drive the set of points towards Pareto optimality and uniform spacing simultaneously. The required (exact) Hessian matrix is derived by a hyper-dual numbers implementation, at very high cost, though, for practical applications.
Target or reference point methods have been used in optimal flow control [11] as a robust way of moving from one Pareto point to another. They require the successive solution of many SOO problems, which is computationally intensive and, as a remedy to this, reduced-order approaches based on Proper Orthogonal Decomposition have been proposed [12].
An alternative to scalarization techniques are the curve continuation strategies which aim at moving from one Pareto point to a neighbouring one by taking into account the local front curvature. Such an approach is introduced in [13] in the form of a prediction-correction scheme. A constrained minimization problem is defined, where equality constraints are imposed to all but one objective function which is to be minimized. Formulating the KKT conditions and satisfying the requirements of the implicit function theorem, a linear system is derived, the solution of which can be used to predict a new Pareto point. Due to the non-linearity of the problem, an iterative correction step is required to bring the point onto the Pareto front.
The corrector can be any constrained optimization solver, such as sequential quadratic programming (SQP). Although this approach has the advantage of controlling the distance between two Pareto points, this depends on the user-defined value of the constraints, a bad choice of which can lead to the SQP problem being unsolvable. Moreover, the need for the Hessian matrix can be a limiting factor regarding performance. In [14], a number of alternatives to computing the Hessian for such a curve continuation approach are presented. These include Hessian approximations using BFGS and the computation of Hessian-vector products, instead of the Hessian itself, based on the truncated Newton method [15].
In this work, a different gradient-based prediction-correction algorithm is developed for efficiently moving on the Pareto front. The prediction phase is inspired by the corresponding step used in [13], while key differences exist in the correction phase. Instead of solving a constrained SOO problem at each correction step, an unconstrained formulation is proposed, where the Pareto front is traced by minimizing the distance (in the objective space) from properly defined target points. The distance function is treated as the objective in the formulation of the KKT stationarity condition, computing the minimum distance solution. The required gradients are obtained using either discrete or continuous adjoint; in the results, Sect. 3, the first case is handled using continuous [16] whereas the second one uses the discrete [17] adjoint. BFGS is selected as the optimization method in the correction step, providing thus a Hessian approximation to be utilized in the subsequent prediction step. After demonstrating the proposed method on an isolated airfoil, its industrial applicability is indicated using a 3D turbomachinery component which is presented in detail.

Method formulation
A MOO problem is defined as are t he objective and constraint functions, B ∶= {b ∈ ℝ N | c(b) = 0, g(b) ≥ 0} denotes the feasible solutions' set, N is the number of design variables and M > 1 is the number of objectives. Since a single b * ∈ B cannot minimize all objective functions f i at the same time, the optimality definition is expanded in the Pareto sense. For b 1 A solution b * is called Pareto optimal if it is non-dominated [5]. Then, the Pareto set P s and Pareto front P f are defined as P s = {b * ∈ B | b * is Pareto optimal} and P f = F(P s ) , respectively.
As discussed in Sect. 1, various methods exist to solve problem (1). In this work, a target point approach is enhanced using a prediction technique. After obtaining an initial point on the Pareto front (upon completion of the so-called "Goto-Pareto" step), an infeasible target point is automatically computed in the objective space and the Euclidean distance of the objective vector to that target vector is minimized, resulting in a new Pareto point. This process is repeated until the Pareto front is adequately populated. Moving between neighbouring Pareto points is accelerated using a prediction-correction technique, in which the prediction is driven by approximate Hessian information.

The "Go-to-Pareto" step
To initialize the method, a first point b * 0 that belongs to the Pareto front must be found. This is achieved by minimizing a weighted-sum of the M objectives subject to the constraints of problem (1), where w j ≥ 0 ∀ j ∈ {1, … , M} . Populating the Pareto front by performing a series of "Go-to-Pareto" steps (3), using different weight ( w ) combinations, is computationally expensive and nonapplicable to non-convex problems; thus, an alternative approach is proposed herein.

The "Move-on-Pareto" steps
Moving along the Pareto front can be achieved by means of appropriately selected target points. This technique can be found in the literature with a wide variety of slightly different formulations [11,18]. At each step, a target point In what follows, points marked with ( * ) belong to the sought Pareto front. T i is an infeasible point, according to the Pareto optimality condition (2). Then, a new (ith) Pareto point is obtained by solving the distance minimization problem (with non-dimensional values for the objective components f j ): subject to the constraints of problem (1). Thus, the solution to problem (5) is geometrically defined as the projection of T i onto the Pareto front. By updating the target position based on Pareto points already found, additional points on the Pareto front can be computed one after the other. In general, this "Move-on-Pareto" technique is not restricted to low objective space dimensions, and it amounts to finding an algorithm for redefining the targets in such a way that the Pareto front is adequately populated.
The simplicity and robustness of the proposed approach are the key ingredients which support its implementation within an industrial workflow. However, the method's efficiency cannot be neglected, since solving problem (5) multiple times is cost intensive. In this paper, this is tackled by adding a prediction step (for the next Pareto point), which provides a good initial guess and accelerates the convergence significantly.
The prediction step is inspired by the Pareto curve continuation method presented in [13]. The main difference lies in the fact that the implicit function theorem is applied to problem (5) instead of (3). Further consideration of constraints is omitted, since only unconstrained optimization problems are addressed in this work. However, there is no gap in the method, since the generic formulation (including constraints) can be derived following [13] and solved using, for instance, the SQP method [14].
Let b * i be the solution to problem (5) for the target point T i as defined in Eq. (4). This should satisfy the KKT conditions: In order for a point to remain on the Pareto front, the total derivative of G(b,f k ) w.r.t. the target (i.e. f k ) should be zero: After solving the linear system of Eq.

the Pareto front is given by a first-order Taylor expansion
It should be noted that, in Eq. (7), the Hessian matrices ∇ 2 b f j are required and, thus, they should be computed or approximated.
The new point b * i+1 is obtained by solving the distance minimization problem (5) from an updated target point T i+1 while using b pred i+1 as an initial guess. This step can be considered as the correction to the previous prediction. Mathematically speaking, a T i+1 consistent with the above prediction is obtained by displacing T i using the perturbation vector [0 ⋯f k ⋯ 0] T . However, for small values of f k , displacing any point that lies along the line defined by T i and F(b * i ) towards that direction gives a different target point, which, when used in problem (5), should result in a similar Pareto solution. In practice, this was verified during the preliminary testing of the method and, thus, defining every new target point T i+1 w.r.t. the previously computed Pareto point F(b * i ) , as in Eq. (4), was found to be the best option, especially when combined with a constant f k . In [11], an alternative way of defining target points is proposed using a tangent approximation, to deal with Pareto fronts of highly varying steepness.

Bi-objective optimization algorithm
In this work, the proposed method is applied only to problems with two objectives ( M = 2 ), as illustrated in Fig. 1, using Algorithm 1. A BFGS-based method was developed in order to solve the SOO problems (3) (i.e. to obtain an initial Pareto point) and (5) and at the same time provide an approximation to the Hessian matrices needed in Eq. (7). The BFGS implementation includes a line search algorithm which aims at accelerating the optimization by satisfying at each step the Armijo condition [19]. This condition assures that the chosen step length sufficiently decreases the function to be minimized. Thus, finding an appropriate step length comes at a cost of up to two additional function evaluations, while reducing the required number of the (more expensive) BFGS steps. The gradients needed to drive the optimization method are obtained by adjoint solvers.  (4): Solve linear system (7) → δb To quantify the method's efficiency improvement due to the addition of prediction steps, in Sect. 3, the first "Moveon-Pareto" step is always performed using exclusively the iterative correction step, while omitting the prediction step used in any subsequent Pareto point search, as shown in Fig. 1 Proposed method in the bi-objective space: a Computation of Pareto point b * 1 by solving the distance minimization problem (5) starting from b * 0 (without prediction); b definition of new target point T i+1 and acceleration of the method using a prediction step after solving the linear system of Eq. (7) Fig. 1a. A comparison is made between the cost of computing b * 1 and the rest of the Pareto points. The second condition in line 6 of Algorithm 1 is employed to ensure that, in case the chosen number of steps n s is too big, the algorithm will terminate when the distance between two subsequent Pareto points becomes smaller than . Herein, = 0.1 | |f2 | | is selected.

Shape optimization of an isolated airfoil
Before moving on to the 3D turbomachinery application, the proposed method is demonstrated in the lift/drag optimization of the isolated NACA0012 airfoil. The initial geometry along with the applied parameterization using two Bezier curves [20] is plotted in Fig. 2. The symmetric baseline shape is parameterized with 14 control points (CPs), 6 out of which remain fixed and the rest are allowed to move in the y direction, resulting in 8 "free" design variables. In the absence of any other constraint (such as a volume constraint), reasonable bounds are selected for these parameters to avoid non-realistic airfoil shapes during the optimization. Then, an unstructured C-type mesh consisting of ∼ 75k hybrid elements is generated. At each design iteration, the mesh is updated to match the varying airfoil geometry using an RBF mesh morpher [21].
Regarding the flow conditions, the freestream flow angle, Mach and Reynolds numbers are: ∞ = 3 • , M ∞ = 0.15 and Re ∞ = 3 × 10 6 , respectively. Flow simulations are performed using a CFD code for unstructured grids [22], running on graphics processing units (GPUs). The steady-state Reynolds-averaged Navier-Stokes (RANS) equations for compressible flow are solved using a vertex-centered finitevolume discretization scheme; turbulence is resolved with the one-equation Spalart-Allmaras (SA) model [23]. Lift and drag gradients w.r.t. the design variables are computed using the continuous adjoint method [16].
The MOO problem is defined as where c L and c D are the lift and drag coefficients of the airfoil. Applying Algorithm 1 with n s = 10 and f 2 = −0.1 , the results plotted in Fig. 3 are obtained. Here, the aim is not to generate the whole Pareto front (not both corner points were sought), but rather to verify the applicability of the adjointbased approach. Additionally, in order to perform the MOO within an acceptable computational cost, the convergence criteria of the SOO problems solved in the correction steps are reasonably relaxed. The cost is measured in terms of the equivalent flow solutions (EFS; one EFS stands for the cost of solving the flow equations using the CFD solver of the case) required for the primal and adjoint simulations, so the "Go-to-Pareto" step involved 30 EFS, the first "Move-on-Pareto" step (without prediction) needed 18 EFS and each subsequent step (using prediction-correction) required on average ∼ 12 EFS, resulting in a total computational cost of 159 EFS for obtaining 11 Pareto points. The cost breakdown demonstrates the benefit of the proposed method in terms of efficiency.
To verify the results, a metamodel-assisted EA-based (MAEA) optimization software (EASY [24], developed by the NTUA) was employed to generate the Pareto front. As shown in Fig. 3, the adjoint-based solution set is in good agreement with the MAEA one (only part of which is shown here), both of which are sufficiently approximating the front. More Pareto points are computed by the MAEA optimization  Fig. 2 Parameterization of baseline airfoil geometry using Bezier curves. The "fixed" CPs are not displaced during the optimization, whereas the "free" ones can move vertically (one design variable per "free" CP)  Fig. 3 Solution of the lift/drag problem for an isolated airfoil: Pareto front approximation by 11 optimal points obtained using the proposed method and comparison with MAEA results. The dashed line stands for the "Go-to-Pareto" step at a slightly higher total cost (175 EFS) though. Due to the relaxed termination criteria, the set of non-dominated solutions obtained by the adjoint method is an approximation to the (exact) Pareto front; this is, however, computed faster than the MAEA front, even for a rather small value of design variables ( N = 8).

Shape optimization of a compressor stator
Herein, the proposed method is applied to the bi-objective optimization of a 3D compressor stator.

Description and CFD setup
This subsection presents the case setup with focus on the steps that have to be considered to perform the CFD simulation. A similar setup has been previously used in [25]. The initial stator geometry is available online [26]. The goal is to optimize this baseline stator (BS) w.r.t. two criteria, which leads to a MOO problem. The first objective to be minimized is the total pressure losses between the inlet and exit, defined by the loss coefficient: where p and p t denote the mass-averaged static and total pressure computed at the inlet (I) and exit (E) of the CFD domain. The second objective to be minimized is the flow angle deviation from the axial direction at the stator exit: where u, v, w denote the Cartesian velocity components and ṁ the mass flow rate. Only the design point operating conditions are considered, with a uniform inlet whirl angle profile of 42 • . These two objectives are contradicting since the higher the turning-which is required here to achieve a more axial outflow-the higher the losses. Moreover, several geometric constraints (such as minimum leading edge (LE) and trailing edge (TE) radii, constant axial chord length, etc.) are taken into account during the optimization process. Their detailed definition, as well as information about the boundary conditions of the case, can be found in [26].
The BS geometry is parameterized using the Rolls-Royce in-house blade design tool [27] and is constructed by 21 individual 2D profile sections radially stacked from hub to tip. A curve in the radial direction is defined by axial and circumferential shifts, and the sections are stacked accordingly, resulting in the final 3D blade shape [28]. For each profile contour definition, a section build-up based on the superposition of a camber-line angle (x) and a thickness distribution T(x) is applied [29]. In the present work, the normalized distributions plotted in Fig. 4 are parameterized by B-splines [20], using 7 CPs for the non-dimensional camber-line angle distribution and 9 CPs for the normalized thickness distribution, but they are kept constant during the MOO. From the 21 available 2D profiles, 7 are selected as design sections and the rest are obtained by spline interpolation. Thus, the 5 design parameters ( r I , r E , I , E andT max ) per design section and 2 stacking parameters for each section's axial and circumferential shift, excluding hub and tip, sum to 45 design variables in total. Using such a CAD-based parameterization, the geometric constraints of the case can a priori be guaranteed by properly bounding the design space. Figure 5 shows the block-structured hexahedral mesh with ∼ 2 million nodes that was generated using an inhouse meshing tool [30], after conducting a mesh refinement study to achieve a mesh-independent computation of the objectives. The maximum y + value of the first series of A steady-state RANS compressible flow solver with the SA turbulence model is used to compute the flow. This is part of the Rolls-Royce in-house CFD suite of codes [32], which have been extensively validated and applied to several industrial applications [33]. The non-linear flow solver uses a vertexcentered finite-volume discretization scheme and the 5-stage Runge-Kutta method [34] (for the pseudo-time-stepping) is accelerated by a block-Jacobi preconditioner and a geometric multigrid technique [35]. Figure 6a illustrates the close-to-surface streamlines after the flow simulation for the BS geometry has converged. It can be seen that corner flow separations occur at both hub and tip which, then, expand along the suction side of the blade. This is common for highly loaded compressor stators and it was also experimentally observed in [36], where this BS geometry was manufactured and measured. These separations contribute to a loss of 0 = 0.0764 , while the flow turning achieved by the BS blade is such that E,0 = 3.74 • . In [36], the presented CFD setup was experimentally verified as an accurate simulation tool for this compressor stator application.
The HYDRA discrete adjoint solver [17] is used to provide the volume sensitivities (i.e. the change in objective function w.r.t. a volume mesh node perturbation). Surface sensitivities (i.e. the change in objective function w.r.t. a surface mesh node perturbation) are computed by applying the inverse operation of a spring-based mesh deformation algorithm to the volume sensitivities [38]. Then, the surface sensitivity map is obtained by projecting each surface node's sensitivity vector onto its corresponding boundary normal. This is visualized for the BS in Fig. 6b where is used as the objective. The sensitivity map shows which areas should be modified to suppress the corner flow separations and, thus, reduce the losses. Here, blue areas should be pulled outwards and red pushed inwards to achieve a decrease in . The black line depicts the zero sensitivity line.
Coupling , where b denotes the blade design parameters. The accuracy of these gradients was verified as in the previous study [37]. The required computation of geometric sensitivities dX S db for each design parameter is done using second-order FD. Although the computational cost of FD for the geometry part scales with the number of design parameters, the total cost for the 45 parameters used in this case is negligible compared to the cost of the primal and adjoint solutions.

Computing the initial Pareto point
The first Pareto point required to initialize the MOO process is computed as described in step 1 of Algorithm 1. Hence, an initial optimal stator ( OPS 0 ) blade is sought, which is obtained by minimizing exclusively w.r.t. the first objective . To perform this SOO for the "Go-to-Pareto" step, an automated workflow was developed within the Rolls-Royce design process for this compressor stator application. This SOO workflow is also employed during the correction steps of the "Move-on-Pareto" process ( Fig. 7), using a different minimization function.
The BFGS optimization method has come up with a leaned (bow-shaped) blade; Fig. 8-(b), which is more loaded at mid-span and less loaded at hub and tip, thus reducing the (still existing) corner flow separations. This is also illustrated in Fig. 9, where the total pressure loss coefficient distribution at an axial cross-section just behind the TE is plotted for the BS (a) and OPS 0 (b) blades. Especially the tip corner separation is almost fully suppressed. Thickness minimization (up to the bound value), LE/TE optimization On the other hand, the exit flow angle deviation from the axial direction rises to * E,0 = 6.7 • , which is mainly driven by its circumferential/whirl component w . The whirl angle distribution at the exit plane of the CFD domain is plotted in Fig. 10 for the BS (a) and OPS 0 (b) blades and is primarily affected by the modifications in E .

Pareto tracing
Once the first Pareto point has been computed, the proposed methodology following Algorithm 1 is employed, while requesting for 10 steps ( n s = 10 ) each of which attempts to reduce E by 0.5 • ( ̂E = −0.5 • ). Figure 7 shows the 11 Pareto optimal points which were obtained in the non-dimensional objective space. These nondominated solutions are in general well distributed along the Pareto front, apart from the final point b * 10 which is further away from the previous point b * 9 due to the "flattening" of the front. An intermediate point could be found (if that were of interest) by decreasing the step size ̂E . The predicted points are close to the optimal solutions and, thus, are omitted from the plot (similar to the target points) for the sake of clarity. The point represented as a red triangle was obtained (at a preliminary stage) by performing an additional "Go-to-Pareto" optimization with [w 1 w 2 ] = [0 1] , and seems to be in line with the shape of the computed Pareto front.
In terms of computational cost, the "Move-on-Pareto" process required in total 159 EFS for the primal and adjoint simulations (including line search steps). This consists of 21 EFS for the first step (without prediction) plus ∼ 15 EFS on average for each subsequent step (using prediction-correction), which is considered to be efficient when compared to the 36 EFS needed for the "Go-to-Pareto" step. A summary of this cost comparison between the components of the method is tabulated in Table 1, where the benefit of incorporating a prediction step is shown.
All the Pareto points in Fig. 7 share minimum values for the design parameters r I , r E and T max of all blade sections excluding hub and tip. The continuous geometry modifications when moving from point to point are attributed to variations in the blade angles I , E (controlling the loading of the blade), as well as to the 3D stacking of the sections (axial and circumferential shift parameters). This is illustrated in Fig. 8, where the geometries corresponding to the Pareto points b * 4 (c), b * 8 (d) and b * 10 (e) are compared to the BS (a) blade. The progressively formed wavy TE is mainly responsible for reducing the exit flow angle to a minimum of * E,10 = 1.98 • , as shown in Fig. 10, where the distribution of the whirl angle component at the exit plane of the CFD domain is plotted for the OPS 4 (c), OPS 8 (d) and OPS 10 (e) blades. The OPS 10 blade achieves an average whirl angle closer to zero.
Finally, the total pressure loss coefficient distribution at an axial cross-section just behind the TE is plotted in Fig. 9 for all the discussed blades. The one corresponding to OPS 10 (e) is slightly worse than the BS one (a) and results in * 10 = 0.077. Performing the presented MOO is clearly beneficial to the design process, since it simplifies the task of the design engineer, who can select one of the computed non-dominated solutions (e.g. the OPS 8 which dominates the BS w.r.t. both objectives: * 8 = 0.0746 , * E,8 = 2.73 • ) according to his/her preference.

Conclusions
In this paper, a new adjoint-based method for the solution of bi-objective optimization problems was developed. The method relies on an iterative prediction-correction approach for efficiently tracing the Pareto front. The main difference from other curve continuation techniques exists in the Fig. 7 Bi-objective optimization of a 3D compressor stator: Pareto front in non-dimensional objective space approximated by 11 welldistributed optimal points obtained using the proposed methodology. The red filled triangle stands for an additional Pareto point, which resulted from a "Go-to-Pareto" optimization with [w 1 w 2 ] = [0 1]  formulation of the KKT conditions, where a distance function is treated as the objective. A BFGS optimization method was implemented to approximate the Hessian information required at each prediction step of the algorithm.
The proposed method was first demonstrated in the lift/ drag optimization of an isolated airfoil and then successfully applied to the bi-objective optimization of a 3D compressor stator. It was found that the BFGS Hessian approximations can further be utilized to accelerate the "Move-on-Pareto" technique. Here, the introduction of the prediction step seems to increase method's efficiency by a factor of ∼ 1.4 (average value for examined cases). As a result, the presented Pareto fronts were obtained in an automatic manner with low computational cost and well-distributed solutions. The extension of the method to multidimensional MOO problems is straightforward, provided that an algorithm to set the targets (for sweeping the Pareto front) is found.