Stratification for Constraint-Based Multi-Objective Combinatorial Optimization

1 Introduction

In Multi-Objective Combinatorial Optimization (MOCO) applications such as scheduling [Iturriaga et al., 2017] or green computing [Zheng et al., 2016], there is more than one objective to be optimized, but there is no pre-defined hierarchy among the objective functions. Hence, there may exist multiple optimal solutions, known as Pareto optimal solutions, each of them favoring certain objectives at the expense of others. Several approaches have been proposed [Jackson et al., 2009] that try to identify the Pareto front, i.e., all Pareto optimal solutions. However, this is known to be very hard for large MOCO instances, and most algorithms are just able to provide an approximation to the Pareto front.

There is a wide plethora of stochastic algorithms that try to approximate the Pareto front [Deb et al., 2000; Xu and Fortes, 2010], and these methods are known for exhibiting a very good performance when the problem instances are not tightly constrained. However, their performance usually deteriorates when the constraint set becomes harder to satisfy.

On the other hand, constraint-based methods are able to deal with large constraint sets, but have a hard time in integrating information from the several objectives to be optimized.

Lately, it has been shown that one can find the Pareto front of a MOCO instance by enumerating the set of Minimal Correction Subsets (MCSs) [Terra-Neves et al., 2017] or the Pminimal models [Soh et al., 2017] of a propositional formula. The main strength of these methods is on being able to quickly satisfy a large constraint set, but in several cases produce lower quality solutions than stochastic algorithms.

Constraint-based solvers for single objective problems are known to use different techniques on how to deal with the objective function. One example is to partition the terms in the objective function according to its weights. To this end, stratification techniques have been proposed [Ans´otegui et al., 2012]. In this work, we improve constraint-based MOCO algorithms by focusing on the optimization of several objective functions. In particular, we first propose a new solver that uses stratification in solving MOCO instances. Our main contributions are as follows: (1) formalization of the usage of stratification in MCS algorithms; (2) incorporation of stratification in solving MOCO instances; (3) a new stratification heuristic and (4) an extensive experimental evaluation on green computing MOCO instances that shows a huge improvement in performance enabled by stratification.

In this paper we start by defining MOCO and MCSs in section 2. Section 3 explains how to integrate stratification in MCS algorithms. Section 4 proposes the use of stratification in a constraint-based MOCO solver. Moreover, different stratification strategies are discussed and a new heuristic is also described. Experimental results showing the effectiveness of the newly proposed techniques are presented in section 5. Finally, the paper concludes in section 6.

2 Definitions

This section starts by describing Pseudo-Boolean Optimization (PBO) and Multi-Objective Combinatorial Optimization (MOCO). Next, Minimal Correction Subsets (MCSs) are defined. Finally, we review how MCSs can be used to find solutions for PBO and MOCO problems.

2.1 Multi-Objective Combinatorial Optimization

Let X = {x1, . . . , xn} be a set of n Boolean variables. A literal is either a variable xi or its complement ¬xi. Given a set of m literals l1, . . . , lm and their respective coefficients ω1, . . . , ωm ∈ N, a Pseudo-Boolean (PB) expression is a weighted sum of literals ∑ ωi · li. Given an integer k ∈ N, a linear PB constraint has the form:

Consider a set F = {c1, . . . , cm} of m PB constraints and a cost function f defined over a set of X Boolean variables. The PBO [Boros and Hammer, 2002] problem consists of finding a complete assignment α : X → {0, 1} that satisfies all constraints in F, denoted α(F) = 1, and minimizes the value of f, denoted f(α). If F is unsatisfiable, then α(F) = 0 for any assignment α. Analogously, given a PB constraint c, α(c) = 1 (α(c) = 0) denotes that α satisfies (does not satisfy) c.

Example 2.1. Consider a PBO instance defined over X = {x1, x2, x3} where F = {(x1 + x2 + x3 ≥ 2)} is the set of PB constraints and f(X) = 4x1 + 2x2 + 3x3 is the cost function. In this case, α = {(x1, 0), (x2, 1), (x3, 1)} is an optimal assignment with a cost of 5. All other assignments to X either do not satisfy F , or result in a higher value of f(X).

