# Exact and approximate determination of the Pareto front using Minimal Correction Subsets

A. P. Guerreiro, J. Cortes, D. Vanderpooten, C. Bazgan, I. Lynce, V.

Manquinho, J.R. Figueira

► To cite this version:

A. P. Guerreiro, J. Cortes, D. Vanderpooten, C. Bazgan, I. Lynce, et al. Exact and approximate determination of the Pareto front using Minimal Correction Subsets. Computers and Operations Research, 2023, 153, pp.106153. 10.1016/j.cor.2023.106153. hal-04792039

HAL Id: hal-04792039

https://hal.science/hal-04792039v1

Submitted on 23 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Exact and approximate determination of the Pareto front using Minimal Correction Subsets

A. P. Guerreiro<sup>a,*</sup>, J. Cortes<sup>a</sup>, D. Vanderpooten<sup>b</sup>, C. Bazgan<sup>b</sup>, I. Lynce<sup>a</sup>, V. Manquinho<sup>a</sup>, J.R. Figueira<sup>c</sup>

## ARTICLE INFO

### Keywords:

* Multi-objective boolean optimization
* Minimal Correction Subsets
* Approximation algorithms

## ABSTRACT

Recently, it has been shown that the enumeration of Minimal Correction Subsets (MCS) of Boolean formulas allows solving Multi-Objective Boolean Optimization (MOBO) formulations. However, a major drawback of this approach is that most MCSs do not correspond to Pareto-optimal solutions. In fact, one can only know that a given MCS corresponds to a Pareto-optimal solution when all MCSs are enumerated. Moreover, if it is not possible to enumerate all MCSs, then there is no guarantee of the quality of the approximation of the Pareto frontier. This paper extends the state of the art for solving MOBO using MCSs. First, we show that it is possible to use MCS enumeration to solve MOBO problems such that each MCS necessarily corresponds to a Pareto-optimal solution. Additionally, we also propose two new algorithms that can find a $(1+\epsilon)$-approximation of the Pareto frontier using MCS enumeration. Experimental results in several benchmark sets show that the newly proposed algorithms allow finding better approximations of the Pareto frontier than state-of-the-art algorithms, and with guaranteed approximation ratios.

## 1. Introduction

In the last decade, new algorithms and tools for Boolean optimization have been developed that rely on successive calls to a highly effective Propositional Satisfiability (SAT) solver. As a result, these new algorithms have been used in several application domains such as timetabling (Asín and Nieuwenhuis, 2014), fault localization in C programs (Jose and Majumdar, 2011), and design debugging (Safarpour et al., 2007), among others (Feng et al., 2017).

Despite its success in Boolean optimization with a single objective function, there are few algorithms for Multi-Objective Boolean Optimization (MOBO) that are based on successive calls to a SAT solver. However, there are some exceptions. For instance, the Guided Improvement Algorithm (GIA) (Jackson et al., 2009), is implemented in the optimization engine of the Z3 solver (Bjørner et al., 2015) for finding Pareto-optimal solutions of Satisfiability Modulo Theories (SMT) instances with multiple objective functions. More recently, new logic-based algorithms have been proposed that are based on the enumeration of Minimal Correction Subsets (MCSs) (Terra-Neves et al., 2017) or P-minimal models (Soh et al., 2017). These logic-based algorithms have shown to be successful at solving MOBO instances when the constraint set is hard to satisfy.

Logic-based algorithms suffer from some drawbacks. Although MCSs are known to provide good approximations in specific contexts, if not all MCSs are enumerated, there is no formal guarantee on the quality of the provided solutions (Marques-Silva et al., 2013). Similarly, while a P-minimal model corresponds to a Pareto-optimal solution (Soh et al., 2017), if not all P-minimal models are found, there is no guarantee that the Pareto-optimal solutions found are representative of the Pareto front (e.g. the found solutions can be skewed to one objective function). Moreover, this approach also requires an expensive encoding of the objective functions into SAT. In some cases, the SAT formula becomes so large that the SAT solver is unable to handle the formula (Eén and Sörensson, 2006; Roussel and Manquinho, 2009).

In this work, logic-based algorithms for MOBO are extended in several ways. First, we start by proving that if one uses the order encoding to encode the objective functions into SAT, then there is a one-to-one correspondence between an MCS and a Pareto-optimal solution. Moreover, we propose new logic-based algorithms to obtain an *a priori* guaranteed $(1+\epsilon)$-approximation of the Pareto frontier. These algorithms are based on approximate encodings of the objective functions that are usually much smaller than the full encoding

<footnote><sup>a</sup> INESC-ID, Instituto Superior Técnico, Universidade de Lisboa, Portugal</footnote>

<footnote><sup>b</sup> LAMSADE (CNRS UMR7243), University Paris-Dauphine, PSL Research University, France</footnote>

<footnote><sup>c</sup> CEGIST, Instituto Superior Técnico, Universidade de Lisboa, Portugal</footnote>

<footnote>* Corresponding author.
*E-mail address:* andreia.guerreiro@tecnico.ulisboa.pt (A. P. Guerreiro).</footnote>
into SAT. Finally, an extensive experimental analysis is performed on the proposed algorithms which includes a comparison with other state-of-the-art logic-based algorithms for MOBO.

The paper is organized as follows. The background information is introduced in Section 2. New ways of encoding the objective functions and their properties are explained in Section 3, these are at the core of the approximation algorithms introduced in Section 4. The algorithms are then experimentally compared against state-of-the art ones in Section 5. Finally, concluding remarks are drawn in Section 6.

## 2. Background: Concepts, definitions, and notation

This section provides the relevant background information for the paper. The formalization of multi-objective Boolean optimization problems is introduced, as well as the main concepts, definitions and corresponding notation used in the remainder of the paper.

### 2.1. Multi-objective Boolean optimization

The Multi-Objective Boolean Optimization (MOBO) problem can be defined as minimizing *p* objective functions defined over a set of *n* (Boolean) variables and *m* constraints as follows:

**Problem 1 (Multi-Objective Boolean Optimization).**

$$
\begin{array}{ll}
\min & f(x) = (f_1(x) = \sum_{j=1}^{n} c_j^1 l_j, \dots, f_k(x) \\
& = \sum_{j=1}^{n} c_j^k l_j, \dots, f_p(x) = \sum_{j=1}^{n} c_j^p l_j) \\
\text{subject to:} & \sum_{j=1}^{n} a_{ij} l_j \ge b_i, & i \in \{1 \dots m\} \\
& l_j \in \{x_j, \bar{x}_j\}, x_j \in \{0, 1\}, & j \in \{1 \dots n\}
\end{array}
$$

In this formulation, constraints are linear Pseudo-Boolean constraints (Roussel and Manquinho, 2009), and each objective function $f_k(x)$ is also a linear expression defined over a set of Boolean variables. Each coefficient, $a_{ij}$, and each right-hand-side, $b_i$, are assumed to be non-negative integers. This is not restrictive since any negative term such as $-a_{ij}x_j$ with $a_{ij} > 0$ can be replaced by $a_{ij}\bar{x}_j - a_{ij}$, where $\bar{x}_j = 1 - x_j$, to obtain an equivalent expression with only non-negative coefficients. For example, $-2x_1 - 3x_2 + 2x_3 \ge -2 \Leftrightarrow 2(1-x_1) + 3(1-x_2) + 2x_3 \ge -2+2+3 \Leftrightarrow 2\bar{x}_1 + 3\bar{x}_2 + 2x_3 \ge 3$. The transformed constraint always has a non-negative right-hand side. Otherwise, it would be trivially satisfied and could be safely removed.

Let $X \subseteq \{0, 1\}^n$ denote the set of feasible solutions in the decision space $(\{0, 1\}^n)$, i.e., the set of solutions that satisfy the problem constraints. The vector-valued function $f: X \to \mathbb{Z}_{\ge 0}^p$ maps a feasible solution to an (objective) vector, or point, in the objective space, $\mathbb{Z}_{\ge 0}^p$, where *p* is the number of objectives. The set $Y \subseteq \mathbb{Z}_{\ge 0}^p$ is the image set of *X*, i.e., $Y = f(X)$, which represents the set of feasible vectors, or feasible points. Hence, to each feasible solution $x \in X$ corresponds an objective vector $f(x) = (f_1(x), \dots, f_p(x))$, where $f(x) \in Y$.

Given two points in the objective space, $z, z' \in \mathbb{Z}_{\ge 0}^p$, z is said to weakly dominate $z'$ if $z_k \le z'_k$ for all $k \in \{1, \dots, p\}$, where $z_k$ and $z'_k$ denote the $k$th coordinate of z and $z'$, respectively. This is represented as $z \le z'$. Point z is said to *dominate* $z'$ if $z \le z'$ and $z \ne z'$. This is represented as $z \le z'$. Point z is said to *strictly dominate* $z'$ if $z_k < z'_k$ for all $k \in \{1, \dots, p\}$. This is represented as $z < z'$. The points z and $z'$ are said to be *incomparable*, or *mutually nondominated*, if neither $z \le z'$ nor $z' \le z$ (see, e.g., Ehrgott, 2005). Weak dominance may be extended to point sets. Given $A, A' \subset \mathbb{R}^p$, A is said to weakly dominate $A'$ if, for each $z' \in A'$, there exists $z \in A$ such that $z \le z'$ (Zitzler et al., 2003).

A feasible solution *x* is said to be an *efficient* (or *Pareto-optimal*) solution if there is no other solution $x' \in X$ such that $f(x') \le f(x)$, in which case $f(x)$ is said to be a *nondominated point*. The set of all efficient solutions is known as the *efficient set* (or *Pareto set*). The corresponding image set is known as the *nondominated set* (or *Pareto front*), which

is here represented by $Y_N \subseteq Y$. Given a MOBO formulation, the goal is (typically) to find the set of all nondominated points and an efficient solution corresponding to each such point. If enumerating all nondominated points is not possible (e.g., it takes too long, or $Y_N$ is too large), an alternative is to find a $(1+\epsilon)$-approximation. Algorithms to find such approximation sets exist in the literature, though mostly are either of theoretical nature, or for a specific problem such as the multiobjective knapsack problem (see, e.g., Herzel et al., 2021b; Bazgan et al., 2009). A set $A_\epsilon \subseteq X$ of feasible solutions is called a $(1+\epsilon)$-approximation, if for any feasible solution $x' \in X$, there exists a solution $x \in A_\epsilon$ such that:

