1 Introduction

Whenever facing a decision, there is often a set of objectives to optimize. For instance, when making a vacation plan with multiple destinations, one wants to minimize both the time spent in airports and the money spent on plane tickets. However, seldom can one obtain a solution that optimizes all objectives at once. It is usually the case that decreasing the value of an objective results in increasing the value of another. This occurs in many application domains [17, 22, 32].

In order to deal with multi-objective problems, we usually cast them into single-objective ones. For example, this can be achieved by defining a linear combination of the objective functions. Other option is to define a lexicographic order of the objectives [24], but this may result in unbalanced solutions where the first function is minimized while the remaining ones have a very high value.

In the multi-objective scenario, we are looking for Pareto-optimal solutions, i.e. all solutions for which decreasing the value of one objective function increases the value of another. After determining the set of all such solutions, known as Pareto front, one can select a representative subset and present it to the user [9].

Frameworks based on stochastic search have been developed to approximate the Pareto front of Multi-Objective Combinatorial Optimization (MOCO) problems [6, 33]. Several algorithms were also proposed based on iterative calls to a satisfiability checker, such as the Opportunistic Improvement Algorithm [8], among others [16]. Additionally, the Guided-Improvement Algorithm (GIA) [26] is implemented in the optimization engine of Satisfiability Modulo Theories (SMT) solver Z3 for finding Pareto optimal solutions of SMT formulas. New algorithms have also been proposed based on the enumeration of Minimal Correction Subsets (MCSs) [30] or P-minimal models [28]. A common thread to these algorithms is that they follow a SAT-UNSAT approach. A path diversification method has also been proposed where unsatisfiable cores are identified in order to cut the path generation procedure [31]. More recently, Maximum Satisfiability (MaxSAT) approaches have been used for MOCO [10, 12], but the proposed algorithms are limited to two objective functions.

In this paper, we propose two new algorithms for MOCO. The first algorithm is a core-guided approach that relies on encodings of the objective functions to effectively cut the search space in each SAT call. Additionally, we also propose a hitting set-based approach where the previous core-guided algorithm is used to enumerate a multi-objective hitting set. Note that these are the first algorithms for MOCO that take full advantage of unsatisfiable core identification over several objectives, as well as the first MOCO algorithm based on an hitting set approach, taking advantage of the duality between Pareto-MCSs [30] and unsatisfiable cores over several objectives. Experimental results show that the new algorithms proposed in this paper are complementary to the existing SAT-based algorithms for MOCO, thus extending the state-of-the-art tools for MOCO based on SAT technology.

The paper is organized as follows. Section 2 defines the MOCO problem and the standard notation used in the remainder of the paper. Next, Sections 3 and 4 describe the new core-guided and hitting set-based algorithms for MOCO. Experimental results and comparisons with other SAT-based algorithms are provided in Section 5. Finally, conclusions are presented in Section 6.

2 Preliminaries

We start with the definitions that fall in the SAT domain. Next, we introduce the definitions specific to solving the MOCO problem.

Definition 1

(Boolean Satisfiability problem (SAT)). Consider a set of Boolean variables \(V=\{x_1, \ldots , x_n\}.\) A literal is either a variable \(x_i \in V\) or its negation \(\lnot x_i \equiv \bar{x}_i\). A clause is a set of literals. A Conjunctive Normal Form (CNF) formula \(\phi \) is a set of clauses. A model \(\nu \) is a set of literals, such that if \(x_i \in \nu \), then \(\bar{x}_i \not \in \nu \) and vice versa.

The truth value of \(\phi \), denoted by \(\nu (\phi )\), is a function of \(\nu \), and is defined recursively by the following rules. First, the truth value of a literal is covered by \(\nu (x_i) = \top , \text { if }x_i \in \nu \), \(\nu (x_i) = \bot , \text { if }\bar{x}_i \in \nu \) and \(\nu (\lnot x_i) = \lnot \nu (x_i)\). Secondly, a clause c is true iff it contains at least one literal assigned to true. Finally, formula \(\phi \) is true iff it contains only true clauses,

$$\begin{aligned} \nu (\phi ) \equiv \bigwedge _{c \in \phi } \nu (c),\quad \nu (c) \equiv \bigvee _{l \in c} \nu (l). \end{aligned}$$
(1)

The model \(\nu \) satisfies the formula \(\phi \) iff \(\nu (\phi )\) is true. In that case, \(\nu \) is (\(\phi \)-)feasible.

Given a CNF formula \(\phi \), the SAT problem is to decide if there is any model \(\nu \) that satisfies it or prove that no such model exists.

Our algorithms require a SAT solver to be used as an Oracle. If the formula is satisfiable, then it returns a satisfiable assignment. Otherwise, the SAT solver returns with an explanation of unsatisfiability, called a core.

Definition 2

(Core \(\kappa \)). Given a CNF formula \(\phi \), we say a formula \(\kappa \) is an unsatisfiable core of \(\phi \) iff \(\kappa \subseteq \phi \) and \(\kappa \vDash \bot \).

Definition 3

(SAT solver). Let \(\phi \) be a CNF formula and \(\alpha \) a conjunction of unit clauses. We call \(\phi \) the main formula and \(\alpha \) the assumptions. A SAT solver solves the CNFFootnote 1 instance of the working formula \(\omega = \phi \cup \alpha \), i.e. decides on the satisfiability of \(\omega \).