In this paper, the negation operator ¬ can also be applied to a given PB constraint c. In this case, ¬c represents that c cannot be satisfied. Observe that ¬c can be easily converted into an equivalent PB constraint by updating the relational operator and the right-hand side. Finally, note that, by definition, all coefficients are required to be positive, but constraints with negative coefficients can easily be converted to equivalent ones with all positive coefficients. The same can be applied for the cost function.

A MOCO [Ulungu and Teghem, 1994] instance is composed of two sets: a set F = {c1, . . . , cm} of constraints that must be satisfied and a set O = {f1, . . . , fk} of cost functions to minimize. In this work, we focus on the special case where c1, . . . , cm are PB constraints and f1, . . . , fk are PB expressions over a set X of Boolean variables.

Given two complete assignments α, α′ such that α ̸= α′ and α(F) = α′ (F) = 1, we say that α dominates α′ , written α ≺ α′ , if, and only if, ∀f∈Of(α) ≤ f(α′ ) and ∃f ′∈Of ′ (α) < f ′ (α′ ). α is said to be Pareto optimal if, and only if, no other complete assignment α′′ exists such that α′′ (F) = 1 and α′′ ≺ α. In MOCO, the goal is to find the set of Pareto optimal solutions, also referred to as Pareto front.

Example 2.2. Let F = {(x1 + x2 + x3 ≥ 2)} be the set of constraints and O = {(2x1 + x2), (2¬x2 + 2x3)} the set of cost functions of a MOCO instance. The costs for each possible satisfiable assignment to F are shown in Table 1. The lines that correspond to Pareto optimal solutions are highlighted in bold, while the other are dominated. For example, {(x1, 1), (x2, 0), (x3, 1)} is not Pareto optimal because it is dominated by {(x1, 0), (x2, 1), (x3, 1)}.

For most problems, finding all Pareto optimal solutions in a reasonable amount of time is extremely hard. Hence, our focus is on finding the best approximation of the Pareto front.

2.2 Minimal Correction Subsets

Let F be an unsatisfiable set of PB constraints. A Minimal Correction Subset (MCS) of F is a subset C ⊆ F such that F \ C is satisfiable and C is minimal, i.e., (F \ C) ∪ {c} is unsatisfiable for all c ∈ C.

Example 2.3. Consider the unsatisfiable set of PB constraints F = {(x1 +x2 = 1), (x1 ≥ 1), (x2 ≥ 1)}. F has three MCSs C1 = {(x1 ≥ 1)}, C2 = {(x2 ≥ 1)} and C3 = {(x1 + x2 = 1)}.

There are several algorithms described in the literature for finding MCSs [Bailey and Stuckey, 2005; Felfernig et al., 2012; Marques-Silva et al., 2013; Menc´ıa et al., 2015]. Optionally, these algorithms can receive as argument a set of hard constraints FH that must be satisfied and a set FS of soft constraints for which we want to find an MCS C ⊆ FS. In this case, C is an MCS if FH ∪ (FS \ C) is satisfiable and FH ∪ (FS \ C) ∪ {c} is unsatisfiable for all c ∈ C. For simplicity, we assume that FH is always satisfiable, but this can be checked using a single call to a satisfiability solver.

MCSs can be used to find approximate solutions of PBO instances. Let F be the set of constraints and f(X) = ∑ ωi · li the cost function of a PBO instance. Let L(f) be the set of all literals in f and L¬(f) the set of clauses built from the negation of the literals in L(f), i.e., L¬(f) = ⋃li∈L(f) {(¬li)}. Applying an MCS algorithm with FH = F and FS = L¬(f), produces an MCS C of L¬(f). We abuse notation and denote as f(C) the cost of C, defined as:

Any assignment that satisfies F ∪ L¬(f) \ C will have a cost of f(C), which provides an approximation to the optimum of the PBO instance. Actually, the PBO problem can be reduced to finding the MCS C ⊆ L¬(f) that minimizes f(C) [Birnbaum and Lozinskii, 2003].

