Dominance Move calculation using a MIP approach for comparison of multi and many-objective optimization solution sets

2020·Arxiv

Abstract

Abstract

Dominance move (DoM) is a binary quality indicator that can be used in multiobjective optimization. It can compare solution sets while representing some important features such as convergence, spread, uniformity, and cardinality. DoM has an intuitive concept and considers the minimum move of one set needed to weakly Pareto dominate the other set. Despite the aforementioned properties, DoM is hard to calculate. The original formulation presents an efficient and exact method to calculate it in a biobjective case only. This work presents a new approach to calculate and extend DoM to deal with three or more objectives. The idea is to use a mixed integer programming (MIP) approach to calculate DoM. Some initial experiments, in the biobjective space, were done to verify the model correctness. Furthermore, other experiments, using three, five, and ten objective functions were done to show how the model behaves in higher dimensional cases. Algorithms such as IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 were used to generate the solution sets, however any other algorithms could be used with DoM indicator. The results have confirmed the effectiveness of the MIP DoM in problems with more than three objective functions. Final notes, considerations, and future research are discussed to exploit some solution sets particularities and improve the model and its use for other situations.

Keywords Multiobjective optimization multicriteria optimization quality indicators performance assessment exact method evolutionary algorithms

1 Introduction

Problems with conflicting objectives arise in most real world optimization cases. These problems are circumscribed by multi or many-objective optimization techniques [4] and exist in different domains such as [8], [5]. The solution sets for

these problems are formed in such a way that each objective represents a trade-off among other competing objectives.

In contrast to the mono-objective optimization, the solutions generated by multiobjective algorithms can be difficult to obtain and compare. Such difficulty grows as the number of candidate solutions and the objectives increase. When there are two or three objectives only, some graphical techniques help to examine the solution set visually. However, when the number of objectives is greater than three, this task is challenging, needing more advanced visualization techniques that can present location, shape, and solution set distribution [11].

Quality indicators are suitable for situations when we need to sum up the solution set [28], taking into account their characteristics. They have been used to compare the outcomes of multiobjective algorithms quantitatively (Section 2). In a recent paper [20], 100 quality indicators were discussed considering some state of the art indicators focusing on which quality aspects these indicators have, as well as its strengths and weaknesses.

In [19], a new quality measure, called dominance move (DoM) was proposed. This measure is able to capture all aspects of solution sets’ quality such as Pareto compliant, spread, uniformity, and cardinality (Section 3). DoM measures the minimum ‘effort’ that one solution set has to make in trying to dominate another set, more specifically the sum of the movement (i.e., Minkowski distance) needed to make a set dominant. The authors presented an exact approach to calculate DoM for the bi-objective case. Despite the low computational cost, the DoM method cannot be applied directly or extended to problems with three or more objectives (Section 4).

DoM presents good quality aspects for an indicator, but its calculation in low computational cost is a challenge [17], [18], and [20]. This work focuses on a DoM formulation based on a mixed integer programming approach (MIP)[14] aiming to overcome this difficulty. The mixed integer formulation is presented, and initial experiments showed that the MIP DoM formulation is correct. More specifically, this paper presents the following contributions: (i) a MIP model to DoM calculation that must be valid for all dimensions and solution sets; (ii) the use of this model to calculate DoM for some common problem sets having three, five, and ten objective functions, including considerations that come from the optimal DoM solution (Section 5); (iii) questions and details about the model behaviour for some test sets, raising future research paths to tackle some problems (Section 6).

2 Problems and deﬁnitions

In general, a multiobjective optimization problem (MOP) includes x decision variable vector from a decision space , and a set of M objective functions. Without loss of generality, a minimization MOP can be simply defined as [25]:

The is constructed by a set of M objective functions, which is a mapping from decision space to vectors in M -dimensional objective space . We are interested in the evaluation of these objective vector (solution) sets, and the comparison relation among them.

Considering two solutions , it is possible to establish a relation in which p is said to weakly dominate q if for 1 M, it is denoted as q. Moreover, if there is at least one objective i on which then it is said that p dominates q, and is denoted as . A solution is called Pareto optimal if there is no that dominates p. The set of all Pareto optimal solutions of an MOP is called Pareto optimal frontier. In the same way, the weak dominance relation can be stated to solution sets:

: The set P weakly dominates Q, i.e. it is denoted as , if every solution q Q is weakly dominated by at least one solution p P .

The comparison between two or more solution sets is of great importance. It can be used to compare the outcomes of multiobjective algorithms or even assist an algorithm during the search process for most proper candidates. It is paramount to make more precise statements in a quality indicator comparison, for example: if one algorithm is better than another, how much better is it? The following definition formalizes the quality indicators [28]:

Definition 2 Quality indicator: An k-ary quality indicator I is a function IR, which assigns each vector of k solution sets () a real value I(, , ..., ).

The quality indicators can be unary, binary, or k-ary, defining value to one solution set, two solution sets, or k solution sets, respectively. In [20], 100 indicators are listed and some are discussed in detail. Four indicator quality facets are analyzed: convergence, spread, uniformity, and cardinality. These facets will be analysed and discussed in the following.

It is expected that the Pareto dominance must be a central criterion in reflect-ing the convergence of solution sets. For two solution sets P and Q, for example, if P weakly dominates Q then I(P, Q) = 0. If P dominates some points of Q and Q does not dominate any point of P, it is reasonable to expect that the indicator prefer P to Q.

The spread of one solution set must consider the region that the set is covering. It involves both the outer portion and inner portion of the set. It must be able to capture the extensity of solution sets accurately.