A query to the solver is denoted by \( \phi \text {-}{} \texttt{SAT}(\alpha )\). The value returned is a pair \((\nu , \kappa )\), containing a feasible model \(\nu \) and a core of assumptions \(\kappa \), i.e. a subset of the assumptions \(\alpha \) contained in some core of \(\omega \). If the working formula \(\omega \) is not satisfiable, \(\nu \) does not exist, and the call returns \((\emptyset , \bullet )\). If \(\omega \) is satisfiable, the call returns \((\bullet , \emptyset )\).

Definition 4

(Relaxing/Tightening a formula). Given \(\phi \), a formula \(\psi \) is a relaxation of \(\phi \) iff \(\phi \vDash \psi \). We also say \(\psi \) relaxes \(\phi \). Conversely, \(\phi \) tightens \(\psi \).

Next we review Pseudo-Boolean formulas and optimization and define the MOCO problem.

Definition 5

(Pseudo-Boolean function, clause, formula (PB)). To any linear function \(\left\{ 0, 1\right\} ^n \rightarrow \mathbb {N}\), given by

$$\begin{aligned} g(\boldsymbol{x}) = g(x_1\ldots x_n) = \sum _i w_i x_i \quad w_i \in \mathbb {N},\quad x_i \in V, \end{aligned}$$
(2)

we call an (integer linear) PB function. Expressions like \(g(\boldsymbol{x}) \bowtie k, \quad \bowtie \;\in \lbrace \le , \ge , = \rbrace \), are called PB clauses. A PB formula is a set of PB clauses. For some model \(\nu : V \rightarrow \left\{ 0, 1\right\} \), let \(\boldsymbol{x}\) be the Boolean tuple \(\nu (V) \equiv (\nu (x_1),\ldots ,\nu (x_n))\). Given a formula \(\phi \), a model \(\nu \) is said (\(\phi \)-)feasible if it satisfies every clause in \(\phi \). The set of Boolean tuples \(Z(\phi ) = \left\{ \boldsymbol{x}= \nu (V)\in \left\{ 0, 1\right\} ^n: \nu (\phi )\right\} \) is called feasible space of the formula \(\phi \), and its elements \(\boldsymbol{x}\) are called feasible points. Any subset of the feasible space is called a \(\phi \)-feasible set.

Definition 6

(Pseudo-Boolean Optimization (PBO)). Let \(\phi \) be a PB formula, and f be a PB function. Then, minimize the value of the objective f over the feasible space \(Z(\phi )\) the formula \(\phi \). That is,

$$\begin{aligned} \text {find } \mathop {\mathrm {arg\,min}}\limits _{\boldsymbol{x}\in Z(\phi )} f. \end{aligned}$$
(3)

Multi-objective optimization generalizes PBO and builds upon a criterion of comparison (or order) of tuples of numbers. The most celebrated one is called Pareto order or dominance.

Definition 7