MCSs can also be used to find the Pareto front of MOCO instances. Let F and O be the constraint and cost function sets, respectively, of a MOCO and L¬(O) = ⋃f∈O L¬(f). Terra-Neves et al. [2017] proved that one can find the Pareto front by enumerating all MCSs of L¬(O).

Example 2.4. Consider F = {(x1 + x2 + x3 ≥ 2)} and O = {(2x1 + x2), (2¬x2 + 2x3)} from example 2.2. In this case, we have L¬(O) = {(¬x1), (¬x2), (x2), (¬x3)}. As a result, there are three MCSs: C1 = {(¬x2), (¬x1)}, C2 = {(¬x2), (¬x3)} and C3 = {(¬x1), (x2), (¬x3)} with costs (3, 0), (1, 2) and (2, 4), respectively. Observe that these MCSs include all the Pareto optimal assignments highlighted in table 1.

3 Stratification

Stratification has been successfully applied to boost the performance of Maximum Satisfiability (MaxSAT)

Algorithm 1: Generic stratified MCS algorithm

Input: FH, FS 1 (P1, P2, . . . , Pk) ←Partition(FS) 2 (S1,C) ← (FH, ∅) 3 for i ← 1 to k do 4 Ci ←MCS(Si, Pi) 5 Si+1 ←Si ∪ (Pi \ Ci) ∪ ⋃c∈Ci {¬c} 6 C ←C ∪Ci 7 return C

solvers [Ans´otegui et al., 2012]. Its use has been suggested by Menc´ıa et al. [2015] for approximating MaxSAT with MCSs, but stratified MCS enumeration has not been properly described in the literature yet. In this section, we describe a generic stratified approach that can be used to find lower cost MCSs faster with any MCS algorithm. Given a cost function f(X) = ∑ ωi · li, the main goal of stratification is to focus the search on satisfying the literals in L¬(f) with larger coefficients. The procedure starts by partitioning L¬(f) into k sets P1, P2, . . . , Pk such that all literals in Pi have larger coefficients than those in Pi+1 for all 1 ≤ i < k. Next, an MCS algorithm is applied considering all hard constraints and just the literals in P1 (the ones with larger coefficients) as soft. In the next iteration, P2 is added to the set of soft constraints. At each iteration, a new partition is added and the procedure ends when all partitions have been considered.

The pseudo-code for the generic MCS stratified approach is presented in algorithm 1. It receives as input a set of hard constraints FH and a set of soft constraints FS and it returns an MCS C of FS. First, it starts by partitioning FS into a set of k partitions such that FS = ⋃ki=1 Pi. At each iteration, Si denotes the set of constraints that are to be satisfied. In the first iteration, all hard constraints are considered (line 2). On the other hand, C denotes the set of unsatisfied soft constraints. Then, for each partition Pi, it computes an MCS Ci of Pi with Si as the set of hard constraints (line 4). Since the constraints in Pi \ Ci are satisfied, these are added to Si+1 (line 5). Moreover, the negation of the constraints in Ci is also added to Si+1 to further restrict the search space in the next MCS call. This is safe since Ci is an MCS of Pi, so satisfying Pi \ Ci entails that the constraints in Ci cannot be satisfied. Finally, the constraints in Ci are added to C (line 6). Next, we prove that, after all k iterations, C is an MCS of FS.

Proof. Suppose that C is not an MCS of FS. The two possible scenarios are: (1) FH ∪ (FS \ C) is unsatisfiable, or (2) there exists c ∈ C such that FH ∪ (FS \ C) ∪ {c} is satisfiable. Consider the first scenario 1 . Ck is an MCS of Pk, considering Sk as hard constraints, so Sk ∪ (Pk \ Ck) is satisfiable. Note that, for any 1 ≤ i, j ≤ k, Ci ⊆ Pi and if c ∈ Ci and c ∈ Pj, then c ∈ Cj [Terra-Neves et al., 2017].

1 Note that we can ignore the negation of the constraints added in line 5, since these are entailed by the remainder of the formula.

Algorithm 2: Stratified MCS-based MOCO algorithm

