Objective space division-based hybrid evolutionary algorithm for handing overlapping solutions in combinatorial problems

1.   Introduction

A multi-objective optimization problem (MOP) can be defined as a problem corresponding to a certain reality, where a decision-maker wishes to optimize several objectives simultaneously, usually in conflict with each other. When solving a MOP of hard complexity, metaheuristic methods are very appropriate. These methods do not guarantee obtaining the exact front of non-dominated solutions, but an approximate set.

The second generation of Multi-Objective Evolutionary Algorithms (MOEAs) has been shown to achieve excellent results solving MOP ^[1,2,3,4]. Stochastic in nature, MOEAs are based on the concept of population of solutions, which gives them great skill in finding multiple non-dominated solutions in solution spaces of diverse nature. The Non-dominated Sorting Genetic Algorithm-II (NSGA-II) ^[5], the Multi-Objective Evolutionary Algorithm based on Decomposition (MOEA/D) ^[6] and the Global Weighting Achievement Scalarizing Function Genetic Algorithm (Global WASF-GA) ^[7], among others, are well-known in the literature. NSGA-II is based mainly on two mechanisms: Pareto dominance as a criterion to converge to the Pareto Optimal Front (POF) and crowding-distance operator as augmenting diversification in the population. MOEA/D explicitly decomposes the MOP into a number of scalar optimization sub-problems which are solved simultaneously through the evolution of a population of solutions. Global WASF-GA, similarly to NSGA-II, uses Pareto dominance as a criterion to converge to the POF but based on the Achievement Scalarizing Function (ASF) proposed by Wierzbicki ^[8]. However, the diversity is enhanced by the use of a set of well-spread projection directions, and by not allowing repeating solutions in the same front. Both MOEA/D and Global WASF-GA use a scalarizing function as fitness function and simultaneously employ several predefined weight vectors for searching new solutions (aggregation methods). Moreover, Global WASF-GA uses two reference points (the nadir and the utopian points) simultaneously in the scalarizing function and most versions of MOEA/D just consider one reference point (the ideal point). Some papers discuss the reference point specification in MOEA/D ^[7,9,10,11].

In MOEAs, diversity maintenance plays an important role in pushing the population towards the POF, but also in avoiding overlapping solutions. In ^[12] it is shown that there is a large number of overlapping objective vectors in each population when MOEAs are applied to multi-objective combinatorial optimization problems with few objectives while this is not the case in the application to MOPs with continuous decision variables and/or many objectives.

In the past, but more recently, papers have been published that apply the Objective Space Division (OSD) to improve MOEAs in the sense of increasing diversity and generate well-distributed solutions along the entire Pareto front ^{[13,14,15,16,17,18,19]}. All of them applied to test problems with continuous decision variables and/or three or more objectives.

The Multi-Objective Knapsack Problem (0-1 MOKP) is a well-known and widely studied combinatorial optimization problem, which has been attempted to be solved by a considerable number of exact and metaheuristic methods ^{[6,20,21,22,23,24,25,26]}. It should be noted that the appearance of overlapping solutions when the bi-objective 0-1 MOKP is solved with a MOEA has a significant negative impact on the diversity of the final approximate front of non-dominated solutions found ^[12].

A hybrid of NSGA-II is proposed in this paper, which divides the objective space into several regions using the nadir solution that is updated at each generation of the algorithm. In addition, a different objective function optimization strategy is used in each region. The underlying idea is to increase the diversity of solutions of the approximate POF reached, allowing that, in each generation of the algorithm, solutions dominated by the extremes of the POF reached can be included in the population of non-dominated ones. Since the proposed algorithm is built from NSGA-II, it cannot be inferred that the proposed methodology is a general solution for EAs, but only for those that have a similar behavior to NSGA-II: concentrated approximate Pareto front with solutions far from the extreme values of the true Pareto front. Preliminary results were presented at the International Conference of Production Research, Americas 2020 ^[27].

This article is structured as follows. The next section briefly explains some multi-objective basic concepts that will make easier to understand the work presented here. Section 3 details the proposed methodology. Section 4 details the experiments carried out to determine the parameters to be used in Section 5 where the results obtained with NSGA-II, MOEA/D, Global WASF-GA and the proposed algorithm are compared. Lastly, Section 6 contains the conclusions.