(Pareto partial order (\(\prec \))). Let \(Y\) be some subset of \(\mathbb {N}^n\). For any \(\boldsymbol{y}, \boldsymbol{y}'\in Y\),

$$\begin{array}{c} \boldsymbol{y}\preceq \boldsymbol{y}'\iff \forall i, \boldsymbol{y}_i \le \boldsymbol{y}'_i, \\ \boldsymbol{y}\prec \boldsymbol{y}'\iff \boldsymbol{y}\preceq \boldsymbol{y}'\wedge \boldsymbol{y}\ne \boldsymbol{y}'. \end{array}$$

We say \(\boldsymbol{y}\) dominates \(\boldsymbol{y}'\) iff \(\boldsymbol{y}\preceq \boldsymbol{y}'\). We say \(\boldsymbol{y}\) strictly-dominates \(\boldsymbol{y}'\) iff \(\boldsymbol{y}\prec \boldsymbol{y}'\).

Given a tuple of objective functions sharing a common domain \(X\), we can compare two elements \(\boldsymbol{x},\boldsymbol{x}'\in X\) by comparing the corresponding tuples in the objective space.

Definition 8

(Pareto Dominance (\(\prec \))). Let \(F: X\rightarrow Y\subseteq \mathbb {N}^n\) be a multi-objective function, mapping the decision space \(X\) into the objective space \(Y\). For any \(\boldsymbol{x}, \boldsymbol{x}'\in X\),

$$\begin{array}{c} \boldsymbol{x}\prec \boldsymbol{x}'\iff F(\boldsymbol{x}) \prec F(\boldsymbol{x}'), \\ \boldsymbol{x}\preceq \boldsymbol{x}'\iff F(\boldsymbol{x}) \preceq F(\boldsymbol{x}'). \end{array}$$

We say \(\boldsymbol{x}\) dominates \(\boldsymbol{x}'\) iff \(\boldsymbol{x}\preceq \boldsymbol{x}'\). We say \(\boldsymbol{x}\) strictly-dominates \(\boldsymbol{x}'\) iff \(\boldsymbol{x}\prec \boldsymbol{x}'\).

Contrary to the single-objective case, the consequence of this comparison criterion is that many different good solutions are mapped to different points in the objective space. Therefore, the solution to the problem is actually a set called Pareto front.

Definition 9

(Fronts). Given a a multi-objective function \(F:X\rightarrow Y\) and a feasible space \(Z\subseteq X\), the Pareto front of \(Z\) is a subset \(P\subseteq Z\) containing all elements that are not strictly-dominated,

$$\begin{aligned} P=\left\{ \boldsymbol{x}\in Z: \not \exists \boldsymbol{x}'\in Z: \boldsymbol{x}'\prec \boldsymbol{x}\right\} . \end{aligned}$$

We call img-front to the subset \(\overline{Y}\subseteq Y\) which is the image of \(P\) by \(F\),

$$\begin{aligned} \overline{Y}\equiv {{\,\mathrm{img\, front}\,}}_{Z} F= \left\{ \boldsymbol{y}\in Y: \exists \boldsymbol{x}\in P: \boldsymbol{y}= F(\boldsymbol{x})\right\} . \end{aligned}$$

Finally, we call arg-front of \(Z\), or simply \({{\,\textrm{front}\,}}\) of \(Z\), to any subset \(\overline{Z}\) of the Pareto Front P which is mapped by \(F\) into \(\overline{Y}\) in a one-to-one fashion

$$\begin{aligned} \overline{Z}= {{\,\textrm{front}\,}}_{Z} F. \end{aligned}$$

Definition 10

(Multi-Objective Combinatorial Optimization (MOCO)). Let \(F: X\rightarrow Y\subseteq \mathbb {N}^n\) be a multi-objective PB function, mapping the decision space \(X\subseteq \left\{ 0, 1\right\} ^n\) into the objective space \(Y\). Let \(Z\subseteq X\) be the feasible space of some PB formula \(\phi \), with variables in V. Then,

$$\begin{aligned} \text {find } {{\,\textrm{front}\,}}_{Z(\phi )} F. \end{aligned}$$
(4)

An instance will be denoted by the triple \(\langle \phi , V,F\rangle \).

Because the solutions of the problems are sets, bounds are now bound sets (Definition 13). In the single objective case, a bound is a value l such that \(\forall y = f (x): l \le y \), or equivalently, \(\not \exists y = f (x): l > y\). This equivalence is broken by the generalization. Each of the previous defining properties of a lower bound gives rise to a differently flavoured comparison of sets (Definitions 11 and 12).

Definition 11

(Set coverage). Let A and B be subsets of some decision space \(X\), equipped with a multi-objective function \(F\). Then, A covers B iff every element of B is dominated by some element of A, i.e. \(\forall b \in B, \exists a \in A: a \preceq b,\) and A strictly covers B iff \(\forall b \in B, \exists a \in A: a \prec b.\)

Definition 12

(Set non-inferiority). Let A and B be subsets of some decision space \(X\), equipped with a multi-objective function \(F\). Then A is non-inferior to B iff there is no element of B that strictly-dominates an element of A, \(\forall a \in A, b \in B: \lnot (a \succ b),\) and A is strictly non-inferior to B iff \(\forall a \in A, b \in B: \lnot (a \succeq b)\).

Note that in the single objective case, non-inferiority and coverage are the same. The next definition correctly generalizes the notion of lower bound.

Definition 13

(Bound sets). \(L \subseteq X\) is a (strictly) lower bound set of \(Z\subseteq X\) iff L (strictly) covers and is (strictly) non-inferior to \(Z\). If L is a lower bound set of \(Z\), we say \(L \preceq Z\). If it is a strictly lower bound set, we say \(L \prec Z\).

One way to generate a lower bound set of some Pareto front is to solve a related problem, where the formula is replaced by a relaxed version (Definition 4).

In our approach, we embed dominance relations into CNF formulas. We are interested in removing from the feasible space solutions that are dominated by some other known feasible solution. In order to do this, we make use of unary counters [3, 13, 14] that have been used to implement efficient PB satisfiability solvers.

Definition 14

(Unary Counter). Let \(f_{i}: \left\{ 0, 1\right\} ^n \rightarrow \mathbb {N}\) be a PB function and set V be an ordered set of variables that parametrize the domain of \(f_{i}\),

$$\begin{aligned} V = \left\{ x_1, \ldots , x_n\right\} , f_{i}(\boldsymbol{x}) = f_{i}(x_1, \ldots , x_n) \end{aligned}$$
(5)

Consider the CNF formula \(\widetilde{\phi }\) with variables \(V\cup O\), where \(O \cap V = \emptyset \) and O contains one variable \(o_{i,k}\) for each value \(k \in \mathbb {N}: \exists \boldsymbol{x}: k = f_{i}(\boldsymbol{x})\). The elements of \(O\) are the order variables. We call the tuple \(\langle f_{i}, V, O, \widetilde{\phi }\rangle \) an unary counter of \(f_{i}\) iff all feasible models \(\nu \) of \(\widetilde{\phi }\) satisfy

$$\begin{aligned} f_{i}(\boldsymbol{x}) \ge k \implies o_{i,k},\quad \boldsymbol{x}= \nu (V). \end{aligned}$$
(6)

3 Core-Guided Algorithm

Although core-guided algorithms for Maximum Satisfiability were initially proposed more than a decade ago [1, 2, 7, 21, 23], there is no such algorithm for MOCO. Hence, our goal is to take advantage of unsatisfiable cores identified by a SAT solver in order to lazily expand the allowed search space.

3.1 Algorithm Description

Fig. 1.
figure 1

Illustration of a run of Core-Guided (Algorithm 1) in the objective space. The img-front is the set \(\left\{ 1,2,3\right\} \). The fence bound \(\lambda \) gets updated at each iteration of the while loop at line 6, starting at \(\textrm{A}\) and ending at \(\varOmega \). The arrows are guided by the core \(\kappa \) (line 19). The green shading represents the evolution of the fence. Darker regions have been fenced for longer. The blue regions are blocked by optimal points. Darker regions are dominated by more points. We will be done in 7 iterations. After verifying that \(\textrm{A}\) is not feasible, we are instructed by the cores k to move along the diagonal twice. We find point 1 fenced. Therefore the associated \(\boldsymbol{x}\) is copied into \(I\) and the dominated region is blocked. We extend \(\lambda \) twice, and find point 2. After moving once more, we find part of the fence blocked, and the point branded with i is never generated. The next movement stations \(\lambda \) at \(\varOmega \). Point 3 is found. The Oracle acknowledges we are done, by returning \(\kappa =\emptyset \) (line 15): she knows that no movement of \(\lambda \) will extend \(I\).

figure a

Algorithm 1 presents the pseudo-code for an exact core-guided algorithm for MOCO. Figure 1 illustrates an abstract execution of the algorithm.

Let \(\langle \phi , V, F\rangle \) be a MOCO instance. Recall that \(\phi \) denotes the set of PB constraints, V is the set of variables and \(F\) denotes the list of m objective functions. First, the algorithm starts by building a working formula with the problem constraints and an unary counter for each objective function (lines 3-4). This is accomplished by the call to . Next, a vector \(\lambda \) of size m is initialized with the lower bound of each objective function (line 5), assumed to be 0 for simplicity.

At each iteration of the main loop, the assumptions \(\alpha \) are assembled from order variables \(o\), chosen with the value of \(\lambda \) in mind (line 7). The call to Footnote 2 returns the next smallest value belonging to the image of the objective i. Given the semantics of the order variables \(o_{i,k}\) (Definition 14), the tuple \(\lambda \) fences the search space, i.e. \(\nu \) satisfies \(\alpha \) only if the corresponding tuple \(\boldsymbol{x}\) satisfies \(F(\boldsymbol{x}) \preceq \lambda \). If the SAT call (line 10) returns a solution (i.e. \(\nu \ne \emptyset \)), \(\boldsymbol{x}\) is stored in and all dominated solutions are removed from I (line 11). Moreover, one can readily block all feasible solutions dominated by \(\boldsymbol{x}\) using a single clause (line 13) [28].

Usually, there are several feasible fenced solutions. This occurs because the algorithm may increase multiple entries of \(\lambda \) at once. In any case, the inner while loop (lines 9-14) collects all such solutions.

When the working formula \(\widetilde{\phi }\) becomes unsatisfiable, the SAT solver provides a core \(\kappa \). If \(\kappa \) is empty (line 15), then the unsatisfiability does not depend on the assumptions, i.e. it does not depend on temporary bounds imposed on the objective functions. At that point, we can conclude that no more solutions exist that are both satisfiable and not dominated by an element of \(I\). As a result, the algorithm can safely terminate (line 16). Otherwise, the literals in \(\kappa \) denote a subset of the fence walls \(\lambda _i\) that may be too restrictive, in the sense that unless we increment them (line 19) no new non-dominated solutions can be found.

3.2 Algorithm Properties

Lemma 1

The img-front \(\overline{Y}\) of \(I\cup Z(\widetilde{\phi })\)(Definition 9) is not changed by the inner loop (lines 9-14).

Proof

Consider some particular iteration of the internal loop. Line 11 and line 13 remove all elements of \(I\cup Z(\widetilde{\phi })\) that are dominated by the feasible point \(\boldsymbol{x}\). Line 11 filters the explicit set \(I\), line 13 filters the implicit set \(Z(\widetilde{\phi })\). Solutions that are strictly dominated by \(\boldsymbol{x}\) cannot be mapped into an element of \(\overline{Y}\). The other solutions \(\boldsymbol{x}'\) that are filtered out must attain the same objective vector attained by \(\boldsymbol{x}\), \(F(\boldsymbol{x}') = F(\boldsymbol{x})\). Because \(\boldsymbol{x}\) is also inserted at line 11, removing \(\boldsymbol{x}'\) will not disturb \(\overline{Y}\).

Lemma 2

At the start of each iteration of the external loop (lines 6-19), every solution in \(I\) is optimal, and no two elements of \(I\) attain the same objective vector.

Proof

We prove this by contradiction. Assume that there is a non-optimal solution \(\boldsymbol{x}\in I\) at the start of the external loop (line 6). In the first iteration, this does not occur because \(I\) is empty. Hence, this can only occur if the inner loop (lines 9-14) finishes with a non-optimal solution \(\boldsymbol{x}\in I\).

The inner loop (lines 9-14) enumerates solutions inside the fence defined by \(\lambda \). We know that \(F(\boldsymbol{x}) \preceq \lambda \) because it is inside the fence and the entries of \(\lambda \) never decrease. If \(\boldsymbol{x}\) is non-optimal, then there must be an optimal solution \(\boldsymbol{x}'\) such that \(F(\boldsymbol{x}') \prec F(\boldsymbol{x}) (\preceq \lambda )\). Hence, \(\boldsymbol{x}'\) is also inside the fence. As a result, \(\boldsymbol{x}'\) must be found before the inner loop finishes, since at each iteration only dominated solutions are blocked (line 13). If \(\boldsymbol{x}\) is found before \(\boldsymbol{x}'\), then \(\boldsymbol{x}\) is excluded from \(I\) (line 11) when \(\boldsymbol{x}'\) is found. Otherwise, if \(\boldsymbol{x}'\) is found first, then \(\boldsymbol{x}\) is not found by the SAT solver (blocked at line 13) because it is dominated by \(\boldsymbol{x}'\). Therefore, we cannot have a non-optimal solution \(\boldsymbol{x}\in I\) at the end of the inner loop or at the start of each iteration of the external loop (lines 6-19). Furthermore, no two elements of \(I\) attain the same objective vector since when a solution \(\boldsymbol{x}\) is found, all other solutions \(\boldsymbol{x}'\) such that \(F(\boldsymbol{x})=F(\boldsymbol{x}')\) are also blocked (line 13).