Input: F , O = {f1, f2, . . . , fl} 1 for i ← 1 to l do 2 (P i 1 , P i 2 , . . . , P iki ) ← Partition(L¬(fi)) 3 (P1 . . . Pp) ←Combine((P 1 1 . . . P 1 k1 ) . . . (P l1 . . . P lkl )) 4 (S1,C) ← (F, ∅) 5 for i ← 1 to p do 6 Ci ←MCS(Si, Pi) 7 Si+1 ←Si ∪ (Pi \ Ci) ∪ ⋃c∈Ci {¬c} 8 C ←C ∪Ci 9 return C

Therefore, we have Sk ∪ (Pk \Ck) = FH ∪ ⋃kj=1 (Pi \Ci) = FH ∪ (FS \ C), contradicting scenario (1). Now let us consider the second scenario. Let Ci be an MCS of Pi such that c ∈ Ci. Then, Si ∪ (Pi \ Ci) ∪ {c} is unsatisfiable. We have Si ∪ (Pi \ Ci) ∪ {c} = FH ∪ ⋃ij=1 (Pj \ Cj) ∪ {c} that is a subset of FH ∪(FS \C)∪{c}, contradicting scenario (2).

Algorithm 1 computes a single MCS C, but it can be used to find another MCS. This can be done by adding a blocking constraint to FH such that at least one constraint in C must be satisfied and re-executing the algorithm. Hence, algorithm 1 can be used to enumerate all MCSs of ⋃ki=1 Pi by blocking previous MCSs in subsequent invocations of the algorithm.

4 Stratification for Multi-Objective Combinatorial Optimization

Stratification for MOCO is not straightforward for two reasons: (1) unlike the single-objective case, a good MOCO algorithm should not only converge fast but also produce an approximation with high diversity in regard to costs; (2) cost functions can be very different in nature. This section describes how stratification can be extended to MOCO in order to find better approximations of the Pareto front and faster.

Algorithm 2 presents the pseudo-code for the usage of stratification for solving MOCO instances. The algorithm has as input a constraint set F and a cost function set O = {f1, . . . , fl}. The algorithm starts by using stratification for each cost function fi (line 2). Hence, the set of soft clauses L¬(fi) is split into ki partitions P i 1 P i 2 . . . P iki . Next, a single partition sequence is built by mixing the partitions of each cost function (line 3), where we have p = ∑li=1 ki. For example, assuming we have two cost functions f1 and f2 with partitions P 1 1 P 1 2 and P 2 1 P 2 2 P 2 3 , a possible sequence could be P 1 1 P 2 1 P 2 2 P 1 2 P 2 3 . Finally, the sequence of partitions is solved as in algorithm 1, generating an MCS of ∪fi∈OL¬(fi).

Recall that the Pareto front can be identified using MCS enumeration. Therefore, using algorithm 2 to enumerate all MCSs enables the generation of the Pareto front of a MOCO instance. However, for most pratical instances, generating the complete Pareto front is not feasible in a reasonable amount of time. Hence, our goal is to use this approach to find the best possible approximation of the Pareto front by generating a diverse set of good quality solutions in a given time limit.

In our algorithm, the partitioning of cost functions is performed separately, since the cost functions can be very different in nature. Observe that if we repeatedly use the same sequence for MCS enumeration, some cost functions can be favored at the expense of others. For example, if our partition sequence is P 1 1 P 2 1 P 2 2 P 1 2 P 2 3 , algorithm 2 will focus on satisfying the high cost literals of L¬(f1) first, which might entail higher costs for f2. Therefore, f1 will be favored at the expense of f2, which violates the diversity requirement. In order to promote diversity, the combination procedure (line 3) should generate a different partition sequence in each call to algorithm 2. In our implementation, partitions are mixed by selecting the next partition of a cost function with uniform probability at each iteration. For instance, assuming we are building a sequence and we already have P 1 1 P 2 1 , the next partition would be P 1 2 or P 2 2 , each having a probability of 1 2 of being selected.

4.1 Partitioning Heuristics