2.   Basic concepts

In terms of maximization, a general MOP can be defined as follows:

where $n\ge 2$ is the number of objective functions to be maximized simultaneously, $p$ and $q$ are the number of equations and inequations, respectively, that define the feasible space $S\subseteq {\mathbf{R}}^{m}$. Thus, $\mathbf{x} = {\left({x}_{1}, {x}_{2}, \dots , {x}_{m}\right)}^{T}\in S$ is a decision vector of $m$ variables.

The set of images of each achievable or feasible solution in the decision space forms the set of achievable or feasible solutions in the objective space $Z = f\left(S\right)$ defined as:

A solution ${\mathbf{z}}^{t} = {\left({z}_{1}^{t}, {z}_{2}^{t}, \dots , {z}_{n}^{t}\right)}^{T}\in Z$ dominate a solution ${\mathbf{z}}^{u} = {\left({z}_{1}^{u}, {z}_{2}^{u}, \dots , {z}_{n}^{u}\right)}^{T}\in Z$ if the following conditions are satisfied:

The set of all non-dominated solutions $z\in Z$ in the objective space is referred to as Pareto Optimal Front (POF) (see Figure 1).

The ideal solution ${z}^{I}$ of a MOP is a solution vector ${\mathbf{z}}^{I} = {\left({z}_{1}^{I}, {z}_{2}^{I}, \dots , {z}_{n}^{I}\right)}^{T}$, constructed with the best objective function values of all the solutions of the objective space (see Figure 1). This solution is generally unfeasible or unattainable.

Meanwhile, the nadir solution ${z}^{N}$ of a MOP is a solution vector ${\mathbf{z}}^{N} = {\left({z}_{1}^{N}, {z}_{2}^{N}, \dots , {z}_{n}^{N}\right)}^{T}$, constructed with the worst objective function values of all the solutions of the Pareto front (see Figure 1). This solution can be feasible or unfeasible.

2.1.   MOEA/D with Zhang's specification of the reference point

Basically, MOEA/D explicitly decomposes the MOP into a number of scalar (mono-objective) optimization sub-problems which are solved simultaneously through the evolution of a population of solutions. Some approaches for reducing a MOP into a set of scalar sub-problems are discussed in ^[6]. One such approach is that of Tchebycheff, which can be formulated as follows:

where ${\mathbf{z}}^{*} = {\left({z}_{1}^{*}, {z}_{2}^{*}, \dots , {z}_{n}^{*}\right)}^{T}$ is the reference point that satisfy, for $\alpha \ge 1$:

This version of MOEA/D will be referred to as MOEA/D Zhang^te or simply Zhang^te.

2.2.   MOEA/D with Ishibuchi's specification of the reference point

For the maximization problem, Ishibuchi et al. ^[9] proposed some innovations to MOEA/D. Namely, they propose to replace Eq (2.2) with a dynamic strategy for specifying the reference point. This strategy is defined as follows:

where ${\mathbf{z}}^{*} = {\left({z}_{1}^{*}, {z}_{2}^{*}, \dots , {z}_{n}^{*}\right)}^{T}$ and ${\mathbf{z}}^{min} = {\left({z}_{1}^{min}, {z}_{2}^{min}, \dots , {z}_{n}^{min}\right)}^{T}$, being ${z}_{i}^{*}$ and ${z}_{i}^{min}$, respectively, the maximum and minimum values reached by the objective function $i$ in the current generation $t$, and ${\alpha }_{t}$ is a parameter that decreases in each generation of the algorithm. It follows from Eq (2.4) that ${\mathbf{z}}^{R} = {\mathbf{z}}^{*}$ when ${\alpha }_{t} = 0$.

Also, in order to make the ${\mathbf{z}}^{R}$ value gradually approaches the value of the reference point ${\mathbf{z}}^{*}$ during the execution of the algorithm, the ${\alpha }_{t}$ value is modified at each generation according to:

where ${\alpha }_{0}$ is the initial value of $\alpha$, ${t}_{max}$ is the maximum number of generations (of the algorithm) and $t$ is the current generation. Note that the ${\alpha }_{t}$ value in the last generation is zero, i.e., ${\mathbf{z}}^{R}$ tends to ${\mathbf{z}}^{*}$ during the execution of the algorithm.