$$f_j(x) \le (1 + \epsilon) f_j(x') \text{ for all } j \in \{1, \dots, p\}$$

where $\epsilon \ge 0$. An approximation ratio of an arbitrary point set, $A \subset \mathbb{R}_+^p$, to a given reference point set, $R \subset \mathbb{R}_+^p$ (e.g., the Pareto front, $R = Y_N$), can be calculated from the (multiplicative) $\epsilon$-indicator (Zitzler et al., 2003):

$$I_{\epsilon}(A, R) = \max_{r \in R} \min_{a \in A} \max_{k \in \{1, \dots, p\}} \left\{ \frac{a_k}{r_k} \right\}$$

If $A_\epsilon$ is a $(1+\epsilon)$-approximation then we have $I_\epsilon(f(A_\epsilon), Y_N) \le 1+\epsilon$. A lower bound set, $L \subset \mathbb{R}^p$, i.e., a set of mutually nondominated points that weakly dominates $Y_N$, can be easier to compute than $Y_N$ and provides an upper bound on the approximation ratio of any point set $A \subseteq Y$, i.e., $I_\epsilon(A, L) \ge I_\epsilon(A, Y_N)$ (Zitzler et al., 2008). With a slight abuse of notation, when referring to a solution $x \in X$, for which we explicitly assign values to the variables, we will refer to it as an assignment $v$ and define it using set notation (e.g., $v = \{x_1 = 1, x_2 = 0\}$). We will also extend the notation to state that $v \in X$, and use $f(v)$. Next, we present an example illustrating some of the introduced definitions.

**Example 1.** Let $f(x) = (2x_1 + x_2 + \bar{x}_3, \bar{x}_1 + x_2 + 2x_3)$ denote the vector with two objective functions to minimize such that the set of constraints $\{x_1 + x_2 + x_3 \ge 2\}$ needs to be satisfied. The four feasible solutions correspond to the following assignments $v^1$ with objective vectors $f(v^1)$:

$$v^1 = \{x_1 = 1, x_2 = 1, x_3 = 0\} \text{ with } f(v^1) = (4, 1)$$

$$v^2 = \{x_1 = 1, x_2 = 0, x_3 = 1\} \text{ with } f(v^2) = (2, 2)$$

$$v^3 = \{x_1 = 0, x_2 = 1, x_3 = 1\} \text{ with } f(v^3) = (1, 4)$$

$$v^4 = \{x_1 = 1, x_2 = 1, x_3 = 1\} \text{ with } f(v^4) = (3, 3)$$

In this case, assignments $v^1, v^2$ and $v^3$ correspond to the three efficient solutions of this problem. Hence, $Y_N = \{(1, 4), (2, 2), (4, 1)\}$. Moreover, $A_\epsilon = \{v^2\}$ is a $(1+\epsilon)$-approximation for any $\epsilon \ge 1$. Observe that this set has smaller cardinality than $Y_N$ and that $I_\epsilon(f(A_\epsilon), Y_N) = 2$. Remark also that $v^4$, which is not efficient, also forms a $(1+\epsilon)$-approximation of small cardinality, but for any $\epsilon \ge 2$.

### 2.2. Boolean optimization problems as maximum satisfiability problems

Boolean optimization problems can be encoded as logic-based problems. Section 2.2.1 introduces the Maximum Satisfiability (MaxSAT) problem, defines Minimal Correction Subsets (MCSs), and explains how an optimal solution to a MaxSAT problem can be found by enumerating all MCSs. Section 2.2.2 describes an encoding of a (single-objective) Boolean optimization problem as a MaxSAT problem, which will be later extended to the multi-objective case in Section 3.

#### 2.2.1. Maximum satisfiability and minimal correction subsets

In propositional logic, a formula is defined over a set of Boolean variables $x_1, \dots, x_n \in \{0, 1\}$. In formulas expressed in Conjunctive Normal Form (CNF), a clause represents a disjunction of literals, where each literal is a variable (e.g., $x_1$) or its negation (e.g., $\bar{x}_1$), and the formula represents a conjunction of clauses. For example, $\phi = \{(x_1 \lor \bar{x}_2), (\bar{x}_1 \lor x_2 \lor x_3)\}$ is a CNF formula with two clauses. A literal $x_j$ is satisfied by an assignment $v$ if $v$ assigns 1 to $x_j$ (i.e., $x_j = 1$), whereas a negated literal $\bar{x}_j$ is satisfied by $v$ if $v$ assigns 0 to $x_j$ (i.e., $x_j = 0$). An assignment
is said to satisfy a clause if it satisfies at least one literal in the clause, and it is said to satisfy a formula if it satisfies all clauses in the formula. A complete assignment that satisfies a given formula is also called a feasible assignment or a model. A formula $\phi$ is said to be satisfiable if there exists an assignment that satisfies $\phi$. Otherwise, $\phi$ is said to be unsatisfiable.

**Example 2.** Consider the following CNF formulas $\phi = \{(x_1 \lor \bar{x}_2), (\bar{x}_1 \lor \bar{x}_2)\}$, and $\phi' = \{(x_1 \lor \bar{x}_2), (\bar{x}_1 \lor \bar{x}_2), (x_2)\}$. Formula $\phi$ is satisfiable and $\nu = \{x_1 = 1, x_2 = 0\}$, denoted more compactly as $\nu = (1,0)$, is a model of $\phi$, whereas formula $\phi'$ is unsatisfiable.

Given a CNF formula $\phi$, the Propositional Satisfiability (SAT) problem consists of finding a satisfiable assignment to $\phi$ or prove that such assignment does not exist. The Maximum Satisfiability (MaxSAT) problem is an optimization version of the SAT problem. In general, a (partial unweighted) MaxSAT formula $\phi$ consists of two sets of clauses, the hard ($\phi^h$) and the soft ($\phi^s$) clauses such that $\phi = \phi^h \cup \phi^s$. Assuming that a feasible assignment for $\phi^h$ exists, the goal is to find an assignment that satisfies all hard clauses and as many soft clauses as possible, i.e., minimizes the number of unsatisfied soft clauses. Observe that usually there is no assignment that simultaneously satisfies all hard clauses in $\phi^h$ and all soft clauses in $\phi^s$.

A correction subset of a MaxSAT formula $\phi$ is a subset of soft clauses $C \subseteq \phi^s$ such that $\phi^h \cup (\phi^s \setminus C)$ is satisfiable. Additionally, if $\phi^h \cup (\phi^s \setminus C) \cup \{c\}$ is unsatisfiable for all $c \in C$, then $C$ is a Minimal Correction Subset (MCS). That is, $C$ is minimal in the sense that none of the proper subsets of $C$ is a correction subset.

A complete assignment $\nu$ is said to correspond to the MCS $C$ if $\nu$ satisfies $\phi^h \cup (\phi^s \setminus C)$. Note that since $C$ is an MCS, then $\nu$ does not satisfy any clause in $C$. Given a MaxSAT formula $\phi$, an optimal solution to $\phi$ corresponds to the MCS $C$ with minimum size. Hence, if we were to enumerate all MCSs of $\phi$ (Previti et al., 2017; Grégoire et al., 2018), then we can simply select an MCS of minimum size as an optimal solution for $\phi$.

**Example 3.** Let $\phi^h = \{(\bar{x}_1 \lor \bar{x}_2 \lor \bar{x}_3), (x_1 \lor x_2), (\bar{x}_1 \lor x_2 \lor x_3)\}$ and $\phi^s = \{(\bar{x}_1), (\bar{x}_2), (\bar{x}_3)\}$ define the set of hard and soft clauses of a MaxSAT formula $\phi$. There are two MCSs, $C^1 = \{(\bar{x}_1), (\bar{x}_3)\}$, and $C^2 = \{(\bar{x}_2)\}$. The assignments $\nu^1 = (1,0,1)$, and $\nu^2 = (0,1,0)$ correspond to MCSs $C^1$ and $C^2$, respectively. Therefore, assignment $\nu^2$, corresponding to $C^2$, is an optimal solution for the MaxSAT instance $\phi$.

Note that linear single-objective Boolean optimization problems can be easily formulated as MaxSAT. For instance, the MaxSAT formula in **Example 3** corresponds to the problem of satisfying $\phi^h$ while minimizing $x_1 + x_2 + x_3$. Observe that the cost function is encoded into the soft clauses $\phi^s$. Arbitrary linear objective functions are formulated using a more general version of MaxSAT, where each soft clause has an associated weight representing the cost of not satisfying it. For example, $2x_1 + 3x_2$ is encoded into the soft clauses $(\bar{x}_1)$ and $(\bar{x}_2)$ with weights 2 and 3, respectively. In fact, any 0-1 Integer Linear Programming (0-1 ILP) can be encoded into MaxSAT by encoding the linear constraints into hard clauses and the objective function as soft clauses. The encoding of linear constraints into CNF formulas has been extensively studied (Bailleux and Boufkhad, 2003; Sinz, 2005; Asín et al., 2011; Ogawa et al., 2013; Hölldobler et al., 2012; Bailleux et al., 2009, 2006; Eén and Sörensson, 2006; Abío et al., 2012).

There exist different ways of encoding the objective function(s) into soft clauses, which has implications on how MCSs and (Pareto-)optimal solutions relate. In the above example, the encoding relies on directly using the decision variables present in the objective function as soft clauses (and their coefficients as the corresponding weights). This idea has already been extended to the multi-objective case, for which it was shown that each Pareto-optimal solution of a MOBO problem corresponds to an MCS (Terra-Neves et al., 2017). Hence, all Pareto-optimal solutions can be found by enumerating all MCSs of the resulting

MaxSAT formula. However, this encoding has the disadvantage that not all MCSs correspond to Pareto-optimal solutions, which means that some dominated solutions may be enumerated as well. One way to overcome this issue is to use a different encoding for which there is a one-to-one correspondence between nondominated points and MCSs. The encoding of objective functions using a unary representation has this property, as seen in the next sections.

### 2.2.2. Encoding the objective function(s) using a unary representation

Objective functions can be encoded in the MaxSAT formulation by adding, as soft clauses, new variables that represent (bounds on) the value of the objective function (see, e.g., Soh et al., 2017). Assume that the objective function $f: \{0,1\}^n \to \mathbb{Z}_{\ge 0}$ is lower and upper bounded by $\ell$ and $u$, respectively, i.e., $\ell \le f(\nu) \le u$ for all $\nu \in \{0,1\}^n$. A way of encoding $f(x)$ to CNF with this approach is to include a new set of Boolean variables, $y_d$, and encode the equivalences $y_d = 1 \Leftrightarrow f(x) < d$ for $d \in D \subseteq \{\ell, \ell+1, \dots, u^+\}$ where, throughout the paper, $u^+$ is assumed to be a value greater than $u$ (i.e., $u^+ > u$). As discussed at the end of this section, such equivalences can be encoded with an additional set of hard clauses. Variables $y_\ell$ and $y_{u^+}$ are not strictly required in the encoding, as $y_\ell$ is trivially false and $y_{u^+}$ is trivially true. To make the analysis smoother, we will assume however that $\ell, u^+ \in D$, which is a weak assumption since $\ell$ and $u$ are supposed to be known lower and upper bounds of the objective function value. The precise definition of $D$ depends on whether we aim at solving our discrete optimization problem exactly or approximately. In the first case $D$ must be (a superset of) all values that can be reached by the objective function, e.g., be the set of all integers between $\ell$ and $u^+$ (in this case we assume $u^+ = u+1$). If so, $D$ will be called complete. In the second case, depending on the required degree of approximation, $D$ should be defined as an appropriate subset of all integers between $\ell$ and $u^+$.

With this unary representation of the objective function, a MaxSAT formulation of any single objective Boolean optimization problem is obtained by defining $\{(y_d) \mid d \in D\}$ as the set of soft clauses, with $D \subseteq \{\ell, \ell+1, \dots, u^+\}$ and $\ell, u^+ \in D$, while constraints and the previous equivalences are represented by hard clauses. Then, solving the single objective Boolean optimization problem amounts to finding the (unique) MCS of the MaxSAT formula. Assuming indeed that the optimal objective function value is $d^*$, the corresponding MCS will be of the form $C = \{(y_d) \mid d \in D \text{ and } d < d^*\}$. Observe first that $C$ always exists since it contains at least $y_\ell$. Let $d'$ be the largest index among variables $y_d \in C$ and $d''$ be the smallest index among variables in $D \setminus C$. Then, any feasible assignment $\nu$ corresponding to $C$ is such that $d' \le f(\nu) < d''$. Hence, the solutions corresponding to $C$ are within an approximation guarantee provided by $d'$ and $d''$. In the particular case where $D$ is complete, any solution $\nu$ corresponding to $C$ is optimal with objective function value $f(\nu) = d'$. Consider the following example where $D$ is complete.

**Example 4.** Let $f(x) = 3x_1 + 2x_2 + 2x_3$ be the function to minimize such that $\{x_1 + x_2 \ge 1, \bar{x}_2 + x_3 \ge 1\}$ must be satisfied. Observing from the definition of $f$ that it can only take values in $D' = \{0,2,3,4,5,7\}$ we can define a complete domain $D = D' \cup \{8\}$ and variables $y_d$, for $d \in D$. Let CNF($f(x)$) denote the CNF formula encoding that $y_d = 1 \Leftrightarrow f(x) < d$. Hence, one can define a MaxSAT formula $\phi = \phi^h \cup \phi^s$ such that $\phi^h = \{(x_1 \lor x_2), (\bar{x}_2 \lor x_3)\} \cup \text{CNF}(f(x))$ and $\phi^s = \{(y_0), (y_2), (y_3), (y_4), (y_5), (y_7), (y_8)\}$. Note that the optimal solution of $\phi$ corresponds to the optimal solution of the original problem. Observe that $\phi$ only has a single MCS corresponding to the optimal solution. In this case, the MCS is $\{(y_0), (y_2), (y_3)\}$ since the optimal solution has a cost of 3.

It remains to explain how to actually encode (into hard clauses) the equivalences $y_d = 1 \Leftrightarrow f(x) < d$, for all $d \in D$. Since, in this paper, $f(x)$ is a linear expression defined over a set of Boolean variables, it
is possible to develop solvers that natively handle these expressions. However, many Boolean solvers perform an encoding of linear Pseudo-Boolean expressions into CNF in order to take advantage of effective state of the art SAT solvers. It is known that such expressions can be represented into CNF using a polynomial size encoding (Bailleux et al., 2009). However, it has been observed that some unary representations of the value of linear expressions are better handled by SAT solvers (Bacchus et al., 2019) even though an exponential-size encoding is required in the worst case. Hence, in this paper, the actual encoding of the aforementioned equivalences is done only after encoding the value of the objective function $f(x)$ into CNF. First, we use an encoding using selection networks (Karpinski and Piotrów, 2019) that has been shown to be more compact than alternative encodings for linear expressions with Boolean variables. Next, the selection networks encoding is extended such that a unary encoding is produced where $y_d = 1 \iff f(x) < d$. Moreover, besides adding the clauses to encode the equivalence $y_d = 1 \iff f(x) < d$ for all $d \in D$, an additional order encoding (Tamura and Banbara, 2008) is also added such that $y_d = 1 \Rightarrow y_{d+1} = 1$ and $y_d = 0 \Rightarrow y_{d-1} = 0$. We refer the interested reader to the literature for further details on this and other encodings (Roussel and Manquinho, 2009; Joshi et al., 2015; Karpinski and Piotrów, 2019, 2020). Finally, recall that in MOBO there are several objective functions. Therefore, each objective function is represented using the described encoding into CNF.

### 3. Pareto MCSs

In this section we investigate the relation between MCSs and the nondominated points of a MOBO problem (see Section 3.1). In the special case where the domains $D_k$ of the objective functions are complete, we prove that MCSs and the nondominated points are in one to one correspondence. A similar result was obtained in Soh et al. (2017), in a different setting, using so-called P-minimal models, instead of MCSs. In the general case, the domains are not complete, meaning that some (or possibly all) values taken by objective function $f_k$ are not present in $D_k$ for $k \in \{1, ..., p\}$. Then, we show that MCSs provide some guarantee on the quality of the associated solutions. This lays the foundation for the generation of $(1+\epsilon)$-approximation sets investigated in Section 3.2.

In the following, let $\ell_k \in \mathbb{Z}_{\ge 0}$ and $u_k \in \mathbb{Z}_{\ge 0}$ denote the lower and upper bound values of the $k$th objective function of a MOBO problem, i.e., $\ell_k \le f_k(x) \le u_k$, for all $k \in \{1, ..., p\}$, and let $D_k \subseteq \{\ell_k, \ell_k + 1, ..., u_k, u_k^+\}$ and $\ell_k, u_k^+ \in D_k$ where $u_k^+ > u_k$. Moreover, let the variables $y_{k,d}$, where $k \in \{1, ..., p\}$ and $d \in D_k$, encode the values of the $k$th objective function (see Section 2.2.2), i.e., $y_{k,d} = 1 \iff f_k(x) < d$. Problem 1 is encoded as a MaxSAT problem, where $\phi^h$ is a set of hard clauses encoding (in CNF) the problem constraints, and encoding each of the objective functions and their corresponding values in $D_1, ..., D_p$. Hence, $\phi^h$ enforces $x \in X$ and $y_{k,d} = 1 \iff f_k(x) < d$, for all $k \in \{1, ..., p\}$ and $d \in D_k$. Finally, the set of soft clauses, $\phi^s$, consists of the set of all literals $y_{k,d}$, i.e., $\phi^s = \{(y_{k,d}) | k = 1, ..., p \text{ and } d \in D_k\}$.

We first give two preliminary simple results related to MCSs.

**Lemma 1.** Let $C \subset \phi^s$ be an MCS for our MaxSAT formulation of a MOBO problem.

(i) If $y_{k,d} \in C$, then $y_{k,e} \in C$ for $e \le d$.

(ii) If $y_{k,d} \notin C$, then $y_{k,e} \notin C$ for $e \ge d$.

**Proof.**

(i) Since $y_{k,d} \in C$, the solution $x$ associated with the current model is such that $f_k(x) \ge d$. Then we necessarily have $f_k(x) \ge e$ for $e \le d$.

(ii) Since $y_{k,d} \notin C$, the solution $x$ associated with the current model is such that $f_k(x) < d$. Then we necessarily have $f_k(x) < e$ for $e \ge d$. $\square$

Let $D = D_1 \times \dots \times D_p$. To any MCS $C$ we associate the two following points:

* its representative point $r \in D$ whose components are $r_k = \max\{d \in D_k | y_{k,d} \in C\}$

* its successor point $r' \in D$ whose components are $r'_k = \min\{d \in D_k | y_{k,d} \notin C\}$, for all $k \in \{1, ..., p\}$.

The point $r$ is named representative due to its relevance, and $r'$ is named successor due to its connection to $r$ as the components $r'_k$ correspond to the successor values of $r_k$ in $D_k$. The following result shows that $r$ and $r'$ correspond to lower and upper bound points of the image of the solutions associated to $C$.

**Lemma 2.** Let $r$ and $r'$ be the representative and successor point of an MCS $C$. Then any assignment $v \in \{0, 1\}^n$ associated to $C$ is such that $r \le f(v) < r'$. Moreover, if domains $D_k$ are complete, for all $k \in \{1, ..., p\}$, then $f(v) = r$.

**Proof.** By definition, $y_{k,r_k} \in C$ implying $f_k(v) \ge r_k$, and $y_{k,r'_k} \notin C$ implying $f_k(v) < r'_k$, for all $k \in \{1, ..., p\}$. When domains $D_k$ are complete, there is no feasible solution $x \in X$, such that $r_k < f_k(x) < r'_k$, implying the equality. $\square$

### 3.1. Finding the nondominated set

Assuming that $D_k$ is complete and $\ell_k, u_k^+ \in D_k$, for all $k \in \{1, ..., p\}$, we show that each MCS corresponds to a nondominated point, and vice-versa. It follows that all nondominated points can be obtained by enumerating all MCSs. We assume that **Problem 1** is feasible. Therefore, all hard constraints can be satisfied. MCSs are computed only with respect to variables $y_{k,d}$ for all $k \in \{1, ..., p\}$ and $d \in D_k$.

**Theorem 1.** Consider a MOBO problem, with feasible point set $Y$, stated using our MaxSAT formulation with complete domains $D_k$ for all $k \in \{1, ..., p\}$. Then nondominated points in $Y_N$ and MCSs of the MaxSAT formulation are in one to one correspondence.

**Proof.** Domains $D_k$ being complete, any point in $Y$ is of the form $(d_1, ..., d_k, ..., d_p)$, where $d_k \in D_k$, for all $k \in \{1, ..., p\}$. We prove that a point $d = (d_1, ..., d_k, ..., d_p)$ belongs to $Y_N$ if and only if

$$C = \{y_1, \ell_1, \dots, y_{1,d_1}, \dots, y_{k, \ell_k}, \dots, y_{k, d_k}, \dots, y_{p, \ell_p}, \dots, y_{p, d_p}\}$$

is an MCS of the MaxSAT formulation. Let $r = d$ and $r'$ be the representative and successor points of $C$, respectively.

⇔: By **Lemma 2**, point $d$, which is the representative point of $C$, belongs to $Y$. Moreover, $d \in Y_N$, since the existence of a feasible solution $x'$ such that $f(x') \le d$, would require $f_k(x') \le d_k$ for all $k \in \{1, ..., p\}$ and $f_q(x') < d_q$ for some $q \in \{1, ..., p\}$, meaning that we could remove variable $y_{q,d_q}$ from $C$, contradicting that $C$ is an MCS. ⇒: Since $d$ is feasible, $C$ is a correction subset. Since $d$ is nondominated, we cannot remove from $C$ any variable $y_{q,d_q}$ for some $q \in \{1, ..., p\}$. By **Lemma 1(i)**, we cannot remove any variable from $C$, which shows that $C$ is an MCS. $\square$

As seen in **Theorem 1**, the representative points of MCSs always correspond to nondominated points when $D_k$ is complete for all $k \in \{1, ..., p\}$. A similar reasoning can be used to show that, in the more general case of $D_k$ defined as an arbitrary subset of $\{\ell_k, \ell_{k+1}, \dots, u_k, u_k^+\}$ where $\ell_k, u_k^+ \in D_k$, the representative points of the MCSs are still mutually nondominated, though are not necessarily feasible, and the correspondence between nondominated points and MCSs becomes one-to-many. That is, an MCS is associated to each nondominated point, but more than one nondominated points may be associated to an MCS (and some dominated points may be as well). The next section explores this more general definition of $D_k$.
<figure><img src="image_4.png" /><figcaption>Fig. 1. An example of a bi-objective problem with different settings of $D_1$ and $D_2$, depending on $\epsilon$. Large and small circles represent the images of feasible and infeasible solutions, respectively. The vertical and horizontal lines represent the values in $D_1$ and $D_2$, respectively, and the shadowed regions illustrate the intervals associated to $D_1$ intersected with those associated to $D_2$. The crosses depict the representative points of the MCSs.</figcaption></figure>

## 3.2. Finding a $(1 + \epsilon)$-approximation set

Having variables $y_{k,d}$ represent, in an exact way, all possible values of the corresponding objective function, allows the enumeration of all nondominated points through MCS enumeration. Although it may be desirable to find the Pareto front, this may require large encodings, i.e., a large set of hard and soft clauses, and the resulting formula may still be difficult to solve. Alternatively, finding a $(1 + \epsilon)$-approximation set may be enough, and even helpful if used as an intermediate step in the search for the Pareto front. In this section we describe two encodings related to the one used in Section 3.1, but which are smaller and (potentially) easier to solve, and for which enumerating all MCSs is equivalent to finding a $(1 + \epsilon)$-approximation set. One of the encodings consists in considering a smaller set of soft clauses (see Section 3.2.1), and in the other the variables $y_{k,d}$ encode an approximate version of the objective function values (see Section 3.2.2).

### 3.2.1. Interval-based approximation

We assume here that each objective function $f_k$ takes positive integer values. An approximate version of the method, with a priori guarantee, can be obtained by modifying the domains $D_k$ of each objective function $f_k$. Assume for instance that we want to generate a $(1 + \epsilon)$-approximation $A_\epsilon \subseteq X$, with $\epsilon \ge 0$. In this case, we define $D_k$, for each $k \in \{1, \dots, p\}$, as follows:

$$D_k = \{d_{k,1}, \dots, d_{k,u}\}$$

with $d_{k,1} = \ell_k$, $d_{k,i} = \max\{d_{k,i-1} + 1, [(1 + \epsilon) \cdot d_{k,i-1}]\}$ for $i \in \{2, \dots, u\}$ where $u$ is the smallest integer such that $d_{k,u} > u_k$. Each value $d_{k,i} \in D_k$, for $i \in \{1, \dots, u-1\}$, defines the beginning of a new interval which is upper bounded by the successor value in $D_k$, i.e., $d_{k,i+1}$. This definition of $D_k$ admits $\epsilon = 0$, in which case $D_k = \{\ell_k, \dots, u_k + 1\}$, which corresponds to the exact version.

Consider the following bi-objective Boolean optimization problem, which is illustrated in Fig. 1(a):

$$ 
\begin{array}{l@{\qquad}l}
\min f_1(x) & = 3\bar{x}_1 + \bar{x}_2 + 2\bar{x}_3 + x_4 + x_5 + 2 \\
\min f_2(x) & = x_1 + 2x_2 + 16x_3 + x_5 + 2 \\
\text{subject to:} & \\
& x_1 + x_2 + x_3 + x_4 + x_5 \ge 3 \\
& x_1 + x_2 + x_3 + x_4 + x_5 \le 4 \\
& x \in \{0, 1\}^n
\end{array}
\qquad (1) 
$$

Fig. 1(a) shows the image $f(v)$ of all assignments $v \in \{0, 1\}^n$, including infeasible ones (small circles), i.e., that do not satisfy the two cardinality constraints. The nondominated points are colored black. In this case, $f_1$ takes values between $\ell_1 = 2$ and $u_1 = 10$, and $f_2$ takes values between $\ell_2 = 2$ and $u_2 = 22$.

In the examples for the particular case of $\epsilon = 0$, we will consider $D_k = \{f_k(x) | x \in \{0, 1\}^n\} \cup \{u_k + 1\}$ which is complete. In Fig. 1(a), these values of $D_k$ are highlighted with dashed lines. For $\epsilon = 0.5$, we have $D_1 = \{2, 3, 4, 5, 7, 9, 13\}$ and $D_2 = \{2, 3, 4, 5, 7, 9, 12, 16, 24\}$. For $\epsilon = 1$, $D_1 = \{2, 4, 8, 16\}$ and $D_2 = \{2, 4, 8, 16, 32\}$, as illustrated by the dashed lines in Fig. 1(b). The values $16 \in D_1$ and $32 \in D_2$ were omitted for convenience. Note that, in Fig. 1(b), there are empty intervals that could be excluded from $D_1$ and $D_2$ while still ensuring the approximation ratio of 2. In the example, setting $D_1$ and $D_2$ to $\{2, 4, 8\}$ and $\{2, 4, 19\}$, respectively, where $D_k \subseteq \{f_k(x) | x \in \{0, 1\}^n\}$ for $k = 1, 2$, would ensure such approximation ratio. However, as the image set of $f_k(x)$ may be much larger than $D_k$, ensuring that $D_k$ is a subset of all values that can be reached by the $k$th objective function may be time consuming and not be advantageous in general, and be advantageous only to particular instances with particular values of $\epsilon$.

For $\epsilon = 0$, each MCS C corresponds to a single nondominated point, which is its representative point. This is not necessarily the case for $\epsilon > 0$. Firstly, although each nondominated point corresponds to a single MCS, to each MCS may correspond multiple feasible points, even dominated ones, all of which are weakly dominated by the representative point of the MCS and are within a $(1 + \epsilon)$ ratio from it. Secondly, the representative point of an MCS does not correspond, in general, to a feasible point. Let $D = D_1 \times \dots \times D_p$, let $r, r' \in D$ be such that $r'_k$ is the successor value of $r_k$ in $D_k$ for all $k = 1, \dots, p$, and let the feasible point $s \in Y$ be such that $r \le s < r'$. In Fig. 1(b), $r$ and $r'$ are the lower-left and upper-right corners of a rectangular orange region, respectively, and $s$ is any point in that region. Any such point $s$ is associated to a model where $y_{1,r_1} = 0, \dots, y_{p,r_p} = 0$ and $y_{1,r'} = 1, \dots, y_{p,r_p'} = 1$. Therefore, $s$ is weakly dominated by $r$ and is within a $(1 + \epsilon)$ ratio from $r$. If $r$ is a representative point of an MCS C, then such feasible points $s$ are the only points corresponding to C. In Fig. 1(b), the points $(2, 21)$, $(3, 21)$, and $(3, 22)$ correspond to the MCS represented by point $(2, 16)$ (orange cross), as their corresponding models include $y_{1,2} = 0, y_{2,16} = 0, y_{1,4} = 1$, and $y_{2,32} = 1$.

Since the set of soft clauses, $\phi^s$, only includes variables $y_{k,d}$ for which $d \in D_k$, enumerating all MCSs and finding a solution associated to each one of them is equivalent to finding a $(1 + \epsilon)$-approximation
set. Additionally, since each nondominated point $s \in Y_N$ is associated to an MCS whose representative point $r \in D$ weakly dominates $s$, the set of the representative points of all MCSs is a lower bound set. In the example of Fig. 1(b), the lower bound set is $\{(2, 16), (4, 4), (8, 2)\}$ which is represented by the set of crosses. Hence, enumerating all MCSs not only provides an approximation set $A_\epsilon \subseteq X$, but also a lower bound set $\mathcal{L} \subset \mathbb{Z}_{\ge 0}^p$, for which $I_\epsilon(f(A_\epsilon), Y_N) \le I_\epsilon(f(A_\epsilon), \mathcal{L}) \le (1+\epsilon)$.

### 3.2.2. Coefficient-based approximation

An alternative approximate version, also with a priori guarantee, can be obtained by modifying the coefficients $w_1, \dots, w_n$ of the objective function $f_k(x) = \sum_{j=1}^n w_j x_j$, where $k \in \{1, \dots, p\}$, and $w \in \mathbb{Z}_{\ge 0}^n$. Consider $W_k$ defined as $W_k = \{\hat{w}_1, \dots, \hat{w}_u\}$ with $\hat{w}_1 = \min\{w_1, \dots, w_n\}$, and $\hat{w}_l = \max\{\hat{w}_{l-1} + 1, [(1+\epsilon) \cdot \hat{w}_{l-1}]\}$ for $i \in \{2, \dots, u\}$ where $u$ is the smallest integer such that $\hat{w}_u \ge \max\{w_1, \dots, w_n\}$. Let $w'_1, \dots, w'_n \in W_k$ denote the modified coefficients where $w'_j = \max\{\hat{w}_j \in W_k | \hat{w}_j \le w_j\}$.

**Lemma 3.** Given the function $f_k(x) = \sum_{j=1}^n w_j x_j$ where $w \in \mathbb{Z}_{\ge 0}^n$ and $x \in \{0, 1\}^n$, consider, for any $\epsilon \ge 0$, the function $f'_k: \{0, 1\}^n \to \mathbb{Z}_{\ge 0}^n$ such that $f'_k(x) = \sum_{j=1}^n w'_j x_j$ where $w'_j = \max\{\hat{w}_j \in W_k | \hat{w}_j \le w_j\}$. Then we have:

$$f_k(x) \le (1 + \epsilon) f'_k(x) \text{ for all } x \in \{0, 1\}^n$$

**Proof.** The lemma is trivially true for $\epsilon = 0$, as $f(x) = f'(x)$. For $\epsilon > 0$, a slightly stronger version of the lemma, with a strict inequality, will be proved. Assume by contradiction that there exists $x \in \{0, 1\}^n$ such that $f_k(x) \ge (1+\epsilon)f'_k(x)$, i.e., $\sum_{j=1}^n w_j x_j \ge \sum_{j=1}^n (1+\epsilon)w'_j x_j$. Then, there must exist an index $j \in \{1, \dots, n\}$ such that

$$(1 + \epsilon) w'_{j} \le w_{j} \qquad (2)$$

Modified weight $w'_j$ corresponds to a weight $\hat{w}_j$ where either $i = u$ or $i < u$. If $i = u$ then $w'_j = w_j(=\hat{w}_u)$ in which case (2) is not possible, leading to a contradiction. If $i < u$ then we have $\hat{w}_i = w'_j \le w_j < \hat{w}_{i+1}$. Then, by definition, $\hat{w}_{i+1}$ is either $[(1+\epsilon)\hat{w}_i]$ or $\hat{w}_{i+1}$. In the first case we get $w_j < [(1+\epsilon)\hat{w}_j] = [(1+\epsilon)w'_j] \le (1+\epsilon)w'_j$ contradicting (2). The second case implies that $w_j = \hat{w}_j = w'_j$ because $w_j$ must be an integer such that $\hat{w}_j \le w_j < \hat{w}_{i+1} = \hat{w}_{i+1}$ also leading to a contradiction, which completes the proof. $\square$