The number of solutions in the set is another quality indicator facet named cardinality. Finally, a good indicator must prefer a set with uniformly distributed points, uniformity, showing an equidistant spacing amongst solutions.

It is plausible to add more details about the indicators due to its usage, one of them being the computational cost. Some indicators present all quality aspects but are hard to compute, notably in high dimensions or in high-cardinality sets. DoM and Hypervolume are well known examples. Other indicator details also deserve to be mentioned, such as the necessity for reference point or set, additional parameters, how to deal with scale, and normalization.

3 Related work

There are many indicators available, and they have been used in numerous situations in the literature [20]. Hypervolume (HV) [7], [24],[1], and inverted generational distance(IGD) [12] are some examples:

– Hypervolume (HV): Let = () be reference points in the objective space that is dominated by all solution sets. Let P be one solution set. The HV value of P with regard to represents the volume of the region which is dominated by P and dominates .

– Inverted generational distance (IGD): Let = () be a reference set of uniformly distributed points on the Pareto front. Considering P as an solution set to the Pareto front, the inverted generational distance between and P is defined as:

Note that d(r,P) is the minimum Euclidean distance from point r to solution set P. This indicator is a measure that represents how far the solution set is from the Pareto front reference. Lower values of IGD represent a better performance. The IGD metric is able to measure both diversity and convergence of P if is large enough [3].

There are other quality indicators such as Generational distance (GD) [12], -indicator [28], KKTPM [6], and others [20] .

An ideal quality indicator must present the four facets. Additionally, it must have a low computational cost, and it does not need a normalization (due to objective scale) and any additional parameters or reference points/sets.

Dominance move (DoM) is an intuitive indicator, and it has the four desirable facets [19]. The first idea presenting DoM came from the performance comparison indicator, PCI [18]. Examining the PCI proposal, it is quite similar to DoM in its essential purpose. PCI, a binary quality indicator, builds up a reference set using two solution sets, P and Q. This reference set is then split up into clusters, and the indicator calculates the movement distance in order to weakly dominate the clusters.

Dominance move is a measure for comparing two sets of multi-dimensional points being classified as a binary indicator. It considers the movement of points in one set to make this set weakly dominated by the other set. The DoM does not need a priori problem knowledge, parameters, or a reference set. However, the computational cost to calculate DoM is prohibitive, and the method to obtain DoM for a number of objective functions higher than two is unknown. DoM can be defined as follows.

Definition 3 DoM: Consider that P and Q are sets of points, with points and points . The dominance move of P to Q, DoM(P, Q), is the minimum total distance of moving points of P, such that any point in Q is weakly dominated by at least one point in P. In fact, the problem is to find from positions such that weakly dominates Q and the total move from to must be minimum [19].

is a set of points that are candidates to dominate Q with some update in one or many objectives that lead to a better distance such as expressed in Equation (3). In such way, it should be noted that with must be generated from p.

The formal expression of DoM can be stated as:

The dominance move indicator is based on the properties of dominance relation among solutions that are trying to dominate each other. These solutions’ efforts scale in a bottom up manner from the solutions to the set relations [19].

The number of possibilities to find is numerous. Any combination of some can dominate Q. Consider from P with some movement updates in the objectives in conjunction with other from original P, for example.

The number of such candidates is detailed in Equation (4), in which g is a group with one or many , and assuming that will be used as a base to be updated, generating that can weakly dominate all g group.

The Equation (4) is, in fact, an upper bound for the number of solutions. It is possible to discard some repeated generated from g and . It is worthy to note that the number of repeated candidates depends on the distributions for each solution sets in g and .

The original DoM formulation proposes a method for the biobjective case [20]. The algorithm can be outlined as:

Step 1: Remove the dominated points in both P and Q, separately. Remove the points of Q that are dominated by at least one point in P.

Step 2: Denote and start the process. Each point of Q in R is considered as a group. For each group of Q, find its inward neighbor ) in R. If the point , then merge r into the group of q, otherwise . If r is not assigned to one group merge the two groups of q and r into one group. Step 3: If there exists no point such that )) (i.e., there is a loop between the points) in any group, then the procedure ends and there is an optimal solution to the case.

Step 4: There is a loop in one or some groups. The procedure replaces these loops with the ideal point. The ideal point is formed of the best of each objective in each point inside the loop or group. Return to step 3 until convergence.

The authors present all definitions, theorems and corollaries to prove that this algorithm is correct in the biobjective case. However, it is stated that there is no method for three or more objectives.

4 The MIP dominance move calculation approach

Our goal is to modify the DoM formulation to deal with problems having three or more objective function. Our DoM calculation proposal is based on the perspective that the problem is, in fact, a particular instance of an assignment problem with two levels and some restrictions [21]. It is considered that to treat the problem, we have to find an assignment of P to Q with the restrictions that each must be assigned to one with the minimum distance. Nevertheless, in classic assignment problems, it is not possible that points change its own features, such as changing positions to alter the distances, for example. In DoM calculation, this issue must be considered.

A simple example to clarify the situation is given in [20]: consider P as {(2.0, 2.0, 2.0), (2.0, 2.2, 1.5), (3.0, 1.6, 1.6)} and Q as {(2.0, 1.2, 2.1), (2.0, 2.1, 1.0), (4.0, 1.5, 1.5)}. The inward neighbor ) of points and are, respectively, and . This create an assignment of P to Q with the minimum DoM(P, Q). Considering that P is fixed: DoM(P, Q) = ) + ) + ) = 1.6. However, if we considered a movement from P to , then would be transformed into = (2.0,1.5, 1.5). In this sense, we can find a better assignment and a lower value of DoM(P, Q) = 1.5. Clearly, other assignments from P to and to Q are capable to generate the same value (in [20], the authors presented another answer with DoM(P, Q) = 1.5.)