The partitioning heuristic is an important component of the stratified approach and can have a significant impact on performance. In general, given a cost function f, we sort the literals in L¬(f) in decreasing order of their respective coefficients. Then, for each literal l ∈ L¬(f), we add l to the current partition P and check whether P is a good partition according to some criteria. If so, P is saved and the next partition is generated using the remaining literals.

Partitioning heuristics differ in the criteria for deciding if P is a good partition. One possibility is to simply split each cost function into a fixed number p of partitions. In the FIXED heuristic, P is considered a good partition if p·|P | ≥ |L¬(f )|. This heuristic is used in Microsoft's SMT solver Z3 [Bjørner et al., 2015] for solving MaxSAT problems. In this paper we propose the Literal-Weight Ratio (LWR) heuristic that is a variant of the heuristic proposed by Ans´otegui et al. [2012]. In LWR, P is a good partition if the ratio between the number of literals in P and the number of distinct coefficients among those literals is above a given threshold β, i.e., |P| |W| > β, where W is the set of coefficients of the literals in P.

In both heuristics, our implementation does not produce partition sets where literals with the same weight are assigned to different partitions, i.e. we handle all literals with the same weight in the same step when generating new partitions.

4.2 Combining Stratification and Division Reduction

It was recently proved that, in the presence of a cost function f of the form f(X) = ∑ki=1 gi(X) hi(X) , f can be replaced by k cost function pairs gi and -hi, as long as gi(X) ≥ 0 and hi(X) > 0 for all 1 ≤ i ≤ k, without sacrificing Paretooptimal solutions of the original problem [Terra-Neves et al., 2018]. This technique is known as division reduction and the set of gi and -hi functions is referred as reduction product.

Suppose that we have a MOCO with cost functions f1 and f2, where f1 is a sum of k divisions and f2 is a regular PB expression. The reduced MOCO will have 2k+1 cost functions, where 2k of them are the reduction product of f1. The stratification process, as described in the previous section, would assign a uniform probability of selection to each function, resulting in an overall probability of 2k 2k+1 of selecting a partition of f1, which is much larger than f2's probability of 1 2k+1 . Therefore, stratification is very likely to favor f1, violating the diversity requirement. To prevent this behavior, the partitioning process must be aware that some cost functions may be the product of division reduction.

We implemented two strategies for combining stratification and division reduction. The first strategy, referred to as MERGE, merges the 2k literal sets L¬(gi),L¬(-hi) of the reduction product of f into a single set L¬(f), i.e., L¬(f) = ⋃ki=1 (L¬(gi) ∪ L¬(-hi)) Then, L¬(f) is partitioned as described in section 4.1, with the exception that the coefficients in the reduction product are considered instead of the ones in f. The second strategy relies on adapting the selection probabilities of reduced cost functions instead of considering them as one. Recall the previous example with f1 and f2, where f1 is a sum of k divisions. Initially, both have a selection probability of 1 2 . After the division reduction of f1, we split f1's probability uniformly among the 2k functions in its reduction product, resulting in a probability of 1 2 · 1 2k = 1 4k for each function, while f2 maintains a probability of 12 . In what follows, we refer to this strategy as SPLIT.

5 Experimental Results

In this section, the performance of the stratified MCS approach for MOCO is evaluated on instances of the Virtual Machine Consolidation (VMC) problem. In VMC, we have several servers with fixed resource capacities and Virtual Machines (VMs) with requirements of those same resources. Each VM must be placed in some server, but server capacities cannot be exceeded and some VMs cannot be placed in the same server. There exists an initial placement, i.e., a VM can be associated with an initial server, incurring a migration cost if the VM is placed in a different one. A migration budget constraint can be used to enforce an upper limit on the migration costs, and is specified as a percentile bp of the total memory capacity of the servers. The goal is to find a placement for all VMs that satisfies the constraints and simultaneously minimizes (1) energy consumption in the data center, (2) migration costs and (3) resource wastage. The latter is a measure of the imbalance of server resource usage.

A detailed description of the VMC problem and the MOCO formulation, can be found in the literature [Terra-Neves et al., 2017; Zheng et al., 2016]. Note that the resource wastage cost function is a sum of divisions, which is handled using division reduction. The evaluation is performed on publicly available VMC benchmarks 2 , based on subsets of workload traces randomly selected from the Google Cluster Data project 3 .