Consider, for example, the objective function $f_1(x) = 2x_1 + 3x_2 + 5x_3 + 7x_4$, and let $\epsilon = 1$, then we have $W_1 = \{2, 4, 8\}$ leading to $f'_1(x) = 2x_1 + 2x_2 + 4x_3 + 4x_4$, and any solution $x \in X$ satisfies $f_1(x) \le 2f'_1(x)$. Let us now prove that a $(1+\epsilon)$-approximation of Problem 1 can be obtained by finding the set of efficient solutions for $f' = (f'_1, \dots, f'_p)$ where functions $f'_k$ are defined as in Lemma 3 for $k = 1, \dots, p$.

**Theorem 2.** Let $A_\epsilon \subseteq X$ be the set of efficient solutions for $f' = (f'_1, \dots, f'_p)$ where functions $f'_k$ are defined as in Lemma 3 for $k = 1, \dots, p$. Then, $A_\epsilon$ is a $(1+\epsilon)$-approximation of Problem 1.

**Proof.** Since $A_\epsilon$ is the set of efficient solutions for $f'$, we have that for any $x \in X$ there exists $x' \in A_\epsilon$ such that $f'(x') \le f'(x)$ and thus such that $(1+\epsilon)f'(x') \le (1+\epsilon)f'(x)$.

Then, by Lemma 3 applied to $x'$, we have $f(x') \le (1+\epsilon)f'(x')$. Finally, by definition of $f$ and $f'$, which ensures that for any $x \in X$ we have $f'(x) \le f(x)$, we have $(1+\epsilon)f'(x) \le (1+\epsilon)f(x)$.