The problem can be modeled as a mixed integer programming approach. Generally, a MIP model (presented in Equation 5) can be described as a set of variables that are non negative, and variables that are integer [15].

Additionally, and are the objective function and left-hand side constraint coefficients, respectively. It is important to note that they are non negative. In the same way, there are the and for integer variables. Lastly, for the constraint set, there is a b right-hand side constant [15].

Using the MIP approach and DoM definition, P and Q are sets of points, with points and points and Q are given in the problem. Otherwise, is a set of points, with with , which are candidates generated from to weakly dominate some (i.e., or a g group), with some update in one or many objectives. It is relevant to note that if already dominates then = weakly dominates Q, resulting in a better distance such as expressed in Equation (3).

The proposal aims to calculate the distances ) and ) in the model and provide some limits for each that could vary due to g. Each g group

is formed by one or many ’s to be dominated.

In the model expressed in Equation (6) to Equation (12), these calculations can be done using MIP model. For each solution vector in P and , there are objectives. DoM was calculated using the auxiliary variables, that are equal to ), and that are equal to ).

The variable expressed in Equation (6), represents the distances between and a possible candidate. The constraints related to are shown in Equation (7). Essentially, it should be greater than zero and less than (always positive) and, finally, it must be greater than the difference in and .

The variable represents the difference for each in an attempt to dominate . The variable constraints are shown in Equation (8). It ensures that can dominate using a binary variable . Simultaneously, it will receive ), as the difference to dominate or 0. We have used a linearization technique to obtain the maximum function. Thus , which assumes 1 to guarantee the maximum value, was introduced; otherwise, the solution will be unfeasible. It is worth mentioning that there is a region in which can assume values, an upper and lower bound, previously defined.

The , a binary variable in Equation (12), was used to guarantee that if a is used in the model, we will know exactly which has generated the candidate. In the same way, the reflects that a generated a and is trying to dominate a . These two binary variables were used to guarantee that at least one is going to be dominated by . The con- straint represented in Equation (10) is derived from the conjunctive normal form as described in Equation (13).

In the same way, the constraints in Equation (11) were obtained from Equation (14). All the variables are binary and the constraints are expressed in Equation (12).

5 Experiments

5.1 Biobjective case

The first test was done to show how MIP DoM calculation addresses the quality indicator facets: convergence, spread, uniformity, and cardinality [20]. The same artificial experiments proposed in [19] to solve DoM in the biobjective case were adopted to verify the model correctness .

Figure 1 presents the artificial experiments for assessing the model correctness. Each facet of quality indicators, convergence, cardinality, uniformity, and spread, was represented in a ‘row’ in Figure 1. Each row was composed of two graphics corresponding to two examples. The examples are slightly different from one another. In each graphic there are two solution sets, P and Q.

Convergence is an important factor to reflect Pareto dominance compliance of sets. This behavior was shown considering two examples in Figure 1, first row of graphics: test = convergence. In Examples 1 and 2, the P and Q are equals, however, in Example 1, the P set was slightly improved in objective f1 , in comparison with Q. In Example 2, Q has some points dominated by P, but not all of them. The MIP DoM values of both graphics reflected the dominance relation.

MIP DoM prefers solutions with more cardinality. In Figure 1, test = cardinality, it could be observed that the solution sets had the same convergence, spread, and uniformity. In Example 1, Q has one point more than P, and in Example 2, P has two more points than Q.

Uniformity indicates the preference for evenly distributed points. The solution sets in Figure 1, test = uniformity, presented this feature. The sets have the same convergence, spread, and cardinality. Set P was distributed uniformly, and Q had a random distribution. In Example 1, Q was distributed in the range of set P, and in Example 2, the distance between neighboring points in set Q increased gradually from bottom to top.

Fig. 1 Artificial experiments proposed in [19] to assess the four facets of quality indicators: convergence, cardinality, uniformity, and spread. For each facet, there was a combination of two graphics. In total, there were four rows of graphics: tests for convergence, test for cardinality, test for uniformity, and test for spread. For each row, Examples 1 and 2 are given, and they are slightly different from one another. In each graphic there are two solution sets, P and Q.

Finally, MIP DoM must exhibit its preference for solutions with better spread. A set with better extensity has a smaller dominance move compared to its competitor. In Figure 1, test = spread, considering Example 1, set Q was generated by shrinking P a little (more concentrated in the middle). In Example 2, set P was distributed uniformly in the range, while Q assumed five bottom right points.

The MIP DoM approach results were equal to the results found in [19], using the biobjective algorithm proposed. The results agreement showed that the DoM MIP model is correct yielding to the same results obtained by its predecessor.

5.2 Multi and Many-objective case

Following the initial and artificial test sets, it would be crucial to validate the MIP DoM model with more than two objectives. Firstly, a test set with three objectives was done and compared with classical indicators, such as HV and IGD, and also with visual graphics to assess the results. Secondly, an attempt to solve problems with five and ten objective functions was executed. In all tests, algorithms such as IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 were used to generate the solution sets. It is important to highlight that the goal was to assess the effectiveness of the proposed MIP DoM model and not to compare algorithm performance, so any other algorithms could have been applied to generate the solution sets.