Lemma 2 establishes a weaker form of anytime optimality. The elements of the incumbent list \(I\) are not necessarily optimal at anytime, but they are optimal immediately after completing the inner loop. It is easy enough to make it anytime optimal. This could be achieved if the algorithm refrains from adding solutions directly to \(I\) in the inner loop and maintain a secondary list, where it stores the solutions that are still not necessarily optimal. This list takes the role of \(I\) inside the inner loop. After completing the inner loop, all elements of the secondary list are optimal, and can be safely transferred to the main list I.

Proposition 1

Algorithm 1 is sound.

Proof

If the algorithm returns, \(Z(\widetilde{\phi }\wedge \alpha ) = \emptyset \). Because \(\kappa \) is empty, no core of the unsatisfiable formula \(\widetilde{\phi }\wedge \alpha \) intersects \(\alpha \), and \(\widetilde{\phi }\) is also unsatisfiable, \(Z(\widetilde{\phi }) = \emptyset \). Using Lemma 1 both at the end and at the start of the course of the algorithm, the img-front of I is the img-front of \(Z(\widetilde{\phi })\), with \(\widetilde{\phi }\) given by line 4. Because the order variables are only restricted by the unary counter formula, the img-front of \(Z(\widetilde{\phi })\) is the img-front of \(Z(\phi )\). Therefore \(I\) must contain an arg-front of the problem. Using Lemma 2, every element of \(I\) is optimal, and there is no pair \(\boldsymbol{x},\boldsymbol{x}'\in I\) such that \(F(\boldsymbol{x})=F(\boldsymbol{x}')\). Therefore, \(I\) is an arg-front of the MOCO instance.