In order to evaluate the quality of the approximations to the Pareto front, two quality indicators are used. The Hypervolume (HV) quality indicator [Zitzler and Thiele, 1999] provides a combined measure of convergence and diversity. Larger values of HV mean that the solution set is composed of solutions of better quality and/or diversity. The Inverted Generational Distance (IGD) indicator [Zhang and Li, 2007] measures the average Euclidean distance, in the cost space, between the Pareto optimal solutions and the approximation returned by the algorithm, and a smaller value is preferred. Since the Pareto front is unknown, we use instead the combination of the solutions produced by all algorithms evaluated.

2 http://sat.inesc-id.pt/dome/

3 http://code.google.com/p/googleclusterdata/

Figure 3: Distribution of number of MCSs for division handling (bp = 100%).

In our approach, any MCS algorithm can be used. However, in order to make a proper evaluation and comparison with the work of Terra-Neves et al. [2017], we also use the ClauseD (CLD) algorithm [Marques-Silva et al., 2013] to compute MCSs. All algorithms are implemented in Java and Sat4j-PB v.2.3.4 [Le Berre and Parrain, 2010] is used as satisfiability checker. Moreover, we merge partitions whenever the satisfiability checker reaches a fixed number of conflicts (200000 in our experiments). This results from the fact that some calls might be too constrained. When this occurs, we relax the formula by adding the next partition. For example, if the satisfiability algorithm reaches the conflict limit when computing an MCS of Pi, the MCS algorithm stops, Pi+1 is added and the search continues for an MCS of Pi ∪ Pi+1. If all partitions have been added, there is no conflict limit.

Each algorithm was ran with a memory limit of 4 GB and a time limit of 1800 seconds. Randomized algorithms were executed with 10 different seeds for each instance, and the analysis is performed using the median values over all executions. Finally, it was observed that algorithms with better performance have lower standard deviation values than algorithms with worse performance. This occurred for both hypervolume and IGD indicators.

5.1 Partitioning Heuristics

First, we evaluate the impact of the partitioning heuristics (section 4.1) on the performance of algorithm 2. Figures 1 and 2 show the distributions of the HV and IGD values, respectively, obtained by the different heuristic configurations with a migration budget of 100%. A point (x, y) in the HV (IGD) distribution plot indicates that the given configuration obtained an HV (IGD) equal to or greater (lower) than y for x instances. For example, the point (50, 0.56) on line 'LWR (β = 15)' in figure 1 indicates that, using the LWR heuristic with β = 15, the algorithm obtained HVs equal to or greater than 0.56 for 50 instances.

Different configuration values were tried for the LWR and FIXED heuristics. In particular, LWR was configured with values of β ∈ {2, 5, 10, 15, 20, 30} and FIXED with p ∈ {10, 20, 40, 80, 160}. In order to improve readability, results are shown solely for representative values. However, all best configurations are included. Additionally, in figure 2, IGD values larger than 1.6 are not displayed. This is done also for readability of the IGD plots.

Observe that LWR is strictly better than FIXED both in terms of HV and IGD, except when β = 2 in the IGD plot. However, with β = 2 was the worst configuration of LWR. For all other values of β, the LWR performance is very similar. The same also occurs for different configurations of the FIXED heuristic. Hence, it was observed that both heuristics are robust, and that the partitioning strategy itself has a much higher impact in performance than parameter fine tuning.

5.2 Handling Division Reduction

This section evaluates the MERGE and SPLIT strategies proposed in section 4.2 for handling the division reduction of the wastage optimization criteria. Both were run using the LWR partitioning heuristic with different values of β.

MERGE and SPLIT showed very similar IGD and HV distributions, i.e. they produced approximations of very similar quality of the Pareto front. However, the number of MCSs found was very different. Therefore, we choose to show the number of MCSs found for representative configurations in figure 3. Note that the number of MCSs is a relevant metric since it can be seen as a measure of progress. The more MCSs an algorithm is able to find, the closer it is to proving that the Pareto front was found instead of an approximation. Figure 3 clearly shows that MERGE is able to find many more MCSs than SPLIT. At best, SPLIT with β = 15 is able to find close to 230 MCSs for one instance, while MERGE is able to find at least as much for close to 20 instances. Actually, MERGE is able to find close to 500 MCSs for some instances.