This version of MOEA/D will be referred to as MOEA/D Ishibuchi^te or simply Ishibuchi^te.

2.3.   NSGA-II

Of recognized efficiency, the general operation of NSGA-II ^[5] is as follows:

1) Define the initial parameters: $N$ (population size), ${t}_{max}$ (maximum number of generations), Pcr (crossover probability), Pmut (mutation probability).

2) Randomly generate an initial population ${P}^{0}$ of size $N$.

3) Assign $t = 0$ and ${P}^{t} = {P}^{0}$.

4) Apply selection by tournament (best rank and crowding) and crossover and mutation operators to create a population of descendants ${Q}^{t}$ of size $N$.

5) Form a combined population ${R}^{t} = {P}^{t}\cup {Q}^{t}$.

6) Classify ${R}^{t}$ on non-dominated fronts: ${\left\{{F}_{j}^{t}\right\}}_{j\ge 1}$.

7) Add each front, in increasing order until a front ${F}_{k}^{t}$ cannot be fully integrated into the new population ${P}^{t+1}$ of size $N$, then classifying the ${F}_{k}^{t}$ front by crowding-distance to add solutions until the size of ${P}^{t+1}$ equals $N$.

8) Assign $t = t+1$ and ${P}^{t} = {P}^{t+1}$.

9) If the stop criterion is not met, go to step 4.

3.   Proposed method: NSGA-II/OSD

The classical NSGA-II -presented in the previous section- fails to approximate the whole POF when addressing combinatorial problems with a low number of objectives. Generally, it only converges to the central part of the POF, leaving large uncovered areas at the borders. This behavior prevents obtaining solutions far away from the center of the POF, which may be interesting for decision-makers. The reason for this behavior is the existence of overlapping solutions, which severely affect the diversity of the population in the objective space. In order to solve this issue, this paper proposes a new algorithm, called NSGA-II/OSD, which hybridizes NSGA-II with an OSD strategy that increases the capabilities of NSGA-II to also reach areas further away from the center of the POF. Specifically, when there are two objectives, the objective space is divided into three regions that are defined by the nadir solution of the current population (see Figure 2), and the individuals in each of these regions are independently classified into Pareto fronts, with a region-specific Pareto optimal criterion. This particular feature allows solutions, that with the classical NSGA-II would have very poor ranking positions, to obtain a better ranking index with NSGA-II/OSD. The final effect of this "decomposed" Pareto ranking is that the diversity of the population is maintained, the presence of overlapping solutions is minimized and the exploration capability of the algorithm is improved. This decomposition strategy is illustrated in Figure 3.

Thus, the hybrid MOEA based on OSD for handling overlapping solutions proposed in this paper is a modified version of the bi-objective NSGA-II that includes an extra parameter $\alpha \left(0\le \alpha \le 1\right)$, which allows the convergence/diversity balance to be adjusted. The general flow chart of this proposed algorithm (see Figure 4) consists of the following steps:

1) Define the initial parameters: $N$, ${t}_{max}$, $\alpha$, Pcr, Pmut.

2) Randomly generate an initial population ${P}^{0}$ of size $N$.

3) Assign $t = 0$ and ${P}^{t} = {P}^{0}$.

4) Apply selection by tournament (best rank and crowding) and crossover and mutation operators to create a population of descendants ${Q}^{t}$ of size $N$.

5) Form a combined population ${R}^{t} = {P}^{t}\cup {Q}^{t}$.

6) Identify the nadir solution in the population ${R}^{t}$ by some method: ${\mathbf{z}}_{t}^{N} = {\left({z}_{1}^{N}, {z}_{2}^{N}\right)}^{T}$.

7) Classify ${R}^{t}$ on non-dominated fronts: ${\left\{{F}_{j}^{t}\right\}}_{j\ge 1}$.

8) If the stop condition is met ($t+1 = {t}_{max}$) go to step 15, else continue.

9) If $t > {\alpha t}_{max}$ go to step 13, else continue.

10) Based on the current nadir solution, the objective space is divided into three disjoint regions ${R}_{1}$, ${R}_{2}$ and ${R}_{3}$ as follows (see Figure 2):