For each problem set, one of the parameters chosen for each algorithm was the population size. The question is closely related to Equation (4) and the cardinality of solution set, one of the quality indicator facet. To have a good approximation set of the Pareto front, in terms of convergence, spread, and uniformity, the number of non-dominated solutions grows exponentially in relation to the dimension of the objective space.

In [22], the shape of the Pareto front was discussed in the niche sizes definition. A useful limit to the number of individuals in the population, given by N = , was provided. The number of individuals in the population was N, M was the number of objectives, and r was the resolution required or the number of points needed to represent the Pareto front. This expression makes it clear how NP, or NQ, must be increased as the problem dimension grows, showing an exponential relation between N and M.

However, other works in many-objective optimization did not strictly follow this rule. In [23], for example, the number of objectives was M = {3, 5, 10} and the population size was N = {105, 126,275}. In [24], an efficient hypervolume calculation was provided and some tests were done with M = {3, 4, 5} and N = {10, .., 200}. Likewise, in CEC’2018, a competition on many-objective optimization [3] established M = {5, 10, 15} and the maximum population size as 240.

Based on Equation (4), which was important in the proposed model (Equations (6) to (12)), we have decided to validate DoM using M = {3, 5, 10} and NP = NQ = {50, 100, 200} indicating the final Pareto front size.

All the experiments were done using Platypus [2] and PyGMO [13] to generate the problem sets and to calculate the indicators (HV and IGD). The number of fitness evaluations was the same for all algorithms, 10000. The model in Equation (6) to (12) was implemented using Python and GUROBI [9] version 9.0.0 build v9.0.0rc2 running on a Linux 64 bits operational system with 12 CPU’s and 16Gb of RAM.

In [10], a review of multiobjective problem sets was presented. To validate our proposed approach, the aim was to select some well-known problem test sets in three dimensions, such as the DTLZ and WFG families. Some test sets were initially chosen based on the characteristics: convexity/concavity, disconnection, multimodality, and degeneracy. DTLZ1, DTLZ2, DTLZ3, and DTLZ7, which are linear, concave/multimodality, concave, and disconnected, respectively, were selected. In the same manner, WFG1, WFG2, WFG3 and WFG9 were also selected taking into account the properties: convex/mixed, convex/disconnected, linear/degenerate, and concave, respectively.

In parallel to the problem set, some algorithms were used to test and compare the MIP DoM with other indicators. Embedded with a similar purpose of that in[16], the NSGAIII and MOEAD algorithms (Pareto-based and decomposition approach, respectively) were firstly chosen, and afterwards, IBEA [26], SPEA2 [27], and NSGAII [1]. All the problem sets tested can be seen in Figure 2.

Fig. 2 Solution sets with NP = NQ = 50 solutions and three objectives, M = 3, generated by IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 algorithms applied to DTLZ1, DTLZ2, DTLZ3, WFG1, WFG2, WFG3 and WFG9 problem sets.

Tables 1 and 2 show two unary quality indicators, the inverted generational distance (IGD) and hypervolume (HV). It is mandatory to have reference sets to calculate these indicators, and this task is not only a challenge one [20] but also sometimes is provided by the user [24]. We decided to use the maximum and minimum values, amongst all algorithm solutions and objectives, for HV and IGD.

It is noted that DTLZ1 and DTLZ3, for example, had the HV value inflated by the presence of the dominance resistant solutions [20], non-dominated solutions with a poor value in one objective but with good values in others. It is also important to observe that IGD and HV are highly sensitive to the reference set. If such reference set does not represent well the Pareto front, then the indicator is compromised. This parameter is crucial to these indicators’ validity.

Table 1 IGD quality indicator for problem sets generated using DTLZ1, DTLZ2, DTLZ3, DTLZ7, WFG1, WFG2, WFG3 and WFG9 for IBEA, MOEAD, NSGAIII, NSGAII and SPEA.

The MIP DoM quality measure is a binary indicator. In Table 3, for example, it is possible to see all the comparisons among the algorithms for the DTLZ family problem set. P was the solution set generated by the algorithm which was trying to dominate, and Q was the solution set which was being dominated.

Using MIP DoM as a unary indicator is still feasible, a straightforward idea was to merely apply a summation by the values for each algorithm in the Table 3 (summing up the elements of a row, for example) .

Table 2 HV quality indicator for problem sets generated using DTLZ1, DTLZ2, DTLZ3, DTLZ7, WFG1, WFG2, WFG3 and WFG9 for IBEA, MOEAD, NSGAIII, NSGAII and SPEA.

For the DTLZ1 test set, which was detailed in Tables 1 and 2, algorithms that present the best IGD were IBEA and NSGAIII. For the HV indicator, it was difficult to compare the algorithms due to inflated values. Using MIP DoM approach and the comparison among algorithms, the IBEA and NSGAIII were both indicated as the best solutions. The results are presented in Table 3. For example, if the MIP DoM values were summed up, the IBEA shows a 2.424 and NSGAIII 4.175. Moreover, when IBEA tried to dominate NSGAIII, it presented a lower value, 1.078, compared to 1.359 when NSGAIII tried to dominate IBEA.

In the DTLZ2 case, IGD indicated IBEA, MOEAD, and NSGAIII algorithms (see Table 1). Considering the HV, IBEA algorithm was the best one (Table 2). Again, looking at Table 3, the best values also indicated IBEA, MOEAD, and NSGAIII algorithms. However, when making two by two comparisons, MOEAD presented better value than IBEA and NSGAIII.