Proposition 2

Algorithm 1 is complete.

Proof

The inner loop will always come to fruition, because in the worst case it will generate every feasible solution dominated by the current \(\lambda \) once, and the feasible space is finite.

If the algorithm does not return for some particular instance, then \(\kappa \) is never empty. In that case, every iteration of the external loop starting at line 6 will increase at least one of the entries of \(\lambda \). Eventually, one entry i must achieve the upper limit of \(f_i\), and the order variable retrieved by \(o_{i,\lambda _i + 1}\) will not exist. Because the evolution of \(\lambda _i\) is monotonous, the assumptions will contain at most \(m - 1\) variables, from that point on. By the same token, the assumptions \(\alpha \) will eventually be empty, and so must be \(\kappa \ \subseteq \alpha \), contradicting the assumption that the algorithm never terminates.

4 Hitting Set-based Algorithm

This section proposes a MOCO solver based on the enumeration of hitting sets. The main idea is to compute a sequence of relaxations \(\psi \) of the formula \(\phi \), and solve the corresponding problems. The front \(\overline{T}\) of the relaxed problem gets incrementally closer to the desired front \(\overline{Z}\), and will eventually reach it.

4.1 Algorithm Description

figure d
Fig. 2.
figure 2

Illustration of a run of the Hitting-Sets (Algorithm 2) in the objective space. The Pareto front is the set \(\left\{ 1,2,3\right\} \), and the feasible solutions are marked by . For each iteration of the main while loop at line 2 we get a narrower lower bound \(\overline{T}\) (line 4), culminating in the solution. We are done in 3 iterations, marked by A, B and . The shading represents the number of iterations whose freshly found points dominate the region. The lighter tone was painted by A, the darker one by all three. We start with the empty formula (line 1) and get A. Because the only point in A is not feasible, we tighten the relaxation (line 13). Iteration B generates one feasible point, 1, which is therefore optimal. Note that the region dominated by 1 can be pruned from now on. The other point is used to tighten the formula once more. Lastly, the lower bound contains the feasible points 2 and 3 in addition to 1, which was already found, and the algorithm stops.

Algorithm 2 contains the pseudo-code for our hitting set-based algorithm for MOCO. Figure 2 illustrates an abstract execution of the algorithm.

The algorithm starts by setting the relaxed formula \(\psi \) to empty (line 1). The main loop that starts at line 2 hones the relaxation until we get the desired result. At each iteration, we solve the current relaxed formula \(\psi \) at line 4. This is accomplished by using some MOCO solver. Because this amounts to computing a lower bound set, the Core-Guided algorithm, previously described, is a good choice for the task. We anticipate that it performs well for problems whose front is in the vicinity of the origin, given that by construction, the focus of its search is biased to that region. Notice that the first relaxation’s arg-front is the set that contains the origin only (assuming all literals in the objective functions are positive). We expect that the first few relaxations will stay close to it.

