{'fraction_non_alphanumeric': 0.02795316496163683, 'fraction_numerical': 0.00851182864450128, 'mean_word_length': 4.752097460062062, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Design optimization using multiple dominance relations

A challenge in engineering design is to choose suitable objectives and constraints from many quantities of interest, while ensuring an optimization is both meaningful and computationally tractable. We propose an optimization formulation that can take account of more quantities of interest than existing formulations, without reducing the tractability of the problem. This formulation searches for designs that are optimal with respect to a binary relation within the set of designs that are optimal with respect to another binary relation. We then propose a method of finding such designs in a single optimization by defining an overall ranking function to use in optimizers, reducing the cost required to solve this formulation. In a design under uncertainty problem, our method obtains the most robust design that is not stochastically dominated faster than a multiobjective optimization. In a car suspension design problem, our method obtains superior designs according to a k‐optimality condition than previously suggested multiobjective approaches to this problem. In an airfoil design problem, our method obtains designs closer to the true lift/drag Pareto front using the same computational budget as a multiobjective optimization.


INTRODUCTION
If an engineering design problem can be expressed as an optimization problem, then optimization algorithms or heuristics can be exploited to efficiently find desirable designs. 1,2 An optimization problem is to find designs that are optimal with respect to particular binary relation * . Common optimization formulations are single-objective minimization, where this binary relation is induced by using less than or equal to in order to compare scalar objectives, and multiobjective optimization, where this binary relation is induced by using Pareto dominance to compare vectors of objectives (and the performance of optimal designs is known as the Pareto front). In these common optimization formulations, a quantity of interest must either be assigned as an objective or a constraint if it is to be included in the problem. Often a quantity does not clearly belong to either of these two categories: it may not be feasible to use additional objectives in a problem that is already computationally expensive to solve, but an appropriate threshold by which to constrain a quantity may not be known. A designer hoping to use optimization for their design problem is faced with a dilemma.
For example, consider transonic airfoil design, in which key design drivers are increasing lift and reducing drag, but other important quantities of interest are trailing edge separation and pitching moment. If all four quantities were used as objectives in a multiobjective optimization, trade-offs would exist between all objectives and the optimizer would be trying to explore a large area of design space. Under a strict computational budget, the optimizer will not be able to explore these trade-offs sufficiently to obtain satisfactory designs. An alternative is to constrain the amount of separation and the pitching moment. The appropriate threshold to use for these constraints may be difficult to determine a priori, since they are physical quantities whose value affects other aspects of the overall system. Too strict a constraint can mean there are no feasible designs, and too loose a constraint is likely to be exploited by the optimizer to give airfoils that are unrealistic for the overall system.
We propose that in such a scenario, it would be useful to search for designs on the Pareto front of separation and moment within designs on the Pareto front of lift and drag. This allows the separation and moment to be taken into account without increasing the number of objectives or constraints in the problem, which would reduce computational tractability. However, with current optimization methods, obtaining such designs would require a multiobjective optimization of lift and drag to be performed first, and then moment and separation to be taken into account in a second step, potentially leading to a bi-level optimization problem.
Consider also the problem of designing a car suspension system to maximize ride comfort and road holding. This is a problem that is known to have many conflicting objectives. Various different multiobjective (and many-objective) optimization formulations have been applied to the design of suspension systems, 4,5 but the most appropriate objectives to use is not well established. 6 It would be desirable to consider all the common objectives without having to resort to a large optimization problem, or choose arbitrary values to constrain them by.
A related type of design problem is one in which we wish to consider vector-valued quantities instead of scalar quantities of interest. For example, when the performance of a system is uncertain, we may want to consider the cumulative distribution function (CDF) of the quantity of interest. A design should be both robust (often quantified by low variance of the quantity of interest) but also give a high probability of good performance. In many cases, minimizing variance can result in designs that are robust, but whose CDFs lie entirely at worse values than other designs; such designs are said to be stochastically dominated. 7,8 Therefore, we would be interested in obtaining the most robust design that is not stochastically dominated, which is equivalent to finding the design with the lowest variance within the set of stochastically nondominated designs.
These examples (and many other design problems) can be formulated as finding designs that are optimal with respect to a binary relation within the set of designs that are optimal with respect to another binary relation. In this article, we develop a method of finding such designs in a single optimization, by proposing an optimization formulation of which they are the solution. As a result, such designs can be obtained more efficiently than existing methods.
In Section 2, we mathematically define the type of design problem we are considering, propose a binary relation whose optimal designs are solutions to this problem, then discuss how this can be used within optimizers to efficiently solve this problem. In Section 3, we apply this formulation to a design under uncertainty problem, to illustrate its relevance to vector-valued performance measures. Then in Sections 4 and 5, we apply this formulation to the two motivating examples discussed in the introduction to illustrate the benefit of being able to solve this type of engineering design problem efficiently. Finally in Section 6, we conclude the article and discuss avenues for future work.

DESIGN OPTIMIZATION USING MULTIPLE DOMINANCE RELATIONS
This section defines the problem formulation, proposes a binary relation for solving the problem in a single optimization, develops a method of using this relation in several optimization algorithms, and discusses how this accelerates the solution of the design problems discussed in the introduction.

Problem definition
We denote a single design as x, which is an element of the set of all admissible designs X (taking into account any bounds or constraints); we refer to X as the design space. We denote the performance of a design q ∈ Q, where Q is the set of values q can take and is referred to as performance space; q is given by a function f (x) ∶ X → Q. We define the performance q such that smaller is better, and so q is to be minimized. The performance, q, can represent any or all information from F I G U R E 1 Illustration of a binary relation ⪯ giving a preorder over design space a simulation analyzing design x and so could contain scalar objectives, one-dimensional functions (eg, pressure as a function of distance over an aerofoil), higher dimensional functions (the pressure as a function of 3D spatial coordinates), and so on. We define an optimization problem (similarly to References 9,10) as the problem of finding the set of optimal designs in X with respect to a binary relation ⪯. A binary relation is a comparison of two designs that determines whether one dominates the other and is usually induced by a comparison of the performance of the designs. An optimal design in X with respect to a relation, denoted x * , is one for which no other design dominates it (and so is also referred to as a nondominated design). Note that for optimal designs with respect to a binary relation to exist, the binary relation must be a preorder † over X. Definition 1. The set of optimal designs in X with respect to ⪯, denoted min(X, ⪯), is given by: (1) Figure 1 illustrates a binary relation ⪯ defining a preorder over the design space given by X = {a, b, c, d, e}. A design from which an arrow originates is dominated by the design at the end of that arrow according to ⪯. In this case, a ⪯ b, d ⪯ a, d ⪯ c, c ⪯ b, and e ⪯ b (also note that because ⪯ is a preorder and so is transitive, d ⪯ b). According to Definition 1, an optimal design with respect to ⪯ is one from which no arrows originates, so in this case X optimal = {d, e}. Note that two designs between which a path along arrows does not exist are incomparable ‡ , so d||e, a||c, a||e, and c||e.
To illustrate further, consider two standard optimization formulations: single-objective and multiobjective optimization. Example 1. Single objective-optimization considers a single scalar value to minimize, that is, q = f (x) ∶ X → R, and compares the performance of designs using the less than or equal to binary relation. This induces the binary relation ⪯ so over

Example 2.
Multiobjective optimization considers m multiple values to minimize, that is, q ∈ R m . The performance q is now a vector of quantities of interest given by m functions, such that Multiobjective optimization searches for optimal designs in X with respect to the (weak) Pareto dominance relation, denoted ⪯ Pareto , which is defined in Definition 2.

Definition 2.
For two designs x and y, x weakly Pareto dominates y, denoted x ⪯ Pareto y, if it is at least as good in all objectives, that is, The performance vectors of the optimal designs with respect to ⪯ Pareto are known as the Pareto front, denoted PF(X):  Figure 2 illustrates how the Pareto dominance binary relation induces a preorder over design space (the same preorder illustrated in Figure 1). Comparing the performance of each design (given by the values of two objectives on the right) using Pareto dominance induces the arrows between designs on the left.
In this article, we are interested in finding the optimal set of designs with respect to a binary relation within the optimal set of designs according to another binary relation. For this two-relation case, this set of designs can be written as: Example 3. For example, consider the following algebraic function, for which q = [f 1 , f 2 , f 3 ], over the 2D design space X = [−5, 5] 2 : where ||.|| indicates the L 2 norm. Let ⪯ 1 indicate Pareto dominance for f 1 and f 2 , and let ⪯ 2 indicate less than or equal to for f 3 . We are looking for the following set: min(min(X, ⪯ 1 ), ⪯ 2 ). This set consists only of the design with the smallest f 3 within the set of designs that are Pareto optimal for f 1 and f 2 , and this design can be found analytically to be x = [3,0]. Some designs (and their corresponding performance) that are optimal with respect to ⪯ 1 are plotted in gray in Figure 3, and the design among them with the lowest f 3 -the single design in min(min(X, ⪯ 1 ), ⪯ 2 )-is highlighted in red; contours of f 3 are given in blue.
In general, we consider that we could nest the formulation of Equation (4) further using n binary relations, denoted ⪯ i , i = 1, … , n, and look for the set given by Definition 3.

Definition 3.
The set of optimal designs in X according to n binary relations denoted ⪯ i , i = 1, … , n, is given by: Definition 3 is a problem formulation in that it specifies we are looking for X optimal , but says nothing about how we might obtain designs in this set. In the following sections, we develop an efficient method of obtaining designs given by Equation (5) in a single optimization.

An overall binary relation using a hierarchy of ranking functions
To use an optimization to find X optimal from Equation (5), we need an overall binary relation, denoted ⪯ (1,2,… ,n) , such that min ( X, ⪯ (1,2,…,n) ) = min ((… min (min (X, ⪯ 1 ) , ⪯ 2 ) …) , ⪯ n ). We propose such a relation by using a ranking function for each of ⪯ k , k = 1, … , n as a measure of how far a design is from true optimal designs, and use lexicographic ordering of these ranking functions: Definition 4. The binary relation ⪯ (1,… ,n) , referred to as the combined binary relation of ⪯ k , k = 1, … , n, is given by: where each R k , k = 1, … , n is a function that satisfies the following two conditions (and is referred to as a ranking function): 1. Comparison of R k gives a necessary condition for ⪯ k , that is, This definition can be written more concisely using recursion of the index k in the notation ⪯ (k,… ,n) : where x ⪯ ∅ y is true ∀x, y ∈ X.
In the Appendix, we give proofs that the combined binary relation in Equation (6) is a preorder, and that min(X, ⪯ (1,…,n) ) is the set X optimal defined in Equation (5)-showing that it achieves our goal for this article.

Ranking of continuous problems and discrete approximations
In order to use the combined binary relation ⪯ (1,2,… ,n) , we need a ranking function for each ⪯ k that satisfies the conditions specified by Definition 4. When ⪯ k is the Pareto dominance binary relation, one such ranking function is the Euclidian distance in Q of the performance of a design from that of its nearest optimal design; another would be the volume enclosed by its performance, a hyperplane perpendicular to each objective, and the Pareto front. These are referred to as distance ranking and volume ranking, respectively, and are illustrated in Figure 4 for a single ⪯ k consisting of Pareto dominance of two objectives in a continuous design space. The issue preventing these ranking functions from being used directly in ⪯ (1,… ,n) is that they cannot be evaluated without knowing the optimal designs with respect to ⪯ k , which are not known a priori (otherwise we have already solved the problem). Note that while scalarization methods 11 -which use a combination of objectives to define a scalar function whose minimizer is an optimal design of Pareto dominance, for example, weighted sum-satisfy condition 1, they do not satisfy condition 2, and so are not suitable to use as ranking functions for our goal.
The purpose of defining ⪯ (1,… ,n) is so that optimizers can be used to efficiently find designs in X optimal defined by the problem formulation in Equation (5). Therefore, for practical purposes we care not about the actual value of a ranking function for a given design, but rather how the ranking function is used to compare designs visited during an optimization to drive the optimizer toward designs in X optimal .
We can evaluate approximations to these rankings by only considering a finite set of designs, denotedX, and treating this as a finite representation of the (potentially infinite) design space X. Therefore, we use the following binary relation within optimizers: whereR k is an approximation to R k obtained by only consideringX. As an optimizer explores design space, these approximations drive the optimizer to evaluate designs closer to designs in X optimal , which in turn will make these ranking approximations more accurate, achieving our goal of finding designs in X optimal efficiently. An approximation to the distance ranking is the nondominated sorting method of Reference 12, where the rank of x is the number of times the current nondominated set needs to be removed from the total set before x would be in the nondominated set. Algorithm 1. DeterminingR k (x) for relation ⪯ k using a distance approximation, given a set of pointsX An approximation to the volume ranking can be obtained by counting the number of other designs which dominate any given design. 13 Algorithm 2. DeterminingR k (x) for relation ⪯ k using a volume approximation, given a set of pointsX These two approximations are illustrated in Figure 5, where each point is labeled by the rank it would be assigned using each of these approximations.
Using a ranking function over a finite approximation to design space is commonplace in population-based evolutionary optimizers (and is why we use the terminology ranking functions), 14 where the current population is used as the finite approximation of design space over which a ranking function is evaluated. For example, assigning fitness in multiobjective optimization by counting the number of dominating solutions in a population is known as Pareto ranking, 13 and the popular NSGA-II algorithm from Reference 15 uses nondominated sorting to assign fitness. A ranking function in these contexts is used as an operator intrinsic to the workings of a particular algorithm. By contrast, here we use a hierarchy of an arbitrary number of ranking functions in order to solve our proposed optimization formulation of Definition 3. Furthermore, we are considering a general context that is independent of any particular algorithm: we first derive the conditions on ranking functions required for continuous problems, and then use approximations on a finite representation of design space. A population within an evolutionary algorithm is just one such representation. Our approach is independent of the particular optimizer used to solve the problem. Indeed, in the entropy-based Bayesian optimization methods of References 16 and 17, a finite number of representer points are used as an approximation to the overall design space, so our approach is well suited to such optimizers.
Using an approximation to a ranking function, we evaluate our proposed binary relation as follows: Then we can obtain the optimal designs within a given finite set of pointsX as follows:

Implementation of the combined binary relation within optimizers
Here we adapt two optimization heuristics in particular to solve our proposed optimization formulation: a multiobjective genetic algorithm (NSGA-II) 15 and a multiobjective particle swarm (MOPS). 18 When using evolutionary algorithms, we could use only the current population in optimizers asX to get approximations to these ranking functions, but this can lead to deterioration-designs visited later in an optimization being dominated by previously visited designs. 19 Therefore in this article, we use the set of every design visited during the optimization asX.
In general the set X optimal from Equation (5) will consist of multiple designs, and so population based optimizers that search for a set of designs in a single run are well suited to solving the optimization formulation considered in this article. Furthermore, engineering design considers a wide range of applications whose design spaces often lack exploitable structure (linearity, convexity, etc), 1 and such optimization heuristics are able to effectively explore such design spaces.
The two chosen optimizers (NSGA-II and MOPS) can be adapted to use the combined binary relation of Equation (8) via the following two modifications: • After the system model has been evaluated for each new candidate population, we evaluate the ranking function of every design visited so far during the optimization for each ⪯ k , k = 1, … , n, using Algorithm 1 or 2. • Whenever a design is normally compared using Pareto dominance, we instead use the combined binary relation of Equation (8).
• In NSGA-II, Algorithm 4 is used when selecting which designs to breed (the selection step) and which designs to keep in subsequent populations (the replacement step). • In MOPS, Algorithm 4 is used when finding the local and global best particles.
The computational cost of ranking each design inX using each of these approximations is O(nMN 2 ), where N is the number of designs inX, n is the number of individual binary relations used, and M is the number of values that define each relation. In many engineering design problems, the computational cost of evaluating f(x) even once will outweigh the cost of doing this ranking, since we can usually only afford N ≃ O(100) to O(1000) evaluations for an optimization.
A limitation of this implementation is therefore memory usage. It relies on building up an approximation to design space (X) over the course of an optimization, so requires all previously visited designs and their performance to be recorded. For practical engineering design problems, the limitation of the number of design points that can be explored is usually dominated heavily by the time taken to run the simulation, and so this should not be a concern. However, if an optimization problem is expected to be run for a very high number of iterations, then storing every previously visited design point may not be practical, and so alternative approaches may be preferable. This is an interesting avenue for future work: which design points should be archived, and can we achieve convergence with a finite archive size?

Illustrative example
We could find designs in X optimal from Equation (5) by doing a standard multiobjective optimization of ⪯ 1 , then choosing designs within this set according to the other binary relations, but doing so can be prohibitively costly when evaluating f(x) is computationally expensive. The purpose of defining ⪯ (1,… ,n) and its approximation in Equation (8) is to accelerate the convergence of optimizers toward designs in X optimal , allowing such designs to be found in feasible computational times in challenging engineering design problems.
To demonstrate the acceleration this provides to optimizers, consider the following 4D version of the test problem used in Example 3, over a design space X = [−5, 5] 4 : where ||.|| indicates the L 2 norm. Similarly to Example 3, the best design (the design in the singleton set X optimal ) is the design in the Pareto front of f 1 and f 2 with the lowest value of f 3 , which can be found analytically to be x = [3, 0, 0, 1]. We consider using the following three binary relations in an optimization to find this design: Pareto dominance of f 1 and f 2 ; Pareto dominance of f 1 , f 2 , and f 3 ; the combined binary relation ⪯ (1,2) with ⪯ 1 indicating Pareto dominance of f 1 and f 2 , and ⪯ 2 indicating less than or equal to using f 3 .
We use a population size of 50 with a computational budget of 1000 designs evaluated, and run each optimizer 20 times using the same set of randomly generated initial design points for each of the three cases. Figure 6 shows convergence in terms of average Euclidian distance from x = [3, 0, 0, 1]. To evaluate the combined binary relation, we use both the ranking function that counts the number of times the current nondominated set would be removed before a design is optimal as an approximation to distance ranking (labeled sorting), as well as the ranking function that counts the number of designs that dominate a design as an approximation to volume ranking (labeled number).
The acceleration in convergence is brought about because fewer designs are incomparable according to the combined binary relation ⪯ (1,… ,n) than just ⪯ 1 , increasing the selection pressure exerted on an optimization algorithm to explore in a particular direction of design space. 20 This is because, given a collection of candidate designs, fewer of them are optimal, so an optimizer will consider fewer search directions in which to search for better designs.
The benefit of the proposed formulation is similar for both ranking functions, therefore for the remainder of this article we just use the number approximation, as it is simpler to implement.
Note that the effect of this increase in selection pressure by reducing the number of incomparable designs is well recognized in the field of many-objective optimization, which has proposed various alternatives to Pareto dominance that extend better to more objectives. [21][22][23][24][25] Often this is achieved by trying to obtain information from a designer to quantify their underlying preferences using techniques from the field of multicriteria decision making, 26-28 such as goals, 29 reference points, 30 trade-offs, 31 and so on. Physical programming 32,33 similarly elucidates information from a designer to convert a multiobjective problem into a single-objective one.
By contrast, this article is not trying to quantify designers' preferences (which can be difficult to do, as exemplified by the vast range of techniques available to try to achieve this); instead we propose that some design problems (such as the motivating examples given in the introduction) are usefully expressed in such a way that cannot be solved efficiently using existing optimization formulations, so we define the combined binary relation ⪯ (1,2,… ,n) to efficiently solve such design problems using optimizers.
Still, a limitation is that the designer still has to decide on which binary relations to use within this formulation, and a hierarchy to place them in. The approach offers more power through more flexibility, but it may be a challenge to set up appropriate optimization problems.

Constraint handling
In the problem definition, we defined X to implicitly include all constraints on the design, but it is straightforward to include constraints in the formulation by setting the first binary relation ⪯ 1 to be any preorder whose optimal designs are the feasible set of designs. This includes many of the constraint handling approaches that have been proposed for evolutionary algorithms. 34 Denoting such a binary relation as ⪯ constraint , and the set of feasible designs as X feasible , our optimal set of designs defined in Equation (5) becomes: Note that ranking methods have previously been proposed to handle constraints in evolutionary algorithms. For example constraints in Reference 35 are handled by ranking populations according to how many other designs are more feasible in more constraints than each design in the population; the constraint handling method in Reference 36 ranks populations according to constraint violation and according to Pareto dominance, and sums these ranks to assign fitness; and the method in Reference 37 ranks populations according to each constraint function separately, then sums these ranks to determine the best infeasible solutions to include in the next population. These methods are specific to constraints with fixed thresholds and particular evolutionary optimizers, whereas in this article, we consider the more general problem of finding the set X optimal defined in Equation (5).
Having formulated an optimization problem to find these designs, and discussed how optimizers can solve this formulation using finite approximations to design space, in Section 3 we apply this formulation to an algebraic design under uncertainty problem to illustrate its relevance to vector-valued performance measures. Then in Sections 4 and 5, we apply this formulation to the two motivating examples discussed in the introduction to illustrate the benefit of being able to solve this type of engineering design problem efficiently.

DESIGN UNDER UNCERTAINTY
In design under uncertainty problems, using a few scalar objectives such as statistical moments can miss out important behavior of the full probability distribution of the quantity of interest. This has motivated the development of several optimization formulations that consider the probability distribution in its entirety. [38][39][40][41] Similarly in Reference 8, where we used multiple dominance criteria whose optimal sets were guaranteed to have a nonempty intersection, giving a trivial overall preorder over design space to use in optimization. For example, if the CDF of a design lies entirely to the left of the CDF of another design, then the first design has a higher probability than the second design of achieving a lower value than any given value of the quantity of interest (the CDF of an uncertain quantity of interest X, evaluated at x, gives the probability of a X taking any value less than x). This is known as stochastic dominance and is not taken into account if only statistical moments are considered, as illustrated in Figure 7, where moments and the CDF are shown for four hypothetical designs. Designs A, B, and D are all on the Pareto front of moments, but design D is guaranteed to give a worse quantity of interest than design A or B. Furthermore, the CDF of design B does not cross that of design A, so is stochastically dominated.
In many design cases, we would not be interested in stochastically dominated designs, but robustness (indicated by low values of the standard deviation of the quantity of interest) is still a desirable quality, so we may want to find the most robust design that is not stochastically dominated. We cannot formulate this problem using the approaches of Reference 8, but can using the approach developed in this article. To demonstrate, we consider the following test problem: where x ∈ [−0.5, 0.5] 2 is the vector of design variables that can be controlled, and u is the realization of an uncontrollable random input, denoted U, which is uniformly distributed over [−1, 1] 2 . We estimate the behavior of g under uncertainty using M = 1000 independent and identically distributed (i.i.d) sampled values from the distribution of U, denoted u i , i = 1, … , M. The function g(x, u) is evaluated for a given x at each sampled value of U to give sampled values of the output, denoted g i and given by From these sampled values we use the following objectives: where g (j) i denotes the jth smallest value out of the sampled values g i , i = 1, … , M. We can then search for optimal designs of ⪯ (1,2) with ⪯ 1 indicating Pareto dominance of f j cdf , j = 1, … , M, and ⪯ 2 indicating lower (or equal) f std . If we were to consider this as a many objective problem, we are dealing with 1002 separate objectives, and so trying to quantify a designer's preferences in order to weight these objectives would be infeasible.
We illustrate the benefit of our approach by running an optimization (particle swarm) with a strict computational budget (150 design points considered) finding optimized designs of three binary relations: Pareto dominance of moments (f mean and f std as objectives), stochastic dominance (f j cdf , j = 1, … , M as objectives), and the proposed combined binary relation ⪯ (1,2) . These represent three methods of finding the design specified in the problem formulation: that with the lowest f std that is stochastically nondominated. Figure 8 gives the optimized designs when just using moments, Figure 9 gives the optimized designs when using stochastic dominance, and Figure 10 gives the optimized designs when using the combined binary relation. In each case we highlight in green the nondominated designs with respect to stochastic dominance (labeled SD Front), the Pareto nondominated designs in red (labeled Pareto Front), and the design on the stochastically nondominated front with the lowest f std out of those visited during the optimization in blue (labeled Best Design).
Using this computational budget, the stochastic dominance-based optimization focused on an area of design space giving moments in the upper left quadrant of Figure 9A, and so has found many designs with poor mean and variance but with a small chance of good performance (with long left-tailed CDFs). This has caused it to be unable to obtain designs close to the true best design under this computational budget. Similarly, the moment-based optimization focused  Figure 8A, where designs with low variance are located (but which we can see from the CDFs are stochastically dominated), also causing it to miss the designs specified in the problem formulation. By contrast, by using the combined binary relation, the optimization was able to focus in around designs of most interest using the given computational budget.
To quantify this behavior, we run 20 of these optimizations, and plot the average distance of the best point found so far by the optimizer to the true best design vs number of function evaluations in Figure 11.
For both optimizers, using the combined binary relation accelerated the convergence of the optimizers toward the best design specified in the problem formulation. Considering the average distance for the moment-based formulation after the budget has been exhausted, using the combined binary relation achieves this average distance in under half the number of function evaluations, and in fact with the particle swarm optimizer we see computational savings of ≃70%.

HALF CAR MODEL
To investigate the usefulness of this formulation for engineering design, we first consider the problem of designing a car suspension system to maximize ride comfort and road holding. This is a problem that is well known to have many conflicting objectives, and which are the most appropriate objectives to use in optimizations is not well established. 6 Therefore, multiobjective optimization approaches have been proposed for the design of suspension systems. 4 The difficulty in choosing objectives effectively has also led to techniques from many-objective optimization being proposed. 5 Here we consider the half-car suspension system design problem from Reference 5, and compare our approach to this design problem with the three-objective optimization approach proposed in that article. This design problem considers two vibrating systems: a half-car suspension system modeled by mass-spring-damper subsystems for the road, car body and tires, and a person located above the center of gravity of the car modeled by a collection of mass-spring-damper subsystems. These two systems are treated as uncoupled, such that the response of the car to the road is used as a separate input to drive the response of the person.
The equations of motion for these two systems can be derived from the diagram given in Figure 12, as detailed in Reference 5. We use the same values for the stiffness, mass, and damping parameters as in Reference 5. We examine the response of the car to driving over a bump in the road represented by a half sinusoid of height 0.015 m, as illustrated in Figure 13.
There are six common objectives to use when designing suspension systems of this type, [4][5][6] given below: where var(.) in this notation refers to variance, z S is the displacement of the car body, ST F is the suspension travel for the front wheel, TD F is the tire deflection for the front wheel, VDV head is the vibration does factor at the head, equal to (∫z 4 1 dt) (1∕4) (where z 1 is the displacement of the head in the body model), and CF head is the crest factor at the head, equal to max(z 1 )∕rms(z 1 ) (where rms denotes the root mean squared). Objectives f 1 , f 4 , and f 5 are measures of the ride comfort, and objectives f 2 , f 3 , and f 6 are measures of the road holding of the vehicle.
In Reference 5, the authors propose only optimizing objectives f 1 , f 2 , and f 3 to keep the problem of a reasonable size, and then using the k-optimality approach of Reference 42 over all six objectives as a postoptimization analysis to determine the most desirable designs to present as the results of the optimization to a designer.
As an alternative, we define three binary relations: Pareto dominance using objectives f 1 and f 2 (denoted ⪯ 1 ), Pareto dominance using objectives f 3 and f 4 (denoted ⪯ 2 ), and Pareto dominance using objectives f 5 and f 6 (denoted ⪯ 3 ). We use optimization to search for designs in the following set: min(X, ⪯ (1,2,3) ).
These binary relations are chosen because the main trade-off for this problem is between ride comfort and road holding, and so the first binary relation uses the two primary objectives that measure each of these two features, then the second and third binary relations are chosen as one more each from the objectives that measure these features.
The design variables for this optimization are given in Table 1, along with their bounds and initial values. We run the NSGA-II algorithm with a population size of 50 for 30 generations, for three optimization formulations, and plot the objectives of the optimized designs in Figure 14. The first case uses Pareto dominance of the three objectives f 1 , f 2 , and f 3 , the second case uses Pareto dominance of the four objectives f 1 , f 2 , f 3 , and f 4 , and the third cases uses our proposed approach as described above (labeled "Search"). The same initial population is used for all three cases.  Figure 14, we can see that the three objective optimization and our approach obtained designs on two separate sections on this Pareto front. However, when looking at the objectives f 3 and f 4 (center figure), it is more evident that our approach was focusing on a smaller area of design space, and has consequently found designs that Pareto dominate those found by the other two cases under objectives f 3 and f 4 . The objectives f 5 and f 6 were not considered by the two standard multiobjective optimizations and were considered at the bottom of the hierarchy in our approach, but the right figure still illustrates how the different cases ended up with designs in different regions of design space.
To investigate the time and frequency response of designs resulting from these optimizations, one design from the middle of the Pareto front for each of the three optimizations are selected. The time series response and frequency spectrum are given forz S in Figure 15 and for z F in Figure 16.
Furthermore, we evaluate the k-optimality of the optimized solutions (an approach proposed in Reference 42, and used in Reference 5 for this half-car problem): a design is k-optimal if it is not Pareto dominated by any other solution for all possible k-element subsets of objectives. The maximum k-optimality value of any solution in blue is 3, whereas the maximum of any solution of the other two fronts is 2. This suggests that a superior design has been obtained using our combined binary relation than the other two optimization formulations.

TRANSONIC AIRFOIL DESIGN
In transonic airfoil design, key design drivers are increasing lift and reducing drag, but other important considerations are trailing edge separation and pitching moment. Low trailing edge separation ensures that the aerofoil is not close to buffet (a dangerous aeroelastic instability in the flow). Low pitching moment ensures that drag induced by the need to trim the aircraft (counteracting this pitching moment using the tail) is also low. As discussed in the Introduction, the appropriate threshold to use if separation and pitching moment were constrained is difficult to know a priori, and using all four objectives in a 4D multiobjective optimization may be too computationally expensive. Therefore, it would be useful to be able to take account of separation and pitching moment in an optimization in such a way that does not require them to be constrained, and that can obtain satisfactory designs under a strict computational budget.
We propose this can be achieved by searching for designs on the Pareto front of separation and moment within the designs on the lift/drag Pareto front. The formulation developed in this article can do so efficiently by using optimizers to search for optimal designs with respect to ⪯ (1,2) , with ⪯ 1 indicating Pareto dominance of lift and drag and ⪯ 2 indicating Pareto dominance of trailing edge separation and pitching moment. In this section, we demonstrate how this formulation can obtain airfoils with better characteristics than a standard multiobjective optimization formulation using the same computational budget.

Setup
We consider designing an airfoil at a fixed angle of attack of 3 • , a free stream Mach number of M = 0.8, and a Reynolds number of 6.5 × 10 6 . We use the open source CFD solver SU2 (along with an Spalart-Allmaras turbulence model) to evaluate the performance of a given airfoil, 43 which has been validated on various aerodynamics problems, making it suitable for the purposes of this study. SU2 directly outputs lift, drag, and moment coefficients, and as a measure of separation, we evaluate the total area under the graph of distance along surface vs surface skin friction whenever skin friction is negative. We parameterize a design space by perturbing a baseline mesh (that is provided with SU2) using a series of Hicks-Henne bump functions on the upper and lower surface, using SU2's built in mesh deformation capability. The mesh for the baseline geometry is shown in Figure 17.  We also constrain the shape of the airfoil such that a box of dimensions 0.5 × 0.07 chord lengths (which just fits inside the baseline geometry) must fit inside any new airfoil, to represent structural requirements.
We use 10 design variables that represent the magnitude of the Hicks-Henne bump functions. The locations on the airfoil and bounds of the 10 design variables are given in Table 2.
We use the particle swarm optimizer with a computational budget of 600 CFD evaluations, and run this optimizer to solve three formulations (using the same initial population in each case): multiobjective optimization of just lift and drag (labeled 2D Pareto), multiobjective optimization of lift, drag, moment, and amount of separation (labeled 4D Pareto), and our combined binary relation using Pareto dominance of lift and drag as ⪯ 1 , and Pareto dominance of moment and separation as ⪯ 2 (labeled combined).

Results
In Figure 18, we plot the designs visited during these three optimization cases as small dots, and we plot the resulting Pareto front of lift and drag as large circles (except in the optimization using of all four objectives, we plot as large circles the Pareto front of lift, drag, separation, and moment). We also plot the minimum drag design from the last few populations of each optimization as a square, and the maximum lift design as a diamond.
Using four objectives in Pareto dominance means that under the given computational budget, the optimizer was not able to converge significantly toward designs on the actual Pareto front of lift and drag, causing it to end up with designs that are dominated in all four objectives.
The main comparison is between the 2D Pareto front of lift and drag and that of the combined binary relation, and there are two principal observations to be made. The combined binary relation focused the search in lift/drag space to only a particular region of the trade-off (lower drag and lower lift), enabling it to converge to better designs in this part of design space, meaning that nearly every design in the Pareto front of lift/drag from this optimization dominates a design on the Pareto front when using just lift and drag in Figure 18A. Second, because we have taken these objectives into account in the combined binary relation, in Figure 18B the designs resulting from this optimization have a better moment and separation than those from the 2D Pareto optimization, also dominating many of those designs (in red) under these two objectives.
This clearly illustrates the advantage of using such a combined binary relation when we have a strict computational budget: we focus the search toward the regions of most interest by taking account of more of the information available to a designer, allowing us to reach the designs of most interest more quickly.
To further investigate the airfoils resulting from these optimizations, we consider in more detail the minimum drag and maximum lift designs from the standard multiobjective optimizations using Pareto dominance, and those from using the combined binary relation. In Figure 19 we plot the airfoil shapes, in Figure 20A we plot the (negative) upper and lower surface pressure distribution, and in Figure 20B we plot the upper surface friction distribution. Figure 19 illustrates that the four selected airfoils have noticeably different shapes, indicating that the optimizers were exploring different regions of design space in each case; Figure 20B illustrates the lower trailing edge separation achieved by the designs found using the combined binary relation compared with those found using Pareto dominance. Furthermore, Figure 20A shows that the minimum drag design found using the combined binary relation has a smaller shock than the other three, which along with the lower skin friction toward the leading edge seen in Figure 20B, explains the improvement in drag and trailing edge separation behavior compared with the other designs. The more desirable pressure and friction distributions exhibited by the designs found using the combined binary relation are a result of the optimization being able to more thoroughly explore this region of design space compared with the Pareto dominance optimization. This in turn is a result of the combined binary relation having fewer incomparable designs than Pareto dominance, driving the optimizer toward this region of design space faster, illustrating the advantage of being able to express design problems in the form of Equation (5) and the utility of the formulation developed in this article for solving such problems efficiently.

CONCLUSION
We have proposed an optimization formulation that searches for optimal designs with respect to a binary relation within the set of optimal designs with respect to another binary relation. In engineering design problems with many quantities of interest, this formulation allows quantities that would otherwise be discarded to be included in an optimization, without reducing the tractability of the problem or constraining the quantities at arbitrary values.
We developed a method of solving this formulation in a single optimization, reducing the computational cost required to find designs that are solutions to this formulation. The method defines a combined binary relation using a hierarchy of ranking functions; the optimal designs with respect to this combined binary relation are solutions to the proposed formulation. We then illustrate how this combined binary relation can be usefully approximated for use within optimizers. This opens up more powerful ways of framing design problems that previously could not be solved efficiently using optimization.
We applied the developed method to three such design problems. First, we applied it to an algebraic design under uncertainty problem: finding the most robust design that is not stochastically dominated. It achieved computational savings of over 70% compared with a standard multiobjective optimization. Second, we applied it to a car suspension design problem, in which it obtained superior designs according to a k-optimality condition than previously suggested multi-objective approaches.
Finally, we applied it to a transonic airfoil design problem in which it is desirable to account for the amount of trailing edge separation and pitching moment but we do not know appropriate values by which to constrain them. Here, given the same computational budget as a standard multiobjective optimizer, the new formulation obtained designs that have both a better lift and drag as well as better moment and separation than the nearest designs found using the multiobjective optimization. This resulted from optimizers using the combined binary relation being able to focus the search onto more desirable regions of the lift/drag Pareto front in order to make better use of the computational budget.
A primary avenue for future work is to investigate the effectiveness of other types of optimizer at solving our proposed formulation, particularly those that already make use of a finite representation of design space (such as entropy search 16,17 ). Furthermore, future work will investigate more efficient single-objective optimizers in the case where we are only searching for a single design (like in the design under uncertainty application). In this case, the combined binary relation gives a total order over design space, and so we should be able to make use of efficient single-objective optimization algorithms to find this design more quickly than population-based heuristics.
In addition, it would be useful to investigate more advanced ways of selecting which designs to include in our approximation of design space when evaluating the approximate ranking functions, since in problems where we can afford more system model evaluations, storing every visited design to use as an approximation of design space might become prohibitively expensive.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.

APPENDIX
Here we show that (a) the combined binary relation (min(X, ⪯ (1,…,n) )) from Equation (6) is a preorder, such that it does not violate transitivity and so min(X, ⪯ (1,…,n) ) exists for any given design space X and (b) min(X, ⪯ (1,…,n) ) is the set X optimal defined in Equation (5). Note that comparison of R k is not a sufficient condition for ⪯ k , that is, R k (x) ≤ R k (y) does not in general imply x ⪯ k y. Condition 2 states within any subset of design space, the ranking of any optimal designs with respect to a given relation must be equal.
The combined binary relation is based on lexicographic ordering, so is a preorder, but we also give a proof here: Lemma 1. The combined binary relation ⪯ (1,… ,n) is a preorder over X, meaning it is reflexive § and transitive ¶ .
of binary relations as follows: x ⪯ (1,…,n) y ⇔ x ⪯ 1 y ∪ ( x || y ∩ x ⪯ (2,…,n) y ) , where x||y indicates x and y are incomparable. However, this binary relation violates transitivity and so it is not a preorder. Therefore, optimal elements with respect to it may not exist and it can give rise to cyclic behavior when used by an optimizer, 9 so is not well suited for use in optimization.