Table 3 MIP DoM values for the problem set of the DTLZ family for comparison among IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 algorithms. It must be noted that P was the solution set generated by the algorithm that was trying to dominate, and Q was the solution set generated by algorithm being dominated. Each solution set had NP = NQ = 50 solutions and M = 3 dimensions.

The result for DTLZ3 was presented in Tables 1 and 2: for IGD, the best algorithms were NSGAII and MOEAD; and IBEA, MOEAD, and NSGAII had the best HV values. Summing up the MIP DoM values among all algorithms (as shown in Table 3), the best algorithms were IBEA and MOEAD. In this case,MIP DoM was in agreement with HV. Comparing IBEA and MOEAD directly, MOEAD presented a lower DoM.

In DTLZ7 problem set, the best HV values were NSGAIII and NSGAII. Considering IGD, the best values were for MOEAD and SPEA2. Applying the same summing up approach as before, and using Table 3, the result was that NSGAIII and NSGAII generated the best candidate solutions. Again, comparing NSGAIII and NSGAII, the best one was NSGAIII with a DoM(P, Q) = 1.012.

The same experiment was done for the WFG family. For WFG1 in Tables 1 and 2, the best algorithm was IBEA for IGD, and HV indicator. In Table 4, the best algorithm was SPEA2, and the second one was IBEA. However, comparing SPEA2 with IBEA, SPEA2 had a lower value of MIP DoM, DoM(SPEA2, IBEA) = 0.246, in contrast with DoM(IBEA, SPEA2) = 0.916. This behavior could be

Table 4 MIP DoM values for the problem set of the WFG family for comparison among IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 algorithms. It must be noted that P was the solution set generated by the algorithm that was trying to dominate, and Q was the solution set generated by algorithm being dominated. Each solution set had NP = NQ = 50 solutions and M = 3 dimensions.

observed in Figure 2. Visually, it is possible to see a similarity between the two solution sets. However, it is relevant to note the graphics’ scale in the axis. The SPEA2 solution set appeared more convergent than IBEA (the spread are similar, but SPEA2 presented minor values in the axis).

In WFG2, the best values for IGD were entirely tied, with a little difference (see Table 1). Considering HV, the best one was the IBEA algorithm (Table 2). Using DoM, detailed in Table 4, it was possible to see the same characteristic as reported by IGD. Nonetheless, if the algorithms were compared two by two, it was possible to indicate IBEA as the algorithm which presented the lower MIP DoM values.

Using WFG3, it was possible to observe little differences for the IGD indicator among the other algorithms. For HV, there was an indication that IBEA had a greater value, but with other values next to it. Assessing MIP DoM in Table 4, there was an indication that SPEA2 had better values. Furthermore, it was still possible to analyze the sets graphically in Figure 2. A visual and qualitative analysis comparing IBEA with SPEA2 graphics could show how DoM acknowledged the spread’s solution set. All the solution sets had a similar spread and cardinality (it is possible to see a kind of straight line formed by points), and SPEA2 was one notable exception.

Finally, for the WFG9 problem set, Tables 1 and 2 showed that for IGD, algorithms IBEA and MOEAD had lower values, but the values were next to each other. For HV, the IBEA algorithm had a better value. Considering DoM, presented in Table 4, all three algorithms were competitive: IBEA, MOEAD, and NSGAIII. Comparing IBEA and MOEAD, for example, the best MIP DoM value was found for IBEA: DoM(IBEA,MOEAD) = 1.932.

The discrepancy for the WFG1 and WFG3 experiments, among DoM, IGD, and HV was something that drew attention. For DoM, the optimal value was obtained using the model described in 6. However, the reference sets for IGD and HV were generated in a simple manner. This might be improved in assessing the convergence among the indicators. Another possibility would be to increase the number of solutions in each set.

DTLZ1 DTLZ2 DTLZ3 DTLZ7 WFG1 WFG2 WFG3 WFG9 Problem set

Fig. 3 A box plot with the time spent by the GUROBI solver [9] to solve the MIP DoM model for the DTLZ1, DTLZ2, DTLZ3, DTLZ7, WFG1, WFG2, WFG3, and WFG9 problem sets with NP = NQ = 50 points and M = 3 objectives.

The GUROBI solver [9] used branch and bound, branch and cut, and other solver capabilities, such as: Gomory cuts, flow cover, clique, and others. The minimum gap established was five percent, and the only additional parameter used was MIPFocus, which helps to improve the best bound during the execution. For some models, the time spent was small. Figure 3 shows a box plot for the results with the medians highlighted in bold. All the models were solved, and the median was less than 35 seconds. The only exception was DTLZ2, which presented the highest median time spent, 1289. Furthermore, there were two cases in high discrepancy, more specifically: IBEA trying to dominate NSGAII, with 89444 seconds, and NSGAIII trying to dominate NSGAII with 35325 seconds. These two cases had some common characteristics of what the model had spent more time than others (high density between the solution sets, for example).

In Figure 4, a piece of similar information was depicted: the number of simplex iterations by the solver. Again, in some cases, such as DTLZ3 and WFG1, the number of iterations was small. The same discrepancy cases, as shown in Figure 3, happened here with DTLZ2.

In this Section, the goal was to verify if the MIP DoM approach could be applied in many-objective scenario. The same problem sets and algorithms from the last experiment were used. It worthy to note that the number of points was established based on some works ([23] and [24], for example).

DTLZ1 DTLZ2 DTLZ3 DTLZ7 WFG1 WFG2 WFG3 WFG9 Problem set

84070 12991934 2 119821 27 155802 200838 155836