11) Solutions from each region ${R}_{i}$ ($i = \mathrm{1, 2}, 3$) are classified independently on different non-dominated fronts: ${\left\{{F}_{j}^{it}\right\}}_{j\ge 1}$. For the ranking of the solutions, different optimization criteria are used for the objective functions ${f}_{1}$ and ${f}_{2}$:

(a) Region ${R}_{1}$: both objective functions, ${f}_{1}$ and ${f}_{2}$, are always maximized.

(b) Region ${R}_{2}$: the objective function ${f}_{1}$ is always maximized. The function ${f}_{2}$ is maximized if $t > \alpha {t}_{max}$ is verified, otherwise ${f}_{2}$ is minimized.

(c) Region ${R}_{3}$: the objective function ${f}_{2}$ is always maximized. The function ${f}_{1}$ is maximized if $t > \alpha {t}_{max}$ is verified, otherwise ${f}_{1}$ is minimized.

12) Combine the fronts obtained in the three regions: ${\left\{{F}_{j}^{t}\right\}}_{j\ge 1}$with ${F}_{j}^{t} = \bigcup _{i = 1}^{3}{F}_{j}^{it}$ (see Figure 3).

13) Add each front, in increasing order until a front ${F}_{k}^{t}$ cannot be fully integrated into the new population ${P}^{t+1}$ of size $N$, then classifying the ${F}_{k}^{t}$ front by crowding-distance to add solutions until the size of ${P}^{t+1}$ equals $N$.

14) Assign$t = t+1$ and${P}^{t} = {P}^{t+1}$, and go to step 4.

15) Add each front, in increasing order until a front ${F}_{k}^{{t}_{max}-1}$ cannot be fully integrated into the new population ${P}^{{t}_{max}}$ of size $N$, then classifying the ${F}_{k}^{{t}_{max}-1}$ front by crowding-distance to add solutions until the size of ${P}^{{t}_{max}}$ equals $N$.

16) End.

It should be noted that when the number $n$ of objectives increases, the number of objective solution space divisions is ${2}^{n}-1$. That is, in the case of $n = 3$, we would have 7 disjoint regions. However, in the case of 0-1 MOKP, the increase in the number of objectives leads to a significant decrease in the number of overlapping objective vectors in each population when moving from two to three or more objectives ^[12].

The computational complexity of NSGA-II/OSD is the same as that of NSGA-II. Note that in the case $\alpha = 0$, NSGA-II/OSD is actually NSGA-II.

4.   Experiments

The implementations of NSGA-II and MOEA/D have been carried out according to ^[5,6], respectively. All algorithms were applied on the 0-1 MOKP with two objectives and with 500 items (MOKP-2.500). The data of the problem and its respective POF can be downloaded from https://sop.tik.ee.ethz.ch/download/supplementary/testProblemSuite.

As mentioned in the Introduction section, there exist a large number of overlapping objective vectors in each population when EMO algorithms are applied to multi-objective combinatorial optimization problems with only a few objectives ^[12]. In particular, in the case of 0-1 MOKP, the increase in the number of objectives leads to a significant decrease in the number of overlapping objective vectors in each population when moving from two to three or more objectives. For this reason, the 0-1 MOKP with two knapsacks represents an ideal test-bed for analyzing the influence of overlapping solutions.

4.1.   Multi-objective knapsack problem

The multi-objective binary knapsack problem (0-1 MOKP) is a classic Operations Research problem that consists of a series of objects or items and a set of knapsacks. Each object has a weight and a profit associated with it, and the knapsacks has a maximum storage capacity associated with them. The issue is then to assign objects to each knapsack in such a way that the total profit yielded by the assigned objects is maximized, meanwhile the knapsack capacity cannot be exceeded by the total weight of the objects. The problem is of NP-hard complexity and can be used to model any real application that fits the model described in Eq (4.1).

Given a set of $m$ items and a set of $n$ knapsacks, the 0-1 MOKP can be stated as

where ${c}_{ij}\ge 0$ is the profit of item $j$ in knapsack $i$, ${w}_{ij}\ge 0$ is the weight of item $j$ in knapsack $i$, and ${c}_{i}$ is the capacity of knapsack $i$. ${x}_{i} = 1$ means that item $i$ is selected and put in all the knapsacks.

4.2.   Experiments: parameters and metrics configuration