Next, for each element \(\boldsymbol{x}\) in \(\overline{T}\) (the Pareto-front of \(\psi \)), we check the \(\phi \)-feasibility of \(\nu : \nu (V) = \boldsymbol{x}\), using the assumptions mechanism, and return a (possibly empty) core of assumptions \(\kappa \). The assumptions \(\alpha _{\boldsymbol{x}}\) built at line 6 are a set of unit clauses whose polarity is inherited from \(\nu \),

$$\begin{aligned} \nu (x_i) \implies x_i \in \alpha , \quad \lnot \nu (x_i) \implies \lnot x_i \in \alpha . \end{aligned}$$
(7)

Assuming \(\phi \) is satisfiable, the returned core \(\kappa \) will be void iff \(\alpha _{\boldsymbol{x}} \wedge \phi \) is satisfiable. In this case, \(\boldsymbol{x}\) corresponds to an optimal solution.

The diagnosis \(\varDelta \) is central for the algorithm. Intuitively, it reports if and why the relaxed problem’s solution is different from the true Pareto solution. We add every non-empty \(\kappa \) to the diagnosis \(\varDelta \) (line 9). In the end, \(\varDelta \) is empty iff every element of the relaxed front \(\overline{T}\) is \(\phi \)-feasible. At that point, we have found a \(\phi \)-feasible lower bound set. All such sets are arg-fronts, and so the algorithm terminates (line 11). Otherwise, if \(\varDelta \) is not empty, then the found cores are added to the relaxed formula \(\psi \) (line 13). This step ensures all tentative points produced in line 4 hit all previously found unsatisfiable cores, and that the algorithm advances in a monotonous fashion towards the solution.

4.2 Algorithm Properties

Given a MOCO instance \(\langle \phi ,V,F\rangle \), the formula \(\phi \) encodes the feasible space \(Z\) implicitly, which in turn defines the desired front \(\overline{Z}\). This is a many to one correspondence, in the sense that there are many different values of \(\psi \) that encode the same Pareto front. It may happen that some of the counterpart instances are easier to solve than the original one, which begs the question: given \(\phi \), can we effectively find a simpler formula \(\psi \) with the same Pareto front? This is the motto of the proposed algorithm. It is done by iteratively honing a relaxed formula (Definition 4).

The main idea is to compute a sequence of relaxations that get incrementally tighter. In that case, the corresponding front \(\overline{T}\) gets incrementally closer to the desired front \(\overline{Z}\),

$$\begin{aligned} \phi{} & {} \implies{} & {} \psi _n{} & {} \implies{} & {} \ldots{} & {} \implies{} & {} \psi _{1},\end{aligned}$$
(8)
$$\begin{aligned} \overline{Z}{} & {} \succeq{} & {} \overline{T}_n{} & {} \succeq{} & {} \ldots{} & {} \succeq{} & {} \overline{T}_1, \end{aligned}$$
(9)

where \(\overline{Z}\) is one of the desired arg-fronts, and \(\overline{T}_i\) is an arg-front of \(\psi _{i}\).

Lemma 3

Consider some multi-objective function \(F: X\rightarrow Y\). Let \(Z, T\) be subsets of \(X\), such that \(T\subseteq Z\). Then, any arg-front of \(T\) is a lower bound set of any arg-front of \(Z\) (Definition 13), i.e. \(T\subseteq Z\implies \overline{T} \preceq \overline{Z}\).

Lemma 3 is true because optimizing over a superset of some feasible space always returns a (non-strict) lower bound set. In a sense, the optimization can only be more extreme when applied to the superset. In particular, the feasible space of a relaxed formula is a superset of the original one. This is why the chain of \(\preceq \) relations in Equation (9) is correct.

Lemma 4

Let \(\phi \) be a formula, \(Z\subseteq X\) be its feasible space and \(F: X\rightarrow Y\) be some multi-objective function. Let L be a lower bound set of the Pareto front of \(Z\). Then, any element \(x \in L\) that is feasible belongs to the Pareto front, \(L \cap Z\subseteq P\). If all elements \(\boldsymbol{x}\in L\) are feasible, then L is an arg-front.

Lemma 4 implies that every lower bound set with only feasible elements must be itself an arg-front (this is an exact analogy with the single-objective case, where lower bound set is replaced by infimum and arg-front by arg-min.) By construction of the diagnosis \(\varDelta \), this is equivalent to the condition used in Algorithm 2 to decide if it can terminate.