Fig. 4 A box plot with the simplex iterations used by the GUROBI solver [9] to solve the DoM model for the DTLZ1, DTLZ2, DTLZ3, DTLZ7, WFG1, WFG2, WFG3, and WFG9 problem sets with NP = NQ = 50 points and M = 3 objectives.

In Table 5, problem sets from the DTLZ family are presented. At this time, the problems had M = 5 objectives and NP = NQ = 100 solutions. The non-dominated points were generated using IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2. For the DTLZ1 problem set, the best algorithms were IBEA and SPEA2, but when they were compared, DoM(IBEA, SPEA2) = 0.525 and DoM(SPEA2, IBEA) = 2.278, IBEA was still better than SPEA2.

Following this idea, for DTLZ2, the IBEA and SPEA2 were more competitive than the others. In this case, two by two comparison using DoM was very tight (1.000 and 1.060, for example). However, IBEA presented a lower value. Considering DTLZ3, again, IBEA and SPEA2 were the best ones, and IBEA exhibited the best value, compared with SPEA2 or any other algorithm. Finally, for DTLZ7, it was possible to observe a curious fact. In this case, SPEA2 already weakly dominated all other solution sets.

In Table 6, problems from the WFG family were tested, with the same proposal as before. For WFG1, the DoM value comparison was tight for IBEA, MOEAD, NSGAIII, and SPEA2. The ranking, using the idea of move dominance summation, was for NSGAIII and MOEAD as the best ones.

For the WFG2 problem, all the lower values reference IBEA algorithm as the best solution set for this problem. Considering WFG3, there was, again, a curious fact: the solution set generated by IBEA almost dominated all other results, the only exception was for MOEAD. Finally, for WFG9, it was clear that the solution set generated by IBEA was lower than the others.

Another relevant fact which is important to observe: some MIP DoM values were repeated, when a solution P was trying to dominate Q, for example. See in Table 6, the problem set WFG1 and P as NSGAII, or the problem set WFG2 and P as MOEAD, or even the problem set WFG9 and P as SPEA2. These cases raised one question: What are the points in P that were updated to dominate Q? In all the last cases, the same points were generated by the same P. For example, for the case with the solution set WFG1 and P as NSGAII, was always the point where was generated, and this was the final solution for all the cases (i.e., in each case the values updated in each objective were different, but all the final solutions used ). The same fact happened with for all last cited cases.

In Figure 5, the number of simplex iterations was presented. There was a discrepancy in DTLZ2, which demanded much more iterations than the others.

Table 5 MIP DoM values for the problem set of the DTLZ family for comparison among IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 algorithms. It must be noted that P was the solution set generated by the algorithm that was trying to dominate, and Q was the solution set generated by algorithm being dominated. Each solution set had NP = NQ = 100 solutions and M = 5 dimensions.

This was when IBEA was trying to dominate SPEA2. The time spent in this extreme case was 10245 seconds.

DTLZ1 DTLZ2 DTLZ3 DTLZ7 WFG1 WFG2 WFG3 WFG9 Problem set

31024 735716 32 364452 448 534245 1899487 1262213

Fig. 5 A box plot with the simplex iterations used by the GUROBI solver [9] to solve the DoM model for the DTLZ1, DTLZ2, DTLZ3, DTLZ7, WFG1, WFG2, WFG3, and WFG9 problem sets with NP = NQ = 100 points and M = 5 objectives.

Table 6 MIP DoM values for the problem set of the WFG family for comparison among IBEA, MOEAD, NSGAIII, NSGAII, and SPEA2 algorithms. It must be noted that P was the solution set generated by the algorithm that was trying to dominate, and Q was the solution set generated by algorithm being dominated. Each solution set had NP = NQ = 100 solutions and M = 5 dimensions.

The final experiment was made to verify if the MIP DoM approach was viable in a higher dimensional objective space . The test was done considering M = 10 objectives and NP = NQ = 200. It was necessary to reduce the number of algorithms. In Platypus [2], there was a parameter to specify the population size; however, for MOEAD and NSGAIII, this parameter was ignored for ten dimensions. For this reason, in Table 7, these algorithms were not presented.

Table 7 shows that the MIP DoM approach could be computed using these test sets. There were some cases in which the model spent 2.05 seconds, with 2556 simplex iterations for one case, or 29.97 seconds, with 45674 simplex iterations, for other cases. These two examples were the fastest ones. On the other hand, there was another one, in which 180468.60 seconds was spent, with 107009076 simplex iterations, presenting the worst case.

For the fastest cases, there were many points in the solution set Q, that already were dominated by some point in P. For these cases, the model and the solver exploited this characteristic, and they were capable of finding the best solution in a reasonable time.

On the other hand, the worst case spent almost 51 hours. In this situation, there was no a priori dominance among the points in each solution set. One thing that can be better investigated is the density between the solution sets when there is a chance for p to change its position to dominate some , considering a g group. In areas with high density, it is plausible that there are many options concerning P and Q, and the branch and bound dynamic become more challenging.

One last interesting fact observed in the solutions generated by the MIP DoM model is related to the number of that was altered to obtain the MIP DoM value. For this final experiment, for example, in which the space was large, the maximum number of that had their objectives updated was six. The minimum number of was one. For these last cases, the MIP DoM recommendation was to change only one point, and that point was capable of dominating all the solutions and generating the best MIP DoM value.

Table 7 MIP DoM values for DTLZ and WFG family problem set for comparisons among IBEA, NSGAII, and SPEA2 algorithms. It must be noted that P was the solution set generated by the algorithm that was trying to dominate, and Q was the solution set generated by algorithm being dominated. Each solution set had NP = NQ = 200 solutions and M = 10 dimensions.