5.3 Comparison with State-of-the-Art

In this section, we compare SCLD, our new stratified approach for MOCO, with the simple CLD enumeration algorithm and the MCSEnumPD variant that incorporates an MCS diversification technique [Terra-Neves et al., 2018]. We also compare with the VMPMBBO [Zheng et al., 2016] algorithm for VMC and with the general purpose evolutionary algorithms NSGAII [Deb et al., 2000] and MOEAD [Zhang and Li, 2007]. SCLD was configured with the LWR partitioning heuristic (β = 15) and the MERGE division reduction strategy. MCSEnumPD, VMPMBBO and NSGAII were configured as suggested in the literature [Zheng et al., 2016; Terra-Neves et al., 2017; 2018]. The only difference is in how we encoded the individuals in NSGAII's population. In this work, the regular integer encoding is used instead of the binary integer encoding, since it produces better results on the VMC problem. Further details on encodings for evolutionary algorithms can be found in the literature [Rothlauf, 2006]. MOEAD was configured with crossover and mutation rates of 0.8 and 0.05 respectively, a population size of 100, a neighborhood size of 20, a 0.9 probability of crossover with individuals from the neighborhood and a maximum of 2 individuals that can be replaced by a single offspring.

Figure 7: IGD distributions (bp = 5%).

Figure 9: IGD distributions (bp = 1%).

Figures 4 and 5 show the HV and IGD distributions of each algorithm for VMC instances with bp = 100%. Observe that the new techniques proposed in this paper allow SCLD to obtain a huge performance improvement when compared with previous state-of-the-art constraint-based algorithms for MOCO, CLD and MCSEnumPD. Moreover, considering the hypervolume indicator, SCLD is competitive with stochastic algorithms NSGAII, MOEAD and VMPMBBO for these instances. Nevertheless, stochastic algorithms are still able to obtain better results measured by IGD. These results are indicative that NSGAII is able to converge faster, but SCLD finds approximations with higher diversity.

Figures 6, 7, 8 and 9 present the results when migration budgets of VMC instances are constrained to 5% and 1%. Observe that NSGAII's, MOEAD's and VMPMBBO's performance degrades considerably as the budget decreases because these algorithms have a much harder time dealing with more tightly constrained instances. On the other hand, constraint-based methods for MOCO thrive in such scenarios. Nevertheless, note that SCLD is still able to improve considerably on both CLD and MCSEnumPD, thus confirming that our approach still improves previous algorithms even on tighter instances. Overall, SCLD is the first constraintbased algorithm to be competitive with stochastic approaches for the VMC instances where bp = 100%, and is the best performing algorithm when bp ≤ 5%.

6 Conclusion and FutureWork

Interest has been surging in solving MOCO formulations using constraint-based methods. Until recently, the state of the art was mostly populated by stochastic methods, which are able to find good approximations fast, but are unsuitable for constrained problems. On the other hand, constraint-based methods are able to quickly satisfy tightly constrained problems, but have a harder time producing good quality approximations. This paper proposes the use of stratification as a means to guide MCS-based MOCO algorithms towards good quality solutions more effectively. Experimental results in a large set of MOCO instances show that stratification provides a very significant performance improvement for MCS algorithms, being now able to compete with (and often outperform) stochastic methods in terms of approximation quality.

As future work, we propose to improve the diversity of the approximation by dynamically adapting the cost function selection probabilities, and to exploit these in order to account for user preferences [Deb and Sundar, 2006]. Moreover, recent work in high-school timetabling problems has shown that significant improvements can be obtained by combining stochastic and constraint-based methods [Demirovic and Musliu, 2017]. Such a combined algorithm for MOCO could reap the benefits of both worlds with none of the weaknesses.

Acknowledgments

This work was supported by national funds through FCT with references UID/CEC/50021/2013 and SFRH/BD/111471/2015.