All experiments used binary coding, uniform crossover with probability 0.8 and mutation probability (bit-by-bit) equal to $1/m$ with $m = 500$. The number of evaluations of the objective functions equal to 400, 000 was used as the stop condition. Three population sizes were used: $N = 50$, 100 and 200. The values examined for the α parameter, defined in Section 3 for the NSGA-II/OSD convergence/diversity balance, were ${\text{a}}_{i} = 0.1i$, and $i = 0, 1, 2, ..., 10$.

In the case of MOEA/D-based algorithms, the tested values of α or ${\alpha }_{0}$ parameter in the calculation of the reference point depend on the version considered. For Zhang^te defined in Subsection 2.3, the α values were fixed at $\alpha = 1.0, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3$, and for Ishibuchi^te defined in Subsection 2.4, the ${\text{a}}_{0}$ values were fixed at ${\text{a}}_{0} = 0.0, \ \ 0.1, \ \ 0.5, \ \ 1.0, \ \ 2.0, \ \ 3.0, \ \ 4.0, \ \ 5.0, \ \ 10.0$. It is worth mentioning that these two variations in the values of α and ${\text{a}}_{0}$ parameters for these versions of MOEA/D represent new results for the literature.

The hypervolume indicator suggested in ^[28] is used as a comparison measure between algorithms. The reference point considered for the calculation of the hypervolume was (0, 0), which guarantees that it is dominated by all the solutions generated during the evolution of the algorithms. In addition, the attainment surface concept introduced in ^[29] is used considering the approach suggested by Knowles ^[30]. You can visit: https://www.cs.bham.ac.uk/~jdk/plot_attainments, if you want to download the source code.

4.3.   NSGA-II/OSD pre-analysis

This subsection analyzes the effect on the convergence and diversity of the solution front reached when varying the α parameter of the proposed NSGA-II/OSD algorithm.

Figure 5 shows the box-and-whisker plots of the hypervolumes obtained by modifying the values of the α parameter, for the three population sizes considered, and the mean hypervolumes obtained in each case. It can be seen that the best hypervolume values are obtained for intermediate values of the α parameter (for further details see Table 1). Figure 6 shows the box-and-whisker plots of the number of solution repairs performed when modifying the values of the α parameter, for the three population sizes considered, and the mean number of solution repairs carried out in each case. It can be seen that the number of repairs carried out decreases as the parameter α increases. Figure 7 shows the 50% attainment surface for $N = 200$. For a more complete comparison, the true Pareto front is also shown. It can be seen that the algorithm gets closer to the true POF when $\alpha = 0$, which corresponds to NSGA-II, but obtains a less extensive approximate POF. As α increases, although the approximate POF becomes more and more extensive, it loses convergence.

According to the results shown in Table 1 and Figure 5, in Figure 7 it is confirmed that NSGA-II/OSD achieves the best 50% attainment surface for the intermediate values of α. In fact, for $\alpha > 0.8$ the approximate POF obtained is substantially different from the true POF. Not all 50% attainment surfaces have been drawn because, for values $0.2 < \alpha < 0.8$, the surfaces stick very close to the surface corresponding to $\alpha = 0.5$. However, all the results are shown in Table 1. Therefore, the choice of $\alpha = 0.5$ allows to obtain a good convergence/diversity balance.

Figure 8 shows the percentage of solutions in the regions ${R}_{2}$ and ${R}_{3}$ for $N = 200$. For $\alpha = 0.0$, which corresponds to the classic NSGA-II, it is observed that the percentage of solutions in these regions decreases rapidly to a value lower than 10%, which explains reaching an approximate POF less extensive than in the other cases (see Figure 7). However, for $\alpha > 0.0$ the rapid decrease in the percentage of solutions in the regions ${R}_{2}$ and ${R}_{3}$ occurs once it is verified that $t > \alpha {t}_{max}$. That is, for $\alpha > 0.0$ the algorithm initially explores the objective space for a percentage of generations, which allows it to obtain a more extensive approximation POF. This is where the correct choice of α allows to obtain a balance between the diversity (extension) and convergence of the approximate POF obtained. Note that, in the case of $\alpha = 1.0$, the percentage of solutions in the regions ${R}_{2}$ and ${R}_{3}$ always remains above 70%; which also helps to explain the greater diversity (and extent) of the approximate POF achieved. However, the trade-off is that this approximate POF is further away from the true POF (see Figure 7).