As can be viewed in Equation (4), the combinatorial nature of the MIP DoM calculation is one reason which has made it difficult to use. The computational time increases exponentially with the size of the solution set. The mixed integer programming is NP hard in general [14]. Despite this fact, the MIP approach has managed to be viable in three objectives with 50 solutions, five objectives with 100 solutions, and ten objectives with 200 solutions for some DTLZ and WFG family test cases.

Concerning the MIP model, some issues need to be mentioned:

i) In the LP relaxation for and xpq(i, j) (in the first two constraints in Equation (12)), the integrality can not be ensured; this problem is related in [15] causing, in some cases, a slow advance in the best bound;

ii) The last problem could be solved if we knew, a priori, the number of that would dominate Q. For the last experiment, for example, it was observed that in a high number of dimensions, a low percentage of points was used to dominate the other set. We believe that imposing a constraint that limits the number of can improve the number of simplex iterations.

iii) The coefficient matrix to bound presented a large amplitude in some cases, and it causes numerical difficult to solvers.

6 Conclusion

The DoM binary indicator considers the minimum move of one set to weakly dominate the other set. Given two solution sets P and Q, the DoM of P to Q is the minimum total distance of moving some points in P in way that any point in Q is weakly dominated by at least one point in P. The indicator is Pareto compliant and does not demand any parameters or reference sets. The handicap about DoM is its combinatorial calculation nature making it difficult to be applied to problems with three or more objective functions.

We explored a MIP approach calculation of DoM, trying to understand if it was a viable way to deal with the problem. The idea was to use the polyhedral method, as it has been proved to be successful in practice with other issues. A comparison with artificial bidimensional examples, the same proposed by the DoM’s authors was made, with the same results. Additionally, some problems in three dimensions showed that MIP DoM results were in agreement and compliant compared with other common indicators used in literature (IGD and HV). Finally, experiments using five and ten dimensions were executed and had shown that MIP DoM model was viable and, in some cases with reasonable computational time. All the experiments have used a fixed number of solutions in each set (i.e., NP = NQ = {50,100, 200}, for M = {3, 5, 10} dimensions). To the best of our knowledge, even with these limitations, an exact method to calculate DoM in three or more dimensions was not known until now.

For some experiments, high variability was found in the time spent by the model resolution. It was observed that this fact was inherent to the distribution or density, involving the two solution sets, P and Q. This issue is something that deserves to be best investigated in order to improve the model, creating, for example, new constraints in the model.

In the MIP approach for DoM, other directions also deserve to be investigated: i) How to efficiently calculate DoM using MIP and exploiting some inherent solution set features; ii) How to define a priori the minimum number of to be moved in an attempt to dominate Q and what benefits it can bring to the MIP model. What the relation between the distributions of P and Q, and the minimum number of to be moved is; iii) How to improve the amplitude of the coefficients’ matrix in the MIP model, and if this is this capable of improving the convergence model time. Such questions were posed here to emphasize possible future research directions. The tests with more dimensions and experiments with other problems test sets are something that must be done, as well.

Finally, MIP DoM is an indicator that can represent all the solution sets quality facets. Moreover, it proposes a more natural and intuitive relation when comparing solution sets using a measure compatible with the Pareto dominance concept, mainly in high dimensional objective space. The MIP DoM, like its predecessor, does not require prior problem knowledge and any other additional parameter. With MIP DoM Model, a viable method to calculate the DoM indicator was proposed to try to incentive its general use for compare solution sets or even to guide some algorithms in the search for solutions with useful quality features.

References

1. Bradford, E., Schweidtmann, A., Lapkin, A.: Efficient multiobjective optimization employing gaussian processes, spectral sampling and a genetic algorithm. Journal of Global Optimization 71 (2018). DOI 10.1007/s10898-018-0609-2

2. Brockhoff, D., Tuˇsar, T.: Benchmarking algorithms from the platypus framework on the biobjective bbob-biobj testbed. In: Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO ’19, pp. 1905– 1911. ACM, New York, NY, USA (2019). DOI 10.1145/3319619.3326896. URL http://doi.acm.org/10.1145/3319619.3326896

3. Cheng, R., Li, M., Tian, Y., Xiang, X., Zhang, X., Yang, S., Jin, Y., Yao, X.: Benchmark functions for the cec’2018 competition on many-objective optimization. Tech. rep. (2018)

4. Chugh, T., Sindhya, K., Hakanen, J., Miettinen, K.: A survey on handling com- putationally expensive multiobjective optimization problems with evolutionary algorithms. Soft Comput. 23(9), 3137–3166 (2019). DOI 10.1007/s00500-017-2965-0. URL https://doi.org/10.1007/s00500-017-2965-0

5. Deb, K.: Scope of stationary multi-objective evolutionary optimization: A case study on a hydro-thermal power dispatch problem. Journal of Global Optimization 41, 479–515 (2008). DOI 10.1007/s10898-007-9261-y

6. Deb, K., Abouhawwash, M.: An optimality theory-based proximity measure for set-based multiobjective optimization. IEEE Transactions on Evolutionary Computation 20(4), 515–528 (2016). DOI 10.1109/TEVC.2015.2483590

7. Deng, J., Zhang, Q.: Approximating hypervolume and hypervolume contributions using polar coordinate. IEEE Transactions on Evolutionary Computation 23(5), 913–918 (2019). DOI 10.1109/TEVC.2019.2895108