To ensure the sequence gets to \(\overline{Z}\) in a finite number of steps, we need more than a string of relaxations. Each entry \(\psi '\) must be strictly tighter than the predecessor \(\psi \).

Lemma 5

Consider Algorithm 2. Let \(\psi \) be the relaxed formula at some iteration, and \(\psi '\) be the relaxed formula at the next iteration. Then,

  1. 1.

    \(\psi \) relaxes \(\psi '\), i.e. \(\psi '\vDash \psi \);

  2. 2.

    Both \(\psi \) and \(\psi '\) relax \(\phi \), i.e. \(\phi \vDash \psi , \phi \vDash \psi '\);

  3. 3.

    \(\psi '\) does not relax \(\psi \), i.e. \(\psi \not \vDash \psi '\);

Proof

Each statement will be proven in turn.

The first is true because \(\psi \subseteq \psi '\), by construction (line 13).

We prove the second by induction on the number of iterations. Initially, \(\psi \) is empty. Therefore, \(\psi \) relaxes any formula, in particular \(\phi \). Assume \(\phi \vDash \psi \) for some iteration. Consider one of the clauses \(\lnot \kappa \) added at line 13. We know that \(\phi \wedge \kappa \) is unsatisfiable. Therefore, \( \phi \wedge \kappa \vDash \bot \implies \lnot (\phi \wedge \kappa ) \vDash \top \iff \phi \vDash \lnot \kappa \). Given the assumption \(\phi \vDash \psi \), we get \(\phi \vDash \psi \wedge \lnot \kappa \). Repeating the process for the other added clauses \(\lnot \kappa _i\), we get \(\phi \vDash \psi \wedge \lnot \kappa _1 \ldots \wedge \lnot \kappa _n \equiv \psi '\).

Assume \(\psi '\) is a relaxation of \(\psi \). Then, any \(\psi \)-feasible model \(\nu \) is also \(\psi '\)-feasible. We will prove there is at least one model that violates this. To start, note that it only makes sense to consider \(\psi '\) if there is some non-empty core \(\kappa \) in the diagnosis \(\varDelta \); otherwise, the algorithm would have terminated before updating \(\psi \) into \(\psi '\). Let \(\kappa \) be one element of \(\varDelta \), generated at line 7 while \(\psi \) is current. Consider the Boolean tuple \(\boldsymbol{x}\in \overline{T}\) used to build the assumptions of the query that generated \(\kappa \). Let \(\nu : \nu (V) = \boldsymbol{x}\). The model \(\nu \) is \(\psi \)-feasible, because it is part of the arg-front of \(\psi \). The model \(\nu \) satisfies \(\kappa \) because \(\kappa \subseteq \alpha _{\boldsymbol{x}}\) and the way \(\alpha _{\boldsymbol{x}}\) is constructed (line 6, Equation (7)). Therefore, \(\nu \) does not satisfy \(\lnot \kappa \). Because \(\lnot \kappa \subseteq \psi '\), \(\nu \) cannot satisfy \(\psi '\), i.e. there is at least one \(\psi \)-feasible model that is not \(\psi '\)-feasible.

Proposition 3

Algorithm 2 is sound.

Proof

By Lemma 5, \(\psi \) relaxes \(\phi \) and therefore \(\overline{T}\) solves a relaxation of the original problem. By Lemma 3, it is a lower bound set of \(\overline{Z}\). When the algorithm returns, all elements of \(\overline{T}\) are feasible. By Lemma 4, \(\overline{T}\) must be an arg-front.

Proposition 4

Algorithm 2 is complete.

Proof

Assume Algorithm 2 never ends, implying \(\overline{T}\) is never completely feasible (i.e. \(\overline{T}\nsubseteq Z\)). The number of relaxed feasible spaces \(T\) is finite. If Algorithm 2 does not end, it will enumerate all of them, never repeating any: at any iteration, the updated relaxed formula effectively blocks the reappearance of any feasible space seen before, because by Lemma 5 the updated value \(\psi '\) strictly tightens \(\psi \). Then, this sequence is necessarily finite, and so must be the number of iterations. But in that case, Algorithm 2 must end, and we have a contradiction.

Consider the sequence whose entries are the value of \(F(\overline{T})\) computed at the beginning of each iteration of the main loop at line 2. The last element of this sequence is the solution. It may happen that for some i, the entries indexed by i and \(i+1\) are the same. Therefore, the sequence may include blocks of contiguous entries that share the same value. In the worst case scenario, there are many different arg-fronts for the same img-front, and the algorithm ends up enumerating all of them without any movement in the objective space. We expect the algorithm will be effective whenever a few of the relaxed problems are enough to get to the full solution. Otherwise, we can end up solving an exponential number of problems.

5 Experimental Results

This section evaluates the performance of the algorithms proposed Footnote 3 in Sections 3 and 4. These algorithms are compared against other SAT-based MOCO solvers.

5.1 Algorithms and Implementation

The Core-Guided algorithm proposed in Algorithm 1 uses the selection delimiter encoding [14] that has been shown to be more compact. Next, the selection delimiter encoding is extended to produce a unary encoding for each objective function. Additionally, an order encoding [29] is also used. We refer the interested reader to the literature for further details on this and other encodings [13,14,15, 27]. Observe that any unary encoding from PB into CNF can be used.

The Hitting-Sets algorithm implements Algorithm 2. This hitting set-based approach uses Algorithm 1 to find the relaxed arg-front (line 4 of Algorithm 2).

The P-Minimal algorithm implements a SAT-UNSAT approach based on the enumeration of P-Minimal models [28]. This algorithm is implemented with the same PB to CNF encoding as the Core-Guided. Finally, the ParetoMCS is based on the stratified enumeration of Minimal Correction Subsets. We used the publicly available implementation of ParetoMCSFootnote 4.

5.2 Experimental Setup and Benchmark Sets

The following MOCO problems are considered: the multi-objective Development Assurance Level (DAL) Problem [5], the multi-objective Flying Tourist Problem (FTP) [22], the multi-objective Set Covering (SC) Problem [4, 28] and the multi-objective Package Upgradeability (PU) Problem [11]. All instances are publicly available from previous research work or were generated from real-world data.

The DAL benchmark set (95 instances) encodes different levels of rigor in the development of a software or hardware component of an aircraft. The development assurance level defines the assurance activities aimed at eliminating design and coding errors that could affect the safety of an aircraft. The goal is to allocate the smallest DAL to functions to decrease the development costs [18].

The FTP benchmark set (129 instances) encodes the problem of a tourist that is searching for a flight travel route to visit n cities. The tourist defines her home city, the start and end of the route. She specifies the number of days \(d_i\) to be spent on each city \(c_i\) (\(1 \le i \le n\)) and also a time window for the complete trip. The problem is to find the route that minimizes the time spent on flights and the sum of the prices of the ticketsFootnote 5.

The SC benchmark set (60 instances) is a generalization of the set covering problem and was used in previous research work [28]. Let X be some ground set and A a cover of X. Each element in A has an associated cost tuple. The goal is to find a cover of X contained in A that Pareto-optimizes the overall cost.

The PU benchmark set (687 instances) were generated from the Package Upgradeability benchmarks [19] from the Mancoosi International Solver Competition [20]. The packup tool [25] was used to generate these benchmarks that contain between two and five objectives to optimize.

All results were obtained on an Intel Xeon Silver 4110 CPU @ 2.10GHz, with 64 GB of RAM. Each tool was executed on each instance with a time limit of 1 hour and 10 GB of RAM memory limit.

5.3 Results and Analysis

Table 1 shows the number of instances whose Pareto front is completely enumerated, for each algorithm and benchmark set. Overall, the new unsatisfiability-based algorithms proposed in the paper completely solve more instances than the ParetoMCS and the P-Minimal algorithms. Note that the ParetoMCS is the one that solves fewer instances since it needs to enumerate all MCSs. The Core-Guided and Hitting-Sets converge faster to the Pareto front due to their UNSAT-SAT approach, while the P-Minimal is slower to converge. Overall, the Core-Guided algorithm is able to solve more instances than the other algorithms.

Table 1. Number of MOCO instances whose complete solution is found and certified per algorithm and benchmark set. Best results are in bold.

All tested algorithms are exact, but in some cases only an approximation of the Pareto front could be found within the time limit. However, the partial solution that is returned may still be valuable. In order to evaluate the quality of the approximations provided by each tool, we use the Hypervolume (HV) [34] indicator. HV is a metric that measures the volume of the objective space dominated by a set of points in the objective space, up to a given reference point. The coordinates of the reference point chosen are the maximal values of each objective. Regions that are not dominated by a reference front are discarded (we combined the results for each algorithm in order to produce the reference front). Larger values are preferred. A normalization procedure is carried out so that the values of HV are always between 0 and 1.

Fig. 3.
figure 3

Comparison of the HV results for each set of instances. Each series is sorted independently, smaller values first. Vertical scale is logarithmical.

Figure 3 shows a cactus plot of the HV for all tools on each benchmark set. The P-Minimal provides better quality approximations of the Pareto front in the DAL (Figure 3a) and PU (Figure 3d) benchmarks since it uses a SAT-UNSAT approach. Hence, it is faster to find an approximation to the Pareto front. Moreover, since some of the instances in these sets have higher optimal values on the objective functions, the Core-Guided and Hitting-Sets take many interactions until they reach the feasible part of the search space. Despite performing an unsatisfiability-based search, Core-Guided and Hitting-Sets algorithms are still able to provide good quality solutions since when these algorithms find solutions, these are in the Pareto front. Moreover, observe that even in these sets of instances, Core-Guided is still able to find all the Pareto front in more instances.

The ParetoMCS is able to provide good quality approximations in the FTP (Figure 3b) and PU (Figure 3d) benchmarks. Note that ParetoMCS does not use an explicit representation of the objective functions. The FTP instances have several large coefficients in the objective functions, but the representation used in Core-Guided is still effective for these instances. Observe that the performance of both algorithms is similar in the FTP dataset.

The Hitting-Sets finds poor approximations for all datasets. A common feature of this algorithm is the need to enumerate many hitting sets before being able to find feasible solutions. Hence, in several instances it is unable to provide good approximations. However, it is still able to prove optimality for more instances in the SC benchmark set than the P-Minimal algorithm.

Overall, the Core-Guided is the best performing algorithm being able to find the complete Pareto frontier in more instances. This is due to the fact that in many cases, it does not need to relax all variables to find solutions in the Pareto front. Moreover, when evaluating the quality of the approximations, it is still able to outperform the other approaches on the FTP and SC benchmark sets, despite applying an unsatisfiability-based approach.

6 Conclusions

This paper proposes two new algorithms for Multi-Objective Combinatorial Optimization (MOCO). The first is a core-guided approach, while the second is based on the enumeration of hitting sets. These are the first SAT-based algorithms that fully integrate these strategies into a MOCO solver.

Experimental results on different sets of benchmark instances show that the new core-guided approach results in a robust algorithm that outperforms other SAT-based algorithms for MOCO. Using unary counters to express Pareto dominance in CNF proved to be an effective way to harness the power of SAT solvers in solving MOCO. The ability to express concepts related to dominance makes the algorithms conceptually simple.

Overall, the new algorithms are able to completely enumerate the Pareto front for more instances than previous SAT-based approaches. Moreover, despite following an unsatisfiability-based approach, the newly proposed algorithms are also able to provide good quality approximations even when they are unable to completely enumerate the Pareto front. Hence, these new unsatisfiability-based algorithms extend the state of the art for MOCO solvers by complementing and improving upon the existing tools based on queries to SAT Oracles.