In addition, Figure 9 (top) shows the evolution of the mean percentage of overlapping solutions in the population for $N = 200$ and for the different values of α considered. It is observed that this percentage increases considerably as soon as it is verified that $> \alpha {t}_{max}$, where $t$ is the current generation and ${t}_{max}$ the maximum number of generations that are carried out, and then decreases again later on. Specifically, it is observed that, for $\alpha < 0.8$, towards the end of the evolution of the algorithm, the mean percentage of overlapping solutions tends to be around 35%. Figure 9 (bottom) clearly shows that the lowest percentage of overlapping solutions, throughout the evolution of the algorithm, is achieved for $\alpha = 1.0$, which explains the greater extension of the approximate POF reached (see Figure 7).

In view of the above and in order to obtain a good convergence/diversity balance, $\alpha = 0.5$ will be set to compare NSGA-II/OSD with other algorithms in Section 5.

4.4.   MOEA/D pre-analysis

The two reference point specifications described in Subsections 2.1 and 2.2 will be compared in this subsection. Figure 10 (left) shows the values of the mean hypervolume when varying the value of α in Zhang^te according to Eq (2.3). It is observed that, for $N = 100$ and $N = 200$, the value of the mean hypervolume clearly increases up to a maximum value (when $\alpha = 1.2$) and then decreases again for larger values of α (see Table 2). Figure 10 (right) shows the values of the mean hypervolume when varying the value of ${\alpha }_{0}$ in Ishibuchi^te according to Eqs (2.4) and (2.5). In this case, it is observed that the value of the mean hypervolume grows continuously as long as $\alpha < 3$, but for $\alpha \ge 3$ the hypervolume gets stagnated (see Table 3). Figure 11 shows that, in general, for both versions of MOEA/D, the number of repairs carried out decreases as the α parameter increases. This trend is clearer when Ishibuchi^te is used (sea Figure 11 (right)).

Figures 12 and 13 show the 50% attainment surface for $N = 200$ yielded by Zhang^te and by Ishibuchi^te, respectively. To get a better idea of the performance of the algorithms, the true POF is also depicted. On the one hand, according to the results shown in Table 2 and Figure 10 (left), in Figure 12 it is confirmed that Zhang^te achieves the best 50% attainment surface when $\alpha = 1.2$. On the other hand, according to the results shown in Table 3 and Figure 10 (right), in Figure 13 it is confirmed that Ishibuchi^te achieves the best 50% attainment surface when ${\alpha }_{0}\ge 3.0$.

In addition, Figure 14 shows the evolution of the mean percentage of overlapping solutions in the population for $N = 200$, when Zhang^te (left) and Ishibuchi^te (right) are used. It can be seen that, broadly speaking, as the α value increases, when Zhang^te is used, the mean percentage of overlapping solutions in the population increases. For example, for $\alpha = 1.2$ this percentage is slightly higher than 50%. Something similar happens with Ishibuchi^te algorithm, if the ${\alpha }_{0}$ value is increased (see Figure 14 (right)). Thus, for ${\alpha }_{0} = 3.0$ this percentage is slightly less than 60%.

Thus, for MOKP-2.500, it seems advisable to use the value $\alpha = 1.2$ in MOEA/D when using the Tchebycheff scalarization approach and Zhang's specification of the reference point, and to use the value ${\alpha }_{0} = 3.0$ in MOEA/D when using the Ishibuchi's specification of the reference point. Furthermore, in both cases, the mean number of repairs is similar (see Figure 11) and of the same order of magnitude as that obtained by NSGA-II/OSD ($\alpha = 0.4$) (see Figure 6 (bottom-right)). Therefore, in Section 5, these values of α and ${\alpha }_{0}$ will be used for comparison purposes.

5.   NSGA-II/OSD vs NSGA-II vs MOEA/D vs Global WASF-GA

This subsection compares the NSGA-II/OSD algorithm proposed in this work with the classical NSGA-II, two MOEA/D-based algorithms and the Global WASF-GA. For NSGA-II/OSD and the two versions of MOEA/D considered (MOEA/D Zhang^te and MOEA/D Ishibuchi^te) the values of α or ${\alpha }_{0}$ were set according to the experiments described in the previous sections, setting $\alpha = 0.5$ for NSGA-II/OSD, $\alpha = 1.2$ for MOEA/D Zhang^te and ${\alpha }_{0} = 3.0$ for MOEA/D Ishibuchi^te. All algorithms were applied on the 0-1 MOKP with two objectives and with 250, 500 and 750 items (MOKP-2.250, MOKP-2.500 and MOKP-2.750, respectively).