8. Gonz´alez-´Alvarez, D.L., Vega-Rodr´ıguez, M.A.: Analysing the scalability of multiobjective evolutionary algorithms when solving the motif discovery problem. Journal of Global Optimization 57, 467–497 (2013)

9. Gurobi Optimization, L.: Gurobi optimizer reference manual (2019). URL http://www.gurobi.com

10. Huband, S., Hingston, P., Barone, L.: A Review of Multi-objective Test Problems and a Scalable Test Problem Toolkit A Review of Multiobjective Test Problems and a Scalable Test Problem Toolkit. IEEE Transactions on Evolutionary Computation 10(2006), 477– 506 (2011). URL http://ro.ecu.edu.au/ecuworks/2022

11. Ibrahim, A., Rahnamayan, S., Martin, M.V., Deb, K.: 3d-radvis antenna: Visualization and performance measure for many-objective optimization. Swarm and Evolutionary Computation 39, 157–176 (2018). DOI 10.1016/j.swevo.2017.09.011. URL https://doi.org/10.1016/j.swevo.2017.09.011

12. Ishibuchi, H., Masuda, H., Tanigaki, Y., Nojima, Y.: Modified distance calculation in generational distance and inverted generational distance. In: A. Gaspar-Cunha, C. Henggeler Antunes, C.C. Coello (eds.) Evolutionary Multi-Criterion Optimization, pp. 110–125. Springer International Publishing, Cham (2015)

13. Izzo, D.: Pygmo and pykep: Open source tools for massively parallel optimization in as- trodynamics (the case of interplanetary trajectory optimization). pp. – (2012)

14. J¨unger, M., Liebling, T.M., Naddef, D., Nemhauser, G.L., Pulleyblank, W.R., Reinelt, G., Rinaldi, G., Wolsey, L.A. (eds.): 50 Years of Integer Programming 1958-2008 - From the Early Years to the State-of-the-Art. Springer (2010). DOI 10.1007/978-3-540-68279-0. URL https://doi.org/10.1007/978-3-540-68279-0

15. Klotz, E., Newman, A.: Practical guidelines for solving difficult mixed integer linear pro- grams. Surveys in Operations Research and Management Science 18, 18–32 (2013). DOI

10.1016/j.sorms.2012.12.001

16. Li, H., Deb, K., Zhang, Q., Suganthan, P., Chen, L.: Comparison between moea/d and nsga-iii on a set of many and multi-objective benchmark problems with challenging difficulties. Swarm and Evolutionary Computation 46, 104–117 (2019). DOI

10.1016/j.swevo.2019.02.003

17. Li, M.: Evolutionary many-objective optimisation: Pushing the boundaries. Ph.D. thesis, Department of Computer Science, Brunel University London, London, UB8 3PH, U.K. (2015)

18. Li, M., Yang, S., Liu, X.: A performance comparison indicator for pareto front approximations in many-objective optimization. In: Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation, GECCO ’15, pp. 703–

710. ACM, New York, NY, USA (2015). DOI 10.1145/2739480.2754687. URL http://doi.acm.org/10.1145/2739480.2754687

19. Li, M., Yao, X.: Dominance move: A measure of comparing solution sets in multiobjective optimization. CoRR abs/1702.00477 (2017). URL http://arxiv.org/abs/1702.00477

20. Li, M., Yao, X.: Quality evaluation of solution sets in multiobjective optimisation: A survey. ACM Comput. Surv. 52(2), 26:1–26:38 (2019). DOI 10.1145/3300148. URL http://doi.acm.org/10.1145/3300148

21. Pentico, D.W.: Assignment problems: A golden anniversary survey. European Journal of Operational Research 176(2), 774–793 (2007). URL https://ideas.repec.org/a/eee/ejores/v176y2007i2p774-793.html

22. Sen, P., Yang, J.B.: Multiple criteria decision support in engineering design (1998)

23. Tian, Y., Cheng, R., Zhang, X., Cheng, F., Jin, Y.: An indicator-based multiobjective evolutionary algorithm with reference point adaptation for better versatility. IEEE Transactions on Evolutionary Computation 22(4), 609–622 (2018). DOI 10.1109/TEVC.2017. 2749619

24. Yang, K., Emmerich, M., Deutz, A., B¨ack, T.: Efficient computation of expected hypervolume improvement using box decomposition algorithms. Journal of Global Optimization 75(1), 3–34 (2019). DOI 10.1007/s10898-019-00798-7. URL https://doi.org/10.1007/s10898-019-00798-7

25. Yuan, Y., Ong, Y.S., Gupta, A., Xu, H.: Objective Reduction in Many-Objective Opti- mization: Evolutionary Multiobjective Approaches and Comprehensive Analysis. IEEE Transactions on Evolutionary Computation 22(2), 189–210 (2018). DOI 10.1109/TEVC. 2017.2672668

26. Zitzler, E., K¨unzli, S.: Indicator-based selection in multiobjective search. In: in Proc. 8th International Conference on Parallel Problem Solving from Nature (PPSN VIII, pp. 832–842. Springer (2004)

27. Zitzler, E., Laumanns, M., Thiele, L.: Spea2: Improving the strength pareto evolutionary algorithm. Tech. rep. (2001)

28. Zitzler, E., Thiele, L., Laumanns, M., Fonseca, C.M., da Fonseca, V.G.: Performance assessment of multiobjective optimizers: an analysis and review. IEEE Transactions on Evolutionary Computation 7(2), 117–132 (2003). DOI 10.1109/TEVC.2003.810758

designed for accessibility and to further open science