It follows that for any $x \in X$ there exists $x' \in A_\epsilon$ such that:

$$f(x') \le (1 + \epsilon) f'(x') \le (1 + \epsilon) f'(x) \le (1 + \epsilon) f(x)$$

establishing that $A_\epsilon$ is a $(1+\epsilon)$-approximation for $f$. $\square$

An implication of Theorem 2 is that $f'(A_\epsilon)$ is a lower bound set of Problem 1. Moreover, note that the lemma and the theorem are general in the sense that they do not depend on the method used to find the approximation, and are not even restricted to problems with

linear constraints. In the particular case of our MaxSAT formulation for MOBO, let $y_{k,d}$ be defined for all $d \in D'_k$ and $k = 1, \dots, p$, where $D'_k$ is complete for $f'_k$, i.e., $f'_k(X) \subseteq D'_k$. Then, by Theorems 1 and 2, the Pareto front for the approximate problem can be obtained by enumerating all MCSs given $\phi^s = \{(y_{k,d}) | k = 1, \dots, p \text{ and } d \in D'_k\}$ which corresponds to a $(1+\epsilon)$-approximation set, $A_\epsilon$, for the original problem.

## 4. Computing $(1+\epsilon)$-approximation sets and tighter approximation ratios

This section proposes the algorithms to solve MOBO problems through the corresponding MaxSAT formulation, and based on the two approximation versions of the unary encodings of the objective function values (see Section 3). Such algorithms are guaranteed to find a $(1+\epsilon)$-approximation set, for any setting of $\epsilon \ge 0$. The base algorithm, which is independent from the approximation version, is described in Section 4.1, and then, exact algorithms based on iterative calls to the base algorithm are proposed in Section 4.2. Section 4.3 discusses the advantages and disadvantages of each of the two approximation versions.

### 4.1. $(1+\epsilon)$-Approximation and lower bound sets

The previous section detailed the ideas of how to encode the values of the objective functions of an optimization problem in such a way that any MCS enumeration algorithm can be straightforwardly used to find a $(1+\epsilon)$-approximation set, $A_\epsilon \subseteq X$, where $\epsilon \ge 0$. This includes the case where $\epsilon = 0$, for which the approximation set is the Pareto front. The algorithm does not have to be aware of the approximation ratio and not even which of the two approximation methods is used. An advantage of these approximation methods, and consequently of the algorithm, is to intrinsically provide a lower bound set, $\mathcal{L} \subset \mathbb{Z}_{\ge 0}^p$, within a ratio lower than or equal to $(1+\epsilon)$. Hence, in practice, the algorithm ensures an approximation ratio tighter than $(1+\epsilon)$, which is upper bounded by $I_\epsilon(f(A_\epsilon), \mathcal{L})$, i.e., $I_\epsilon(f(A_\epsilon), Y_N) \le I_\epsilon(f(A_\epsilon), \mathcal{L}) \le (1+\epsilon)$.

This section describes the MCS enumeration algorithm and gives examples of its application to three scenarios: (1) with $\epsilon = 0$; (2) with $\epsilon > 0$ using the interval-based approximation; and (3) with $\epsilon > 0$ using the coefficient-based approximation.

#### 4.1.1. Algorithm

Algorithm 1, named MCS-Approx, shows the pseudocode of the algorithm based on MCS enumeration to compute a $(1+\epsilon)$-approximation set $A_\epsilon$ (stored in $\mathcal{A}$), as well as a lower bound set $\mathcal{L}$. In this case, $\mathcal{A} \subseteq X \times f(X)$ is a set of pairs, where each element $(v, f(v)) \in \mathcal{A}$ represents a feasible solution $v$ and the corresponding image under $f$, $f(v)$. For simplicity, $\mathcal{A}$ will be loosely referred as a $(1+\epsilon)$-approximation set. The algorithm works for any of the two approximation versions presented in the previous section. If $\epsilon = 0$, then Algorithm 1 computes the Pareto front, and in a such case $\mathcal{L}$ also corresponds to the Pareto front.

The input of Algorithm 1 is the following. Let $f$ be the vector of objective functions and $f'$ be the modified vector as defined in Lemma 3 if the coefficient-based approximation is used, otherwise $f' = f$. We assume that $f'_k$ is lower and upper bounded by $l'_k$ and $u'_k$, respectively, for all $k \in \{1, \dots, p\}$. Let $\phi^h$ be the set of hard clauses that encode the constraints of Problem 1 and the approximate objective functions $f'_k$ (as described in Section 2.2.2) for all $k \in \{1, \dots, p\}$. The set $D_k$ denotes the domain of $f'_k$. In the case of the interval-based approximation (see Section 3.2.1), $D_k$ depends on $\epsilon$, the greater $\epsilon$ is, the smaller is $|D_k|$. For the coefficient-based approximation (see Section 3.2.2), $D_k$ is complete with respect to $f'_k$, e.g., $D_k = \{l'_k, \dots, u'_k + 1\}$.

First, Algorithm 1 verifies if $\phi^h$ is satisfiable. If so, the set of soft clauses $\phi^s$ is created by adding a unary clause, $(y_{k,d})$, for each of the values $d \in D_k$ for all $k = 1, \dots, p$ (line 2). Then, while there are MCSs to
<figure><img src="image_5.png" /><figcaption>Fig. 2. Example of the approximation sets returned by Algorithm 1 for the problem in (3), for different settings of $\epsilon$ and different approximation methods.</figcaption></figure>

**Algorithm 1:** MCS-Approx algorithm: Find a $(1 + \epsilon)$-approximation set $(\epsilon \ge 0)$, where $\phi^h$ encodes the problem constraints and the unary representation of the approximate objective functions, $f' = (f'_1, \dots, f'_p)$, whose domains are given in $D = (D_1, \dots, D_p)$.

**Input:** $f, f', \phi^h, D$
**Output:** $A$ $(1 + \epsilon)$-approximation set $A$ and a lower bound set $\mathcal{L}$
1 &nbsp;&nbsp;&nbsp;&nbsp; if $\text{SAT}(\phi^h)$ then
2 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\phi^s \leftarrow \{(y_{k,d} | k = 1, \dots, p \text{ and } d \in D_k)\}$;
3 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $A \leftarrow \mathcal{L} \leftarrow \{\}$;
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{st} \leftarrow \text{true}$;
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **while** ($\text{st} = \text{true}$) **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ($\text{st}, v, C) \leftarrow \text{MCS}(\phi^h, \phi^s)$; // MCS $C$ and a corresponding assignment $v$ ($\text{if st is true}$)
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** ($\text{st} = \text{true}$) **then**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $A \leftarrow A \cup \{(v, f(v))\}$; // $f(v)$ is a feasible point of Problem 1
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** ($k = 1$) **to** ($p$) **do**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $r_k \leftarrow \max\{d \in D_k | y_{k,d} \in C\}$; // $r$ is the representative point of the MCS $C$
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{L} \leftarrow \mathcal{L} \cup \{(r_1, \dots, r_p)\}$; // $(r_1, \dots, r_p) \subseteq f'(v) \subseteq f(v)$ holds
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $c \leftarrow y_{1,r_1} \lor \dots \lor y_{p,r_p}$;
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\phi^h \leftarrow \phi^h \cup \{c\}$; // Block the MCS (i.e., the region weakly dominated by $(r_1, \dots, r_p)$)
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **return** $(A, \mathcal{L})$; // Returns the approximation set and the lower bound set

find, the algorithm iteratively: (i) finds an MCS $C$ and a corresponding (optimal) assignment $v$ (line 6); (ii) updates the set of solutions found by adding the assignment $v$ and the corresponding image in objective space $f(v)$ to $A$ (line 8); (iii) updates the lower bound set by adding the representative point $r$ of $C$ (line 11); (iv) blocks the MCS $C$ (line 13). This blocking clause ensures that the next MCS found is such that there is a $q \in \{1, \dots, p\}$ such that $y_{q,r_q} = 1$. Hence, it prevents the algorithm from finding a solution weakly dominated by $r$ and, therefore, prevents it from visiting the MCS $C$ again. Note that, typically, all literals in $C$ would have been added to the clause, but since $y_{k,j} = 1$ implies $y_{k,j} = 1$ for $i < j$ then the remaining literals are redundant.

### 4.1.2. Example

Fig. 2 shows examples of approximation sets returned by Algorithm 1, depending on the algorithms settings, i.e., the chosen approximation ratio and the approximation method, for the following problem:

$$ \left\{ \begin{array}{l} \min f_1(x) = 3x_1 + 3\bar{x}_2 + x_3 + 2x_4 + 1 \\ \min f_2(x) = 4\bar{x}_1 + 5\bar{x}_2 + 5\bar{x}_3 + 7\bar{x}_4 + 1 \\ \text{subject to } x \in \{0, 1\}^n \end{array} \right. \qquad (3) $$

For the sake of simplicity, this unconstrained problem is used as example. Nevertheless, the reasoning used next is the same for constrained problems since the MCS algorithm ensures that only feasible solutions are computed.

Fig. 2(a) shows all feasible points and, in black, the nondominated points which are all returned by the exact version. In such a case, the algorithm input is $\epsilon = 0$, $f = f'$, and $D_1 = \{1, \dots, 11\}$ and $D_2 = \{1, 5, 6, 8, 10, 11, 12, 13, 15, 17, 18, 22, 23\}$ which are complete, i.e., $f_1(X) \subset D_1$ and $f_2(X) \subset D_2$. Algorithm 1 finds a new nondominated point per iteration. For example, in one of the iterations, line 6 returns the MCS $C = \{y_{1,1}, \dots, y_{1,4}, y_{2,1}, \dots, y_{2,10}\}$, the (only) corresponding model $v = (0, 0, 1, 1)$ and $st = \text{"true"}$. The latter indicates that a new MCS was found ($\phi^h$ is still satisfiable). The corresponding assignment $v$ and $f(v) = (4, 10)$ are added to $A$ in line 8. The representative point of $C$, which in this case is $(4, 10)$, is computed and added to the lower bound set in lines 9–11. To prevent the algorithm from returning any other solution weakly dominated by $r$, lines 12–13 block the MCS $C$ by constructing and adding the clause $(y_{1,4} \lor y_{2,10})$ to the current set of hard clauses, $\phi^h$. Hence, any subsequent solution found must be strictly lower than 4 in objective 1 or strictly lower than 10 in objective 2. The algorithm stops when $st = \text{"false"}$ is returned in line 6 which indicates that there are no more MCSs and, consequently, no more feasible points to be found. In the example, this will happen only after all 6 nondominated points are found: $Y_N = \{(1, 22), (2, 17), (3, 15), (4, 10), (7, 5), (10, 1)\}$.

Fig. 2(b) shows an example of a 2-approximation set (black dots), and the corresponding lower bound set (orange crosses) returned by Algorithm 1 considering the interval-based approximation. The gray dots with black outline represent other points that could have also been returned. In such a case, the algorithm input is $\epsilon = 1$, $f = f'$, and $D_1 = \{1, 2, 4, 8, 16\}$, $D_2 = \{1, 2, 4, 8, 16, 32\}$. The shaded regions show the intersection of the intervals associated to $D_1$ and $D_2$, which contain all feasible points. In this case, as there are only four MCSs, Algorithm 1 executes four iterations, and in each one it finds a solution (which may not be an efficient solution) within the $(1 + \epsilon)$ ratio, and a point of the lower bound set that dominates it. For example, when the MCS $C = \{y_{1,1}, y_{1,2}, y_{1,4}, y_{2,1}, y_{2,2}, y_{2,4}\}$ is found, only the solutions mapping to either point (7, 5) or point (7, 6) correspond to $C$. In the example, we assume that the solution mapping to the latter point is the one returned in line 6. Then, the point added to the lower bound set is (4, 4) because
this is the representative point of C, which is computed in lines 9–10. By blocking the MCS C in line 13 with the clause $(y_{1,4} \lor y_{2,4})$, solutions mapping to a point weakly dominated by (4, 4) become infeasible under $\phi^h$. In the example, the algorithm returns in $\mathcal{A}$ the approximation set $A_\epsilon$ and the corresponding mapping to $f(A_\epsilon) = \{(1, 22), (3, 15), (7, 6), (10, 1)\}$, and the lower bound set $\mathcal{L} = \{(1, 16), (2, 8), (4, 4), (8, 1)\}$, which provide a *posteriori* guaranteed approximation ratio of $I_\epsilon(A_\epsilon, \mathcal{L}) = 1.875 \ge I_\epsilon(A_\epsilon, Y_N) = 1.5$.

Fig. 2(c) shows an example of a 2-approximation set (black dots) returned by Algorithm 1 considering the coefficient-based approximation. The gray dots with black outline represent other points that could have also been returned. In this approximation version with $\epsilon = 1$, the approximate objective functions are: $f'_1(x) = 2x_1 + 2x_2 + 1x_3 + 2x_4 + 1$ and $f'_2(x) = 4\bar{x}_1 + 4\bar{x}_2 + 4\bar{x}_3 + 4\bar{x}_4 + 1$.

In this case, $D_1 = \{1, ..., 9\}$, $D_2 = \{1, 5, 9, 13, 17, 18\}$ are complete with respect to $f'_1$ and $f'_2$, respectively (instead of $f_1$ and $f_2$). In Fig. 2(c), the arrows represent, for each solution $x \in \{0, 1\}^n$, the mapping of $f(x)$ to $f'(x)$, and the blue squares represent all feasible points of $f'$ whereas the darkest ones represent the nondominated points under $f'$. For each of these dark blue squares, i.e., for each nondominated point under $f'$, a solution $v$ mapping into it and the corresponding evaluation under $f$ (the original objective functions) are added to $\mathcal{A}$. The representative point, which in this case must be $f'(v)$, is added to the lower bound set. Blocking the MCS in this case is equivalent to blocking the region weakly dominated by $f'(v)$. In the example, the algorithm returns in $\mathcal{A}$ the approximation set $A_\epsilon$ and the corresponding mapping to $f(A_\epsilon) = \{(1, 22), (2, 17), (5, 13), (7, 6), (10, 1)\}$, and the lower bound set $\mathcal{L} = \{(1, 17), (2, 13), (4, 9), (6, 5), (8, 1)\}$, i.e., the set of black points, and the set of dark-blue squares, respectively. In this case, the *posteriori* guaranteed approximation ratio is $I_\epsilon(A_\epsilon, \mathcal{L}) \approx 1.4 \ge I_\epsilon(A_\epsilon, Y_N) = 1.3$.

### 4.1.3. Enumerating the efficient set

Although Algorithm 1 only finds one solution for each feasible point returned in $\mathcal{A}$, the algorithm could be adapted to find all solutions mapping to each such point (assuming $\epsilon = 0$). This is achieved by replacing the blocking clause in line 12 by $p$ clauses blocking, in the objective space, the region dominated by the representative point $r$ of the last MCS found, and a clause blocking, in the decision space, the assignment $v$. Thus, the representative point $\hat{r}$ of the next MCS found cannot be dominated by $r$ (note that $\hat{r} = r$ is allowed) and the corresponding assignment must be different from the previous ones. Hence, all assignments associated to each MCS would be enumerated with this modified version. For $\epsilon = 0$, this results in enumerating the whole efficient set. For $\epsilon > 0$, at least one efficient solution would be enumerated per MCS but non-efficient solutions could be enumerated as well, which makes this modified version of lesser value for $\epsilon > 0$.

## 4.2. Enumerating all nondominated points through re-approximations

The value of $\epsilon$ has impact on the size of the encoding of the objective functions and on the search process, that depends on which of the two described approximation methods is used. In the interval-based approximation, $|D_k|$ represents the number of intervals into which the range of objective values $\{\ell_k, ..., u_k\}$ is split, for $k \in \{1, ..., p\}$. The larger $\epsilon$ is, the smaller is $|D_k|$, whereas $|D_k|$ tends to $u_k - \ell_k + 1$ as $\epsilon$ tends to zero. Hence, setting a variable $y_{k,d}$ to 1 allows to skip larger regions of the objective space when $\epsilon$ is large, than when it is small. In the coefficient-based approximation, the coefficients of the approximate function $f'_k$ tend to 1 and the upper limit $u'_k$ tends to $n$, as $\epsilon$ grows. Consequently, the size of the encoding of the objective functions, and therefore the computational burden thereof, is expected to decrease when increasing $\epsilon$.

This section describes two algorithms taking advantage of the above observations, one using coefficient-based approximations (Section 4.2.1) and another using interval-based approximations (Section 4.2.2). The main idea is to start with a large $\epsilon$, then to iteratively

**Algorithm 2:** CoRe algorithm: Enumerate all nondominated points for $f(x)$ subject to the problem constraints encoded in $\phi$, through re-approximations using the coefficient-based approximation method. The initial approximation ratio is $\epsilon$.

Input: $f, \phi, \epsilon$
Output: An optimal set $\mathcal{A}$
1 & A & {}
2 repeat
3 & $\phi^h \leftarrow \phi$
4 & $\mathcal{L} \leftarrow \{$
5 & for $(k=1)$ to $(p)$ do
6 & $f'_k \leftarrow \text{CoeffApprox}(f_k, 1+\epsilon)$ // Update approximation of $k$-th objective function
7 & $\phi^h \leftarrow \text{EncodeObjFunction}(f'_k)$
8 & $\phi^h \leftarrow \phi^h \cup \phi^h$
9 & $D_k \leftarrow \text{ComputeD}(f'_k, 1)$ // Compute $D_k$ for $f'_k$
10 & for $d \in D_k$ do
11 & $\left[ (y_{k,d}, \phi^h) \leftarrow \text{EncodeLT}(k, d) \quad \text{// Encode} \atop \text{y}_{k,d} = 1 \Rightarrow f'_k(x) < d \right.$
12 & $\left. \phi^h \leftarrow \phi^h \cup \phi^h \right)$
13 & for $(v, f(v)) \in \mathcal{A}$ do
14 & $\phi^h \leftarrow \phi^h \cup \{(y_{1,f'_k(v)}, \dots, y_{p,f'_k(v)})\}$ // Block region
15 & $\left. \mathcal{L} \leftarrow \mathcal{L} \cup \{f'_k(v)\} \quad \text{// weakly dominated by } f'_k(v) \right.$
16 & $(\mathcal{A}', \mathcal{L}') \leftarrow \text{MCS-Approx}(f, f', \phi^h, \mathcal{D})$
17 & $\mathcal{A} \leftarrow \text{Non-dominated}(\mathcal{A} \cup \mathcal{A}')$ // $\mathcal{A}$ stores a $(1+\epsilon)$-approximation set
18 & $\mathcal{L} \leftarrow \text{Non-dominated}(\mathcal{L} \cup \mathcal{L}')$ // $\mathcal{L}$ is a lower bound set
19 & $\epsilon \leftarrow \text{UpdateRatio}(\epsilon)$
20 until $f = f'$ or StopCriterion($\epsilon$)
21 return $\mathcal{A}, \mathcal{L}$

use Algorithm 1 to find a $(1+\epsilon)$-approximation for progressively smaller values of $\epsilon$, until $\epsilon$ is small enough. The approximation ratio of the approximation set is improved from iteration to iteration. If $\epsilon = 0$ in the last iteration, then an optimal set is returned.

### 4.2.1. Coefficient-based re-approximations

Algorithm 2 describes a procedure to sequentially update an approximation set $\mathcal{A}$ while improving the corresponding approximation ratio. Given the objective functions $f$, a set of clauses $\phi$ representing the problem constraints, and an initial value of $\epsilon$, the algorithm runs until $\mathcal{A}$ represents an optimal solution set. In each iteration, the algorithm performs four main steps: (1) it encodes the (approximate) objective functions in lines 5–12 given the current value of $\epsilon$; (2) a region in objective space within a $(1+\epsilon)$ ratio of $\mathcal{A}$ is blocked in lines 13–14; (3) in lines 16–18, $\mathcal{A}$ and $\mathcal{L}$ are updated using Algorithm 1 to search for solutions mapping to points outside the blocked region, and to make sure that $\mathcal{A}$ becomes a $(1+\epsilon)$-approximation and $\mathcal{L}$ the respective lower bound set; (4) the value $\epsilon$ is updated (it is decreased) at the end of each iteration, in line 19.

In the first step, the approximate version, $f'_k$, of the objective function $f_k$ is computed first for each $k \in \{1, ..., p\}$ in line 6, then the set of clauses $\phi^h$ encoding $f'_k$ (e.g., as explained in Section 2.2.2) is computed in line 7, and is added to $\phi^h$. The set $D_k$ is computed in line 9, where the second argument of the procedure ComputeD indicates the interval-approximation factor for computing $D_k$. In this case (when the coefficient-based approximation is used), this factor must be 1 and therefore, $D_k$ is complete. The procedure EncodeLT $(k, d)$ in line 11 computes the set of clauses $\phi^h$ that encodes the implication $y_{k,d} = 1 \Rightarrow f'_k(x) < d$. These clauses, $\phi^h$, are then added to $\phi^h$ in line 12.

The second step is skipped (only) in the first iteration of the repeat- until loop as $\mathcal{A}$ is empty, i.e., no solutions are known. Let $\epsilon'$ be the value of $\epsilon$ from the previous iteration, hence, $\epsilon' > \epsilon$. At the beginning of a new iteration, the set $\mathcal{A}$ is a $(1+\epsilon')$-approximation. Although not
<figure><img src="image_6.png" /><figcaption>Fig. 3. Example of the first two iterations of Algorithm 2: (a) illustration of the problem solved in the first iteration for $\epsilon = 3$, (b) the region blocked in the beginning of the second iteration for $\epsilon = 1$, and (c) the problem solved after.</figcaption></figure>

all the efficient solutions are, in general, expected to be within the updated, and tighter, $(1+\epsilon)$ ratio from the solutions stored in $\mathcal{A}$, some may be. To take advantage of this, for each solution $v$ stored in $\mathcal{A}$, the region weakly dominated by $f'(v)$ is blocked in line 14. The union of the blocked regions is represented by $\mathcal{L}$ in line 15, i.e., $\mathcal{L} = f'(\mathcal{A})$. At this point, $\phi^h$ encodes the problem constraints, the $p$ approximate objective functions $f'_1, \dots, f'_p$, the $y$ variables needed, and the objective-space region that is blocked. Hence, $\phi^h$ is satisfiable only if there is a solution that satisfies the problem constraints, $\phi$, and which maps to a point outside the blocked region.

In the third step, Algorithm 1 is executed in line 16 and returns the set $\mathcal{A}'$ that represents a $(1+\epsilon)$-approximation with respect to the set of feasible solutions mapping to points outside the blocked region, and the respective lower bound set, $\mathcal{L}'$. Hence, $\mathcal{A}' \cup \mathcal{A}$ is a $(1+\epsilon)$-approximation set, and $\mathcal{L}' \cup \mathcal{L}$ is a lower bound set. The set of solutions that are not dominated by any other solution in $\mathcal{A}' \cup \mathcal{A}$ is then stored in $\mathcal{A}$ in line 17. Similarly, only the mutually nondominated points in $\mathcal{L}' \cup \mathcal{L}$ are stored in $\mathcal{L}$ in line 18.

Finally, in the last step, the approximation ratio $\epsilon$ is updated (line 19). Different update strategies may be used, where the only requirement is to decrease the value of $\epsilon$. Setting $\epsilon$ to zero ensures that all points in the Pareto front are enumerated at the end of the next iteration. The algorithm will stop when all approximate functions are equal to the corresponding original objective functions ($f = f'$), which will happen if $\epsilon = 0$ but may also happen for small enough values of $\epsilon$ greater than zero. Note that the algorithm may be stopped earlier using a user-defined criterion such as when a time limit or a desired approximation ratio is reached.

Fig. 3 illustrates two iterations of Algorithm 2 for the problem defined in expression (3). Fig. 3(a) shows the feasible points under $f'$ given $\epsilon = 3$ (blue squares) where $f'_1(x) = x_1 + x_2 + x_3 + x_4 + 1$ and $f'_2(x) = 4\bar{x}_1 + 4\bar{x}_2 + 4\bar{x}_3 + 4\bar{x}_4 + 1$. Since they are all nondominated, Algorithm 1 finds one solution mapping to each one of these points (black circles), which are returned in $\mathcal{A}'$ in line 16. Note that, (7, 13) is weakly dominated by (7, 6) and thus, it is discarded in line 17. The resulting set $\mathcal{A} = \{(0, 0, 0, 0), (1, 22), ((0, 1, 0, 0), (4, 17)), ((1, 0, 1, 1), (7, 6)), ((1, 1, 1, 1), (10, 1))\}$ is a $(1+3)$-approximation set. Assuming that $\epsilon$ is then set to 1 in line 19, the next iteration encodes the new approximate objective functions, which are now $f'_1(x) = 2x_1 + 2x_2 + x_3 + 2x_4 + 1$ and $f'_2(x) = 4\bar{x}_1 + 4\bar{x}_2 + 4\bar{x}_3 + 4\bar{x}_4 + 1$. For each solution $v$ stored in $\mathcal{A}$, $f(v)$ and $f'(v)$ are represented by a diamond and square, respectively, in Fig. 3(b) and they are linked by an arrow. The shaded region in Fig. 3(b) is the region blocked in lines 13–14. It is lower bounded by the images of

the assignments in $\mathcal{A}$ according to the newly defined $f'$, namely points (1,17), (3,13), (6,5), and (8,1). Hence, all solutions $v$ such that $f'(v)$ is in the shaded region are now considered infeasible. Fig. 3(c) shows the remaining feasible points under $f'$ (squares) and the points under $f$ associated to them (circles), where the solutions mapping to black circles are returned in $\mathcal{A}'$ in line 16. At the end of the iteration, the set $A_\epsilon$ represented in $\mathcal{A}$, where $f(A_\epsilon) = \{(1, 22), (2, 17), (5, 13), (7, 6), (10, 1)\}$, is a $(1+1)$-approximation and the associated lower bound set is $\mathcal{L} = \{(1, 17), (2, 13), (4, 9), (6, 5), (8, 1)\}$.

Note that, in the example, the approximate objective function $f'_2$ does not change for the two settings of $\epsilon$. In this case, there is no need to repeat lines 7–12 for $k=2$. The algorithm can be easily adapted to avoid repeating the encoding of approximate functions that do not change.

### 4.2.2. Interval-based re-approximations

Algorithm 3 is similar to Algorithm 2 in the sense that it follows an identical scheme to iteratively improve an approximation set $\mathcal{A}$ and the corresponding approximation ratio, until $\mathcal{A}$ represents an optimal set. Similarly, the input of Algorithm 3 are a set of objective functions, $f$, the problem constraints encoded as a set of clauses $\phi$, and an initial value of $\epsilon$. The algorithm has a preprocessing step in lines 2–4 where all objective functions are encoded. Then, in each iteration, the algorithm also performs four main steps: (1) determine the objective functions intervals by computing $D_1, \dots, D_p$, and encode the associated $y$ variables, in lines 6–10; (2) block the region weakly dominated by $\mathcal{A}$, in lines 11–16; (3) search for a $(1+\epsilon)$-approximation set outside the blocked region, and update $\mathcal{A}$ in lines 18–19; (4) update (i.e., decrease the value of) $\epsilon$ in line 21.

Unlike Algorithm 2, Algorithm 3 directly encodes the original objective functions $f_1, \dots, f_p$ and does it just once at the beginning of the algorithm (line 3). The corresponding clauses are added to $\phi^h$ (line 4). In this case, only the sets $D_1, \dots, D_p$ will (possibly) change from iteration to iteration as $\epsilon$ decreases. Therefore, the first main step of the repeat-until loop is to update the sets $D_1, \dots, D_p$ according to the current value of $\epsilon$ (line 7) and the $y$ variables associated to the intervals they represent (line 9). Note that, the second argument of ComputeD is $1+\epsilon$, so that the values in $D_1, \dots, D_p$ represent the approximation intervals. The second main step is to block the region weakly dominated by each new solution $v$, stored in $\mathcal{A}'$ (lines 12–16). After line 17, $\mathcal{L} = f(\mathcal{A})$ represents the blocked region. Firstly, the procedure EncodeLT is called (line 14) for the value of each component of $f(v)$, stored in $d_k$, because $f_k(v)$ is not necessarily in $D_k$ and thus, $y_{k,d_k}$ may not have
**Algorithm 3:** IntRe algorithm: Enumerate all nondominated points of *f* subject to the problem constraints encoded in *φ*, through re-approximations using the interval-based approximation method. The initial approximation ratio is *ε*.

Input: $f, \phi, \epsilon$
Output: An optimal set $\mathcal{A}$
1 & $\mathcal{A} \leftarrow \mathcal{A}' \leftarrow \phi^h \leftarrow \{\}$
2 for ($k=1$) to ($p$) do
3 & $\phi' \leftarrow EncodeObjFunction(f_k)$ // Encode $k$-th
4 & $\phi^h \leftarrow \phi^h \cup \phi'$ // objective function
5 repeat
6 & for ($k=1$) to ($p$) do
7 & & $D_k \leftarrow ComputeD(f_k, 1+\epsilon)$ // Compute $D_k$ for $f_k$
8 & for $d \in D_k$ do
9 & & $\left[ (y_{k,d}, \phi'') \leftarrow EncodeLT(k, d) \right.$ // Encode
10 & & $\left. \quad \phi^h \leftarrow \phi^h \cup \phi'' \right]$ // $y_{k,d} = 1 \Rightarrow f_k(x) < d$
11 & for $(v, f(v)) \in \mathcal{A}'$ do
12 & & for ($k=1$) to ($p$) do
13 & & & $d_k \leftarrow f_k(v)$
14 & & & $(y_{k,d_k}, \phi'') \leftarrow EncodeLT(k, d_k)$
15 & & $\phi^h \leftarrow \phi^h \cup \phi''$
16 & & $\phi^h \leftarrow \phi^h \cup \{(y_{1,d_1} \vee \dots \vee y_{p,d_p})\}$ // Block region weakly
17 & & $\phi^h \leftarrow \phi^h \cup \{(y_{1,d_1} \vee \dots \vee y_{p,d_p})\}$ // dominated by $f(v)$
18 & $\mathcal{L} \leftarrow \{f(v) | (v, f(v)) \in \mathcal{A}\}$
19 & $(\mathcal{A}', \mathcal{L}') \leftarrow \text{MCS-Approx}(f, f, \phi^h, \mathcal{D})$
20 & $\mathcal{A} \leftarrow \text{Non-dominated}(\mathcal{A} \cup \mathcal{A}')$ // $\mathcal{A}$ stores a
21 & $\epsilon \leftarrow \text{UpdateRatio}(\epsilon)$ // ($1+\epsilon$)-approximation set
22 until $\mathcal{A}' = \{\}$ or $\text{StopCriterion}(\epsilon)$
23 return $\mathcal{A}, \mathcal{L}$

been encoded in the previous for loop (in line 9). Then, the clause $(y_{1,d_1} \lor \dots \lor y_{p,d_p})$ is permanently added to $\phi^h$ (line 16) to block the region weakly dominated by $f(v)$ in the remaining iterations.

Similarly to Algorithm 2, Algorithm 1 is executed in line 18 to find a $(1+\epsilon)$-approximation with respect to the solutions outside the blocked region, $\mathcal{A}'$, and the associated lower bound set, $\mathcal{L}'$. Then $\mathcal{A}$ and $\mathcal{L}$ are updated with $\mathcal{A}'$ and $\mathcal{L}'$, respectively, in lines 19–20. Hence, $\mathcal{A}$ becomes a $(1+\epsilon)$-approximation and $\mathcal{L}$ the associated lower bound set. The value of $\epsilon$ is then updated in line 21. If Algorithm 1 returns an empty set as $\mathcal{A}'$ (in line 18), then $\mathcal{A} = \mathcal{L}$ is the Pareto front, and the algorithm terminates. As in Algorithm 2, Algorithm 3 may be interrupted when a user-defined criterion is met.

Fig. 4 illustrates two iterations of Algorithm 3 for the problem in expression (3). With $\epsilon$ initially set to 3, $D_1 = \{1, 4, 16\}$ and $D_2 = \{1, 4, 16, 64\}$ in the first iteration (line 7). Fig. 4(a) shows all feasible points (circles), the points returned in $\mathcal{A}'$ in line 18 (black circles), i.e., $\{(3, 15), (10, 1)\}$, and the points in the respective lower bound set $\mathcal{L}' = \{(1, 4), (4, 1)\}$ (orange crosses). At the end of the first iteration, the value of $\epsilon$ is updated to 1 (line 21). In the second iteration, after computing $D_1 = \{1, 2, 4, 8, 16\}$ and $D_2 = \{1, 2, 4, 8, 16, 32\}$, the regions weakly dominated by points found in the previous iteration are blocked in lines 11–16. Such regions are shaded in Fig. 4(b)) where the dashed lines represent the values $d_1, \dots, d_p$ computed in line 13 and used in line 14 to encode the $y$ variables needed to construct the blocking clauses. Fig. 4(c) shows the remaining feasible points outside the blocked region, the point set $\{(1, 22), (7, 6)\}$ returned in $\mathcal{A}'$ (black circles) and the respective lower bound set $\mathcal{L}' = \{(1, 16), (4, 4)\}$ (orange crosses) computed in line 18. At the end of this iteration, $\mathcal{A} = \{(3, 15), (10, 1), (1, 22), (7, 6)\}$ is a $(1+1)$-approximation set, and $\mathcal{L} = \{(3, 15), (10, 1), (1, 16), (4, 4)\}$ is the associated lower bound set. Note

that, since the region weakly dominated by $\{(3, 15), (10, 1)\}$ is already blocked, in the third iteration new clauses will be created to block only the regions weakly dominated by $(1, 22)$ and $(7, 6)$.

### 4.3. Remarks

Algorithm 1 finds a $(1+\epsilon)$-approximation set, provided that (some of) the values of the objective functions were previously encoded using some unary representation. The (re-approximation) Algorithms 2 and 3, named CoRe and IntRe, respectively, are responsible for such encoding, and iteratively call Algorithm 1. Both are flexible in the sense that they can start searching for a $(1+\epsilon^s)$-approximation set for any $\epsilon^s \ge 0$, then iteratively refine $\epsilon^s$, and stop when a $(1+\epsilon^d)$-approximation set is found, given a desired value of $\epsilon^d$. If $\epsilon^d = \epsilon^s$, then the algorithm stops after just one call to Algorithm 1. CoRe and IntRe may also be terminated when a given time budget is reached. In any case, they return in $\mathcal{A}$ the best solutions found so far, and as long as the first iteration terminates within the time budget, the algorithms provide a lower bound set, $\mathcal{L}$. This delimits the location of the Pareto front, and provides an approximation factor of $I_\epsilon(\mathcal{A}, \mathcal{L})$ which is possibly tighter than $(1+\epsilon^c)$, where $\epsilon^c$ is the value of $\epsilon$ in the last completed iteration. Hence, these algorithms are useful to a Decision Maker interested in learning about the location of the Pareto front, or in finding the best approximation possible, as soon as possible or, within a time limit.

The anytime performance of the algorithms, i.e., the quality of the obtained approximation set at any time during their execution, depends on multiple factors such as the problem instance, the approximation method, the starting $\epsilon$ value and how it is updated, how many times Algorithm 1 and the SAT solver are called, and so on. For example, a small initial value for $\epsilon$ may lead to a good approximation set at the end of the first iteration, but finding it may not be possible within the time limit. On the other hand, a larger initial value for $\epsilon$ will possibly lead to a worse first approximation, but which is found faster, allowing the algorithm to progressively find better approximation sets.

The runtime of Algorithm 1 is also influenced by the size of the approximation set found. Papadimitriou and Yannakakis (2000) showed that, for a fixed number of objectives, $p$, there is a $(1+\epsilon)$-approximation set of polynomial size in $n$ and $1/\epsilon$ for any $\epsilon > 0$. This is the case even with different $\epsilon_k$ values for each objective function $k \in \{1, \dots, p\}$ and with the $i$th one set to $\epsilon_i = 0$ (Herzel et al., 2021a). The proposed algorithms could be easily adapted to consider this case, even with different update strategies for each $\epsilon_k$. This indicates that the number of MCSs enumerated by Algorithm 1 may be polynomial, even though $Y_N$ may have exponential size. In fact, the interval-based approximation ensures such polynomial-size approximation sets because it relies on the same hypergrid partitioning of the objective space as in the proofs of the aforementioned articles. Hence, with an initial $\epsilon > 0$, IntRe (and possibly CoRe for some instances) can first search for an approximation set of polynomial size, giving an idea of the location and extent of the Pareto front, and then use the time left to enumerate as many nondominated points as possible (by updating $\epsilon$ to zero). The lower bound set could also be updated along the run to tighten the known approximation factor. In contrast, with an initial $\epsilon = 0$, the algorithms may only be able to return a (possibly not well distributed) subset of the Pareto front with no known approximation factor when $|Y_N|$ is large (e.g., exponential in $n$).

There are a few differences distinguishing the re-approximation algorithms that stem from the approximation method that each one uses. With the interval-based approximation, the set of clauses that encode the objective functions is the same independently of the value of $\epsilon$. One advantage of this is that the objective functions are encoded just once in IntRe and afterwards, the $y_{k,d}$ variables (i.e., the unary representation of the objective values) are encoded as needed, depending on $\epsilon$. However, this implies that a large encoding may be needed even if $\epsilon$ is large. Another advantage of IntRe is that as new
<figure><img src="image_7.png" /><figcaption>Fig. 4. Example of the first two iterations of Algorithm 3: (a) illustration of the problem solved in the first iteration for $\epsilon = 3$, (b) the region blocked in the beginning of the second iteration for $\epsilon = 1$, and (c) the problem solved after.</figcaption></figure>

solutions are found, the region weakly dominated by those points can be permanently blocked.

A disadvantage of the coefficient-based approximation is that, since the approximate objective functions may change from iteration to iteration, they have to be re-encoded at the beginning of each iteration, and new clauses have to be added to block the region dominated by the points in $\mathcal{A}$ using the new $y_{k,d}$ variables associated with the newest encodings. Therefore, as $\epsilon$ tends to zero, the total amount of time spent in these re-encoding and blocking steps may be more time consuming in CoRe than in IntRe. An advantage of the coefficient-based method is that the encoding of the approximate objective functions is potentially much smaller than for $\epsilon = 0$ (and than in IntRe). The size of the encoding depends on three factors: the method used to obtain the unary encoding, the coefficients of the function being encoded, and the function upper bound. Although the influence of the coefficients may not be easy to quantify, the smaller they are the smaller is the function upper bound, and the more equal coefficients in the function the smaller is the number of different values the function can take (e.g., if they are all equal, then the function can only take up to *n* different values). Hence, in general, the larger $\epsilon$ is, the smaller the encoding is expected to be.

Finally, note that, for some problem instances, running CoRe for any setting of $\epsilon$ is equivalent to solving the exact case, $\epsilon = 0$, and therefore, the approximation set returned could be of exponential size. For example, when all the coefficients of the objective functions are equal, or are powers of $(1+\epsilon)$. In such cases, IntRe may be preferred. It is also important to note that if $c^{max}$ is the largest coefficient in the objective functions, then running CoRe with $(1+\epsilon) > c^{max}$ is equivalent to running it with $(1+\epsilon) = c^{max}$ because in either case, the non-zero coefficients of the approximate functions will all be one. Additionally, if the *p* approximate functions are all equal, then the corresponding approximation sets will have size one.

## 5. Computational experiments

This section evaluates the performance of the algorithms proposed in Section 4. In particular, we show some preliminary results on the anytime performance of IntRe and CoRe (see Section 5.2) to highlight their advantages and potential, and then we compare them with state-of-the-art solvers (see Section 5.3).

### 5.1. Experimental setup

All experimental results were obtained on a server with processor Intel(R) Xeon(R) CPU E5-2630 v2 @ 2.60 GHz with 64 GB of memory. The IntRe and CoRe algorithms were implemented in C++ by extending Open-WBO (Martins et al., 2014) and adapting it for the multi-objective case. Glucose SAT solver (version 4.1) was used in MCS-Approx and the CLD algorithm (Marques-Silva et al., 2013) was used to compute MCSs. In the implementations we have some preprocessing steps to reduce the size of the encodings of the objective functions. For example, if, at some point during the execution, the SAT solver infers the value of some variable, i.e., that the variable holds that value for any feasible solution, then this variable can be excluded when (re-)encoding the objective functions. Another example is the computation of an upper bound on the feasible values of each objective function, tighter than $u_k$, and encoding the function only up to that value. Another relevant implementation aspect is the computation of a complete domain $D_k$. A straightforward way is to simply consider the full range {$l_k, ..., u^+$}, but this set may be very large with respect to $f(X)$. To avoid such a computation, we first use a dynamic programming approach to infer all values that $f_k$ may take. If, at some point, such a computation becomes too time consuming because there are too many different values, then we simply compute the full range (in linear time in the range size).

The algorithms IntRe and CoRe were tested with different input parameters and different update strategies of the approximation ratio. In the plots shown in this section, the name of the algorithms is followed by a tuple indicating such settings, which may include one or two components. The first component corresponds to the value of $(1+\epsilon)$, and the second, if present, indicates the value by which $\epsilon$ is divided in each iteration. For example, IntRe(2, 10) corresponds to running IntRe with a starting value of $\epsilon = 1$ (i.e., $2 = 1+\epsilon$), and in the second iteration $\epsilon$ will be set to $1/10$, in the third it will be set to $0.1/10$, and so on. Another example is CoRe(4) which corresponds to running CoRe with $\epsilon = 3$, and terminating the algorithm after the first iteration. Since IntRe and CoRe are equivalent when $\epsilon$ is initially set to 0, which corresponds to encoding the objective functions in an exact way and call Algorithm 1 just once to enumerate all nondominated points, this particular case will be referred as MCS-Approx(1).

The following MOBO problems were considered: the Multiobjective Set Covering Problem (MSCP) (Bergman and Ciré, 2016; Soh et al., 2017), and the Multiobjective Development Assurance Problem (MDAP) (Bieber et al., 2011). MSCP is a generalization of the classical
<figure><img src="image_8.png" /><figcaption>Fig. 5. Performance of IntRe and CoRe algorithms along the run, at time $t \in [0, 3600]$ seconds, on a given MCSP instance: (a) the hypervolume indicator of the set of solutions known; (b) the value of the current setting of $(1 + \epsilon)$.</figcaption></figure>

set covering problem that consists of deciding which subsets of a ground set $X$, and pre-specified in a set $A$, to select such that the total cost associated to selecting the subsets is minimized given that all elements of $X$ are covered (i.e., each element is in at least one of the selected subsets). MDAP encodes different levels of rigor of the development of a software or hardware component of an aircraft. The Development Assurance Level (DAL) 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 in order to decrease the development costs.<fnref n="1" />

## 5.2. Anytime performance example

To show the potential of IntRe and CoRe, we show some results on the anytime performance of these algorithms on a 2-objective instance of MSCP, named “2scp43A”,<fnref n="2" /> for which $n = 200$ and $m = 40$, and where $l_1 = 0$, $u_1 = 21075$, $l_2 = 0$, and $u_2 = 21366$. We compare, against the exact version (i.e., $\epsilon = 0$), the performance of IntRe and CoRe with different initial value of $\epsilon \in \{1, 100\}$, which is divided by 10 at each iteration until $\epsilon < 1e^{-4}$, in which case $\epsilon$ is set to 0. The algorithm runs until it finishes enumerating the Pareto front, or the algorithm achieves the time limit of 3600 seconds or the memory limit of 8 GB. For Fig. 5, we recorded the following information: each new feasible solution found and the time when it was found, the time when a new iteration of IntRe/CoRe started, and to which value $\epsilon$ was set. Fig. 5(a) shows the quality of the set of solutions found along the run according to the hypervolume indicator (Zitzler and Thiele, 1998). The indicator was computed with the reference point set to $r = (1, 1)$ and after normalizing the objective values of each point by dividing them by the maximum of the respective objective function plus one, i.e., by 21076 and 21367, in the case of the first and the second objectives, respectively. Fig. 5(b) shows the value of $(1 + \epsilon)$ along the run.

Fig. 5(a) shows that, in this instance, the algorithms that first search for an approximation of the Pareto front can find better sets of solutions sooner than the algorithm enumerating all nondominated points from the start. For example, the quality of the set of solutions found by IntRe(101, 10) and CoRe(2, 10) were always better than MCS-Approx(1) for any time $t \in [0, 3600]$. Fig. 5(b) shows that, for IntRe, the value

of $\epsilon$ rapidly decreased to 1, and therefore, somewhat justifies the high hypervolume values of the sets of solutions found by IntRe early in the run.

Fig. 6 shows in more detail the sets found by each call of Algorithm 1 within each algorithm tested. The symbols identify the iteration and the setting of $1 + \epsilon$ for that iteration. The opaque symbols represent the points found in that iteration, and the faded symbols represent the lower bound set at the end of the iteration. For example, in Fig. 6(a), IntRe(2, 10) found the four blue circle points in the first iteration for which $1 + \epsilon$ was 2, it found the green triangles filled white points in the second iteration where $1 + \epsilon = 1.1$, and found the yellow-star points in the third iteration for which $1 + \epsilon = 1.01$. The three faded blue circles represent the lower bound set at the end of the first iteration, and the faded green triangles represent the lower bound set at the end of the second iteration. Note that the $(1+\epsilon)$-approximation set at the end of the $i$th iteration is the set of points found in the first $i$ iterations (excluding the dominated ones). Hence, at the end of the second iteration of IntRe(2, 10), the set $A$ contains the points represented by triangles, plus the points represented by circles and which are not dominated by the triangles.

Fig. 6(c) illustrates the nondominated points found by MCS-Approx(1). Note that the algorithm did not finish, and therefore, the Pareto front may be incomplete. The other plots show, from iteration to iteration, how both IntRe and CoRe narrow down the region where the Pareto front is located, and show that, even by setting $\epsilon$ to large values, allows to exclude a large portion of the objective space. Note that, in this instance, CoRe(101, 10) found points closer to the Pareto front than IntRe(101, 10) did when considering the same iteration. However, since more points were needed to guarantee the approximation factor by CoRe, it took more time to achieve the same hypervolume values that IntRe achieved.

These results illustrate the potential of IntRe and CoRe. In particular, they show that both IntRe and CoRe can outperform algorithms that focus on finding points from the Pareto frontier from the start, with respect to the quality of the set, and at any point in time. However, the observed behavior might not generalize for other types of problems or problem instances. This is discussed in more detail in the next section.

## 5.3. Algorithm comparison

The proposed algorithms are evaluated next on multiple instances of MSCP (see Section 5.3.1) and of MDAP (see Section 5.3.2). The size of the formulas needed in each of the proposed algorithms with different settings is compared. Moreover, the quality of the approximation sets obtained within a time budget is evaluated and compared against the quality of the approximation sets obtained with state-of-the art algorithms.

<footnote n="1">¹ The benchmark instances are available online at https://www.lifl.fr/LION9/challenge.php. Although they define a lexicographic order to the objective functions, in the context of this paper, we ignore this order and compute the Pareto frontier.</footnote>

<footnote n="2">² Instance from https://github.com/vOptSolver/vOptLib/blob/master/SCP/scp1998.md</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Illustration of the points found by each algorithm at each iteration. The name of the algorithm is in the title of the plot, and the legend indicates the symbol used to represent the points found in each iteration, and the setting of $(1 + \epsilon)$ on that iteration.</figcaption></figure>

### 5.3.1. Multiobjective set covering problem

For MSCP, we used the instances<fnref n="3" /> by Bergman and Ciré (2016). An instance was randomly selected for each triple $(n, m, p)$ of $n \in \{100, 150, 200\}$ variables, $m \in \{n-20, n-30, ..., n-90\}$ constraints, and $p \in \{3, 4\}$ objectives for a total of 48 instances. In these instances, each constraint involves 5 non-zero coefficients, and each coefficient in the first objective function is 1, and is an integer between 1 and 100 in the remaining objective functions. Additionally, we randomly generated a (new) set of 30 highly constrained instances. Such instances are similar to the previous ones, with the same characteristics and parameters, but the number of constraints is $m \in \{n, 2n, 5n, 10n, 20n\}$. The IntRe and CoRe algorithms were tested in these instances with different settings (different starting value of $\epsilon \in \{0.1, 1, 10, 100\}$) and the most relevant to the discussion are shown in Figs. 7 and 8. These algorithms are compared against state-of-the-art ones: an algorithm based on binary decision diagrams (Bergman and Ciré, 2016), and an algorithm based on the $\epsilon$-constraint method (Tamby and Vanderpooten, 2021), which are here referred to as BDD and TamVan, respectively.

The results in the figures evaluate the algorithms outputs and total encoding sizes considering a time limit of 10 min per instance and a memory limit of 8 GB. The figures show the empirical cumulative distributive function over the number of instances for: the hypervolume indicator of the output sets (Figs. 7) the total number of clauses for the encoding of the objective functions (Figs. 8(c) and 8(f)); the a priori approximation ratio given by the $\epsilon$-indicator value for the output sets considering the respective lower bound set, if known (Figs. 8(a) and

8(d)); and a posteriori approximation ratio given by the smallest $\epsilon$-indicator value obtained considering the Pareto front (if known) and the lower bound sets returned by all algorithms (Figs. 8(b) and 8(e)). Hence, for each algorithm, the value $t$ for the second axis at value $i \in \{1, ..., 48\}$ of the first axis, indicates that the evaluation for $i$ instances was less or equal than $t$, and was greater than $t$ for $48 - i$ instances. Therefore, the closer the line is to the top and to the left of the plot, the better.

Fig. 7 shows the hypervolume indicator after normalizing the objective values of the points in the output sets, and computed considering the reference point $r = (1, ..., 1)$. The normalization was performed by considering the largest values for each coordinate among all output sets from all algorithms and multiplying them by 1.1. This ensures that, when considering $r = (1, ..., 1)$ as the reference point, all points in the output sets contribute to the hypervolume indicator. Regarding the low-constrained data set, Fig. 7(a) shows that CoRe(2) and CoRe(11) already perform better than MCS-Approx(1), i.e., when CoRe searches for an approximation set and stops after the first iteration, the returned approximation sets (for $\epsilon = 1$ and $\epsilon = 10$) have greater hypervolume value than the set of nondominated solutions found by MCS-Approx(1) in these instances. The reason for this seems to be related to the fact that the true approximation ratio (i.e., the $\epsilon$-indicator value with the Pareto front as the reference set) of the returned approximation sets seems to be very close to 1 even though the initial $\epsilon$ value may be high. See in Figs. 8(a) and Fig. 8(b), that for approximately 20 instances, the a posteriori approximation ratio is much smaller than the a priori ratio and that the y-axis values of most of the lines in Fig. 8(a) are now very close to 1 in Fig. 8(b). Since the approximation with a priori guarantee was found within 52 s in average for these 20 instances, then by allowing the CoRe algorithm to continue to run for tighter approximation ratios

<footnote n="3">³ Instances at http://www.andrew.cmu.edu/user/vanhoeve/mdd/code/multiobjectiveCP2016.tar.gz</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Hypervolume indicator of the output sets obtained with: (a) CoRe algorithms; and (b) IntRe algorithms on low-constrained MSCP instances; and (c) CoRe and IntRe on high-constrained MSCP instances.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Results on the approximation factor (first two columns) of the output sets, and the number of clauses needed to encode the objective functions (last column) by the CoRe algorithms (first row) and by the IntRe algorithms (last row) on the low-constrained MSCP instances.</figcaption></figure>

the quality of the returned sets may be improved, and this is visible in the case of CoRe(11, 10). The hypervolume values achieved with CoRe indicate that, even though the points in the returned sets may not belong to the Pareto front, such sets dominate regions of the objective space that are not dominated by MCS-Approx(1).

The IntRe did not perform as well as CoRe in these instances (see Figs. 7(a) and 7(b)), even though it was able to warrant an approximation ratio for more instances (36), and faster (10 seconds in average) than IntRe (see Fig. 8(d)), but these ratios were not as good as those warranted by CoRe (see Figs. 8(a) and 8(b)). One reason for this observation may be the fact that CoRe is considering the exact representation of the first objective function (because all coefficients are 1) and therefore, the approximation ratio for the first objective is

one. Another advantage of CoRe over IntRe in these instances was the much smaller number of clauses needed by the former to encode the objective functions, and even though one of them is encoded in an exact way, as can be observed in Figs. 8(c) and 8(f).

Although BDD solved, within the time limit, about half of the low-constrained instances (see Fig. 8(a)) including all instances with $n = 100$, for the other half and for all of the high-constrained instances (see Fig. 7(c)) it was not able to finish within the time and memory limits. For this reason, in the overall, CoRe (11,10) performed much better than BDD. The TamVan algorithm outperformed all other algorithms in the low-constrained instances. However, for 22 out of the 30 high-constrained instances, TamVan did not find any Pareto-optimal
<figure><img src="image_12.png" /><figcaption>Fig. 9. Results on the: (a) hypervolume indicator; (b)–(c) approximation ratio; and (d) number of clauses of IntRe algorithms on the MDAP instances.</figcaption></figure>

solution. For such instances CoRe(11) shows the best performance. However, guaranteeing the *a priori* approximation ratio is more challenging and requires more than 10 min for such instances. This is the reason why no results are reported for CoRe(11, 10). The IntRe(11, 10) is able to guarantee the *a priori* ratio for 5 instances, but neither IntRe(11) nor IntRe(11, 10) outperforms CoRe(11).

### 5.3.2. Development assurance problem

The algorithms were ran on 48 instances<fnref n="4" /> of MDAP, which contain 7000 to 18000 variables, 20000 to 70000 constraints, and 7 objective functions. Since all coefficients of all objective functions are either -1 or 1 in these instances, only IntRe was tested and compared against a state-of-the-art logic-based algorithm (Terra-Neves et al., 2017), referred here as ParetoMCS.

Fig. 9 shows plots analogous to those shown in Section 5.3.1. It shows IntRe with two initial settings for $\epsilon$: 0.1 and 100. Though 100 is a very large value for $\epsilon$, this allowed IntRe to find, for 47 instances and in 108 s in average, an approximation set within this approximation ratio (and sometimes providing a much tighter one). This was much more difficult to achieve for an initial setting of 1.1 (see IntRe (1.1) and IntRe (101) in Figs. 9(b) and 9(c)). Consequently, IntRe (101,10) was able to really improve upon the approximation sets found by IntRe (101) and was able to find $(1+\epsilon)$-approximation sets for 18 instances, and to outperform MCS-Approx, ParetoMCS, and TamVan (see Figs. 9(a) to 9(c)). The only downside of IntRe is that it needed almost the same number of clauses as MCS-Approx (1) to encode the objective functions. Regarding TamVan, note that it does not search for nondominated points in any particular order, even though it could to some extent which would improve its performance. However, there would still be 18 instances for which it would not find any nondominated point within the time limit.

It is not always obvious whether IntRe or CoRe will perform better than the other for a given problem instance, nor with which parameter settings, though for very large objective functions, i.e., many variables and large coefficients, CoRe may be preferable due to the smaller number of clauses required. However, the results shown in this section indicate that they can be competitive against the state-of-the-art, particularly for highly-constrained instances, while warranting that the returned approximation set is within some approximation ratio, provided that they complete, at least, one iteration within the time limit. Moreover, depending on how the value of the desired approximation ratio is updated within IntRe and CoRe, the algorithm is able to tighten the *a priori* approximation ratio along the run. Additionally, the results also show that it can be helpful to first search for an approximation set before searching for the whole Pareto front. This allows to ignore a larger portion of the objective space sooner, and to find better approximation sets, with respect to the hypervolume indicator, than those found by the algorithms that search for nondominated points from the beginning of its execution.

## 6. Conclusions

This paper shows that with a unary encoding of the objective function values, there is a one-to-one correspondence between the nondominated points and Minimal Correction Subsets (MCSs). This implies that, with such encoding, any MCS enumeration algorithm can be used to enumerate all the nondominated points. This paper also proposes two approximation versions of such encoding which ensure that MCS enumeration algorithms provide a $(1+\epsilon)$-approximation set of the Pareto front, where $\epsilon \ge 0$. Additionally, two algorithms based on such approximations and on the idea of sequentially tightening $\epsilon$ until the set of solutions found has reached a desired approximation factor, or corresponds to the Pareto front, were proposed. Even if interrupted sooner, e.g., due to reaching a pre-defined time limit, the algorithms return the set of solutions found, and if at least the first iteration is completed, then they also provide a lower bound set and, consequently, an approximation factor tighter than the parameterized one. In fact, the proposed approximation based simply on modifying the objective function coefficients can be used with any MOBO solver to obtain a $(1+\epsilon)$-approximation set.

The preliminary results presented here show that such algorithms based on re-approximations can provide better anytime performance than the version dedicated to searching only nondominated points ($\epsilon = 0$) while requiring smaller encodings. These algorithms find particular relevance for instances for which the Pareto front is very large, or for which enumeration algorithms are not able to find the full Pareto front within a reasonable amount of time. These results show the potential of using such approximation versions of the encoding and the idea of iteratively finding better and better approximation sets. It remains unclear for which type of instances each approximation is the best for and with which parameter settings. This analysis is left for future work.

### CRediT authorship contribution statement

A. P. Guerreiro: Conceptualization, Methodology, Software, Validation, Investigation, Writing – original draft, Visualization. J. Cortes: Software, Writing – original draft. D. Vanderpooten: Conceptualization, Methodology, Writing – original draft. C. Bazgan: Conceptualization, Writing – review & editing. I. Lynce: Writing – review & editing, Supervision, Funding acquisition. V. Manquinho: Conceptualization, Methodology, Writing – original draft, Supervision, Project administration, Funding acquisition. J.R. Figueira: Conceptualization, Writing – review & editing, Supervision, Project administration, Funding acquisition.

### Data availability

Data will be made available on request.

<footnote n="4">⁴ https://www.cristal.univ-lille.fr/LION9/sampleB.gz</footnote>
## Acknowledgments

We would like to thank Satya Tamby for providing us the source code of TamVan algorithm, and promptly adapting it to our needs. This work was supported by Portuguese national funds through FCT - Fundação para a Ciência e Tecnologia, I.P., Portugal under projects UIDB/50021/2020, PTDC/CCI-COM/2156/2021, PTDC/CCI-COM/ 31198/2017 and DSAIPA/AI/0044/2018. The work by José Rui Figueira was also supported by Portuguese national funds through the FCT, Portugal under the project UIDB/00097/2020.

## References

Abío, I., Nieuwenhuis, R., Oliveras, A., Rodríguez-Carbonell, E., Mayer-Eichberger, V., 2012. A new look at BDDs for Pseudo-Boolean constraints. J. Artif. Intell. Res. 45, 443-480.

Asín, R., Nieuwenhuis, R., 2014. Curriculum-based course timetabling with SAT and MaxSAT. Ann. Oper. Res. 218 (1), 71-91.

Asín, R., Nieuwenhuis, R., Oliveras, A., Rodríguez-Carbonell, E., 2011. Cardinality networks: A theoretical and empirical study. Constraints 16 (2), 195-221.

Bacchus, F., Järvisalo, M., Martins, R., 2019. MaxSAT evaluation 2018: New developments and detailed results. J. Satisf. Boolean Model. Comput. 11 (1), 99-131.

Bailleux, O., Boufkhad, Y., 2003. Efficient CNF encoding of Boolean cardinality constraints. In: Principles and Practice of Constraint Programming. In: LNCS, vol.2833, Springer, pp. 108-122.

Bailleux, O., Boufkhad, Y., Roussel, O., 2006. A translation of Pseudo Boolean constraints to SAT. J. Satisf. Boolean Model. Comput. 2 (1-4), 191-200.

Bailleux, O., Boufkhad, Y., Roussel, O., 2009. New encodings of pseudo-boolean constraints into CNF. In: International Conference on Theory and Applications of Satisfiability Testing. In: LNCS, vol.5584, Springer, pp. 181-194.

Bazgan, C., Hugot, H., Vanderpooten, D., 2009. Implementing an efficient fptas for the 0-1 multi-objective knapsack problem. European J. Oper. Res. 198 (1), 47-56.

Bergman, D., Ciré, A.A., 2016. Multiobjective optimization by decision diagrams. In: Rueher, M. (Ed.), Principles and Practice of Constraint Programming - 22nd International Conference. CP 2016, Toulouse, France, September 5-9, 2016, Proceedings, In: Lecture Notes in Computer Science, vol.9892, Springer, pp. 86-95.

Bieber, P., Delmas, R., Seguin, C., 2011. Dalculus - theory and tool for development assurance level allocation. In: Flammini, F., Bologna, S., Vittorini, V. (Eds.), Computer Safety, Reliability, and Security - 30th International Conference. SAFECOMP 2011, Naples, Italy, September 19-22, 2011. Proceedings, In: Lecture Notes in Computer Science, vol.6894, Springer, pp. 43-56.

Bjørner, N., Phan, A., Fleckenstein, L., 2015. v Z - an optimizing SMT solver. In: International Conference on Tools and Algorithms for the Construction and Analysis of Systems. pp. 194-199.

Eén, N., Sörensson, N., 2006. Translating Pseudo-Boolean constraints into SAT. J. Satisf. Boolean Model. Comput. 2, 1-26.

Ehrgott, M., 2005. Multicriteria Optimization, second Ed. Springer.

Feng, Y., Bastani, O., Martins, R., Dillig, I., Anand, S., 2017. Automated synthesis of semantic malware signatures using maximum satisfiability. In: 24th Annual Network and Distributed System Security Symposium. NDSS 2017, San Diego, California, USA, February 26 - March 1, 2017, The Internet Society.

Grégoire, É., Izza, Y., Lagniez, J., 2018. Boosting MCSes enumeration. In: Lang, J. (Ed.), Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence. IJCAI 2018, July 13-19, 2018, Stockholm, Sweden, ijcai.org, pp. 1309-1315.

Herzel, A., Bazgan, C., Ruzika, S., Thielen, C., Vanderpooten, D., 2021a. One-exact approximate Pareto sets. J. Global Optim. 80 (1), 87-115.

Herzel, A., Ruzika, S., Thielen, C., 2021b. Approximation methods for multiobjective optimization problems: A survey. INFORMS J. Comput. 33 (4), 1284-1299.

Hölldobler, S., Manthey, N., Steinke, P., 2012. A compact encoding of Pseudo-Boolean constraints into SAT. In: KI 2012: Advances in Artificial Intelligence. In: LNCS, vol.7526, Springer, pp. 107-118.

Jackson, D., Estler, H., Rayside, D., et al., 2009. The Guided Improvement Algorithm for Exact, General-Purpose, Many-Objective Combinatorial Optimization. Technical Report MIT-CSAIL-TR-2009-033, MIT.

Jose, M., Majumdar, R., 2011. Cause clue clauses: Error localization using maximum satisfiability. In: Programming Language Design and Implementation. ACM, pp. 437-446.

Joshi, S., Martins, R., Manquinho, V.M., 2015. Generalized totalizer encoding for pseudo-boolean constraints. In: Pesant, G. (Ed.), Principles and Practice of Constraint Programming - 21st International Conference. CP 2015, Cork, Ireland, August 31 - September 4, 2015, Proceedings, In: Lecture Notes in Computer Science, vol.9255, Springer, pp. 200-209.

Karpinski, M., Piotrów, M., 2019. Encoding cardinality constraints using multiway merge selection networks. Constraints 24 (3-4), 234-251.

Karpinski, M., Piotrów, M., 2020. Incremental encoding of Pseudo-Boolean goal functions based on comparator networks. In: Pulina, L., Seidl, M. (Eds.), Theory and Applications of Satisfiability Testing - SAT 2020 - 23rd International Conference. Alghero, Italy, July 3-10, 2020, Proceedings, In: Lecture Notes in Computer Science, vol.12178, Springer, pp. 519-535.

Marques-Silva, J., Heras, F., Janota, M., Previti, A., Belov, A., 2013. On computing minimal correction subsets. In: International Joint Conference on Artificial Intelligence. pp. 615-622.

Martins, R., Manquinho, V., Lync, I., 2014. Open-WBO: A modular MaxSAT solver. In: International Conference on Theory and Applications of Satisfiability Testing. In: LNCS, vol.8561, Springer, pp. 438-445.

Ogawa, T., Liu, Y., Hasegawa, R., Koshimura, M., Fujita, H., 2013. Modulo based CNF encoding of cardinality constraints and its application to MaxSAT solvers. In: International Conference on Tools with Artificial Intelligence. IEEE, pp. 9-17.

Papadimitriou, C., Yannakakis, M., 2000. On the approximability of trade-offs and optimal access of web sources. In: Proceedings 41st Annual Symposium on Foundations of Computer Science. pp. 86-92.

Previti, A., Mencía, C., Järvisalo, M., Marques-Silva, J., 2017. Improving MCS enumeration via caching. In: Gaspers, S., Walsh, T. (Eds.), Theory and Applications of Satisfiability Testing - SAT 2017 - 20th International Conference. Melbourne, VIC, Australia, August 28 - September 1, 2017, Proceedings, In: Lecture Notes in Computer Science, vol.10491, Springer, pp. 184-194.

Roussel, O., Manquinho, V.M., 2009. Pseudo-Boolean and cardinality constraints. In: Biere, A., Heule, M., van Maaren, H., Walsh, T. (Eds.), Handbook of Satisfiability. In: Frontiers in Artificial Intelligence and Applications, vol. 185, IOS Press, pp. 695-733.

Safarpour, S., Mangassarian, H., Veneris, A.G., Liffiton, M.H., Sakallah, K.A., 2007. Improved design debugging using maximum satisfiability. In: Formal Methods in Computer-Aided Design. IEEE Computer Society, pp. 13-19.

Sinz, C., 2005. Towards an optimal CNF encoding of Boolean cardinality constraints. In: Principles and Practice of Constraint Programming. In: LNCS, vol.3709, Springer, pp. 827-831.

Soh, T., Banbara, M., Tamura, N., Le Berre, D., 2017. Solving multiobjective discrete optimization problems with propositional minimal model generation. In: International Conference on Principles and Practice of Constraint Programming. Springer, pp. 596-614.

Tamby, S., Vanderpooten, D., 2021. Enumeration of the nondominated set of multiobjective discrete optimization problems. INFORMS J. Comput. 33 (1), 72-85.

Tamura, N., Banbara, M., 2008. Sugar: A CSP to SAT translator based on order encoding. In: Proceedings of the Second International CSP Solver Competition. Citeaseer, pp. 65-69.

Terra-Neves, M., Lync, I., Manquinho, V., 2017. Introducing Pareto minimal correction subsets. In: International Conference on Theory and Applications of Satisfiability Testing. Springer, pp. 195-211.

Zitzler, E., Knowles, J.D., Thiele, L., 2008. Quality assessment of Pareto set approximations. In: Branke, J., Deb, K., Miettinen, K., Slowinski, R. (Eds.), Multiobjective Optimization, Interactive and Evolutionary Approaches [Outcome of Dagstuhl Seminars]. In: Lecture Notes in Computer Science, vol.5252, Springer, pp. 373-404.

Zitzler, E., Thiele, L., 1998. Multiobjective optimization using evolutionary algorithms - a comparative case study. In: Eiben, A.E., Bäck, T., Schoenauer, M., Schwefel, H. (Eds.), Parallel Problem Solving from Nature - PPSN V, 5th International Conference, Amsterdam, the Netherlands, September 27-30, 1998, Proceedings. In: Lecture Notes in Computer Science, vol.1498, Springer, pp. 292-304.

Zitzler, E., Thiele, L., Laumanns, M., Fonseca, C.M., da Fonseca, V.G., 2003. Performance assessment of multiobjective optimizers: An analysis and review. IEEE Trans. Evol. Comput. 7 (2), 117-132.