The data of the problems and its respective POF (of MOKP-2.250 and MOKP-2.500) can be downloaded from https://sop.tik.ee.ethz.ch/download/supplementary/testProblemSuite.

All comparisons used binary coding, uniform crossover of probability 0.8 and mutation probability (bit-by-bit) equal to $1/m$, where $m$ is the number of items. In addition, the population size ($N = 200$) and the number of generations ($G = 2, 000$) were fixed.

Tables 4-6, and Figures 15-26 show a detailed comparison between the five algorithms. Figures 15(left), 16(left) and 17(left) show the box-and-whisker plots based on the hypervolume approximation metric for MOKP-2.250, MOKP-2.500 and MOKP-2.750, respectively. From these figures, it is observed that NSGA-II/OSD obtains the best median and the lowest dispersion for MOKP-2.250, and it performs similarly to MOEA/D Zhang^te for MOKP-2.500 and MOKP-2.750. This is because the number of overlapping solutions of 0-1 MOKPs decreases as the number of items considered increases (see Table 5 and Figures 15(right), 16(right) and 17(right)). Moreover, for MOKP-2.500 and MOKP-2.750 can be seen that the average percentage of overlapping solutions in the population, of practically all generations, is lower with NSGA-II and NSGA-II/OSD. Finally, the best performance of NSGA-II/OSD and MOEAD algorithms over the NSGA-II and Global WASF-GA algorithms is also shown in Figures 18(left), 19(left) and 20(left), that present the average hypervolume evolution per generation for MOKP-2.250, MOKP-2.500 and MOKP-2.750, respectively.

From Table 5 it is extracted that at the end of the evolution, NSGA-II and NSGA-II/OSD have similar values of the average percentage of overlapping solutions, although with better hypervolume value for NSGA-II/OSD (see Table 4 and Figures 15(left), 16(left) and 17(left)). This is because in the first half of the generations, the average percentage of overlapping solutions is much lower for NSGA-II/OSD than for NSGA-II, which justifies that almost all the time NSGA-II/OSD has better average hypervolume value than NSGA-II.

Regarding the number of repairs required by the algorithms, Figures 21-23 show that, in general, the algorithm that has required the least number of repairs is NSGA-II/OSD, followed by both versions of MOEA/D. NSGA-II and Global WASF-GA have required significantly more repairs than the other algorithms.

Figures 24-26 show the 50% attainment surfaces. For a more complete comparison, the true Pareto front is also shown when is possible (MOKP-2.250 and MOKP-2.500). On the one hand, according to the results shown in Table 4 and Figure 15(left), in Figure 24 it is confirmed that NSGA-II/OSD achieves the best 50% attainment surface for MOKP-2.250. On the other hand, according to the results shown in Table 4 and Figures 16(right) and 17(right), in Figures 25 and 26 it is confirmed that NSGA-II/OSD achieves wider diversity of solutions in the 50% attainment surfaces, meanwhile MOEA/D-based algorithms reach better convergence.

Finally, a Wilcoxon rank-sum test ^[31] has been performed to make pairwise comparisons between the algorithms to study the significance of the results obtained. A significance level of 5% has been taken and the null hypothesis was "the two algorithms have the same average indicator", that is, the Wilcoxon test was used to determine if the difference between both means was statistically significant. If the p-value is less than 0.05, the null hypothesis is rejected and, then, we can conclude that the means of the two algorithms are comparable. In this case, a new Wilcoxon rank-sum test with one-sided alternative "greater" (or "less") and the same significance level is performed, that allows us to conclude which algorithm obtains better results than the other. The results are shown in Table 6 where a symbol ∆ indicates that the algorithm in the row has reached a better average hypervolume than the algorithm in the column, a symbol ∇ indicates that the algorithm in the row has reached a worse average hypervolume than the algorithm in the column, and a symbol - indicates that the comparative is not statistically significant. The results obtained indicate that, in general, NSGA-II/OSD and MOEA/D Zhang^te achieve better mean hypervolume than the other algorithms, MOEA/D Ishibuchi^te achieves better mean hypervolume than NSGA-II and Global WASF-GA, and NSGA-II achieves better mean hypervolume than Global WASF-GA. However, it should be noted that NSGA-II/OSD performs clearly better than MOEA/D Zhang^te when solving MOKP-2.250 (p-value = 2.537e-16), the difference between the two algorithms is not significant solving MOKP-2.500 (p-value = 0.4853) and NSGA-II/OSD performs better than MOEA/D Zhang^te when solving MOKP-2.750 (p-value = 0.007124). These results are in consonance with the above interpretations of Figures 15-17.

6.   Conclusions

A hybrid MOEA implemented from the NSGA-II and based on objective space division for handling overlapping solutions, called NSGA-II/OSD, is presented. A parameter $\alpha \in \left[0, 1\right]$ allows to balance between convergence and diversity of the front of non-dominated solutions achieved.

A preliminary study over the 0-1 MOKP-2.500 and using two quality indicators: the hypervolume metric measure and the attainment surface concept, allows us to establish that NSGA-II/OSD obtains the best mean hypervolumes for α values around 0.5. Clearly, for $0.2 < \alpha < 0.8$ the proposed algorithm obtains very good performance and better results when compared to NSGA-II. Note that in the case $\alpha = 0$, NSGA-II/OSD is actually NSGA-II.

A similar preliminary study was carried out to analyze performance of MOEA/D with the Tchebycheff scalarization approach, with both Zhang's specification of the reference point (Zhang^te) and Ishibuchi's specification of the reference point (Ishibuchi^te), over the 0-1 MOKP-2.500. In the first case, it was concluded that $\alpha = 1.2$ is a good choice. In the second case, ${\alpha }_{0}\ge 3.0$ performs well. In this case, for the comparisons, we set ${\alpha }_{0} = 3.0$ because the mean number of repairs of the solutions, in order to manage unfeasible solutions, is similar to Zhang^te with $\alpha = 1.2$.

Finally, the NSGA-II/OSD algorithm proposed in this work (with $\alpha = 0.5$) was compared with NSGA-II, MOEA/D Zhang^te ($\alpha = 1.2$), MOEA/D Ishibuchi^te (${\alpha }_{0} = 3.0$) and Global WASF-GA using the number of evaluations of the objective functions as stop criterion in all algorithms. In general, NSGA-II/OSD and MOEA/D Zhang^te achieve better mean hypervolume than the other algorithms, MOEA/D Ishibuchi^te achieves better mean hypervolume than NSGA-II and Global WASF-GA, and NSGA-II achieves better mean hypervolume than Global WASF-GA. However, it should be noted that NSGA-II/OSD performs clearly better than MOEA/D Zhang^te when solving MOKP-2.250, the difference between the two algorithms is not significant solving MOKP-2.500 and NSGA-II/OSD performs better than MOEA/D Zhang^te when solving MOKP-2.750.

Moreover, NSGA-II/OSD achieves the best 50% attainment surface for MOKP-2.250 and, for MOKP-2.500 and MOKP-2.750, it achieves a greater diversity of solutions in the 50% attainment surfaces, while the MOEA/D-based algorithms achieve better convergence to the true POF. This decrease in the performance of NSGA-II/OSD compared to MOEA/D-based algorithms, in terms of convergence, when the number of variables increases, can be attributed to the fact that, in general, when the number of variables increases, the number of overlapping solutions generated by MOEAs decreases, and their impact on the diversity of the approximate final front of non-dominated solutions found decreases as well.

NSGA-II/OSD has been proved to perform well with a well-known multi-objective combinatorial optimization problem that generates a large number of overlapping solutions when solved with MOEAs. As a future work, we pretend to study its behavior with other multi-objective combinatorial optimization problems that behave similarly in terms of overlapping solutions. Another interesting line of future work could be to compare the proposed algorithm with other diversity maintenance strategies.

Acknowledgments

The authors are grateful for partial support from the following sources: National Scientific and Technical Research Council CONICET Grant PIP Nº 11220150100777 and Universidad Nacional del Sur Grant PGI 24/J086.

Conflict of interests

All authors declare no conflicts of interest in this paper.