Network Models for Multiobjective Discrete Optimization
David Bergman∗1 , Merve Bodur†2 , Carlos Cardonha‡3 , and Andre A. Cire∗∗4

arXiv:1802.08637v1 [math.OC] 23 Feb 2018

1

Department of Operations and Information Management, University of Connecticut
2
Department of Mechanical and Industrial Engineering, University of Toronto
3
IBM Research
4
Department of Management, University of Toronto Scarborough

February 26, 2018

Abstract
This paper provides a novel framework for solving multiobjective discrete optimization problems with an arbitrary number of objectives. Our framework formulates these problems as
network models, in that enumerating the Pareto frontier amounts to solving a multicriteria
shortest path problem in an auxiliary network. We design tools and techniques for exploiting
the network model in order to accelerate the identification of the Pareto frontier, most notably
a number of operations to simplify the network by removing nodes and arcs while preserving
the set of nondominated solutions. We show that the proposed framework yields orders-ofmagnitude performance improvements over existing state-of-the-art algorithms on four problem
classes containing both linear and nonlinear objective functions.
Keywords. Multicriteria decision making; multiobjective discrete optimization; 0-1 programming; Pareto frontier; network models

1

Introduction

Multiobjective optimization is the study of algorithms for optimization problems associated with two
or more objective functions. Such problems arise naturally in prescriptive decision making, where
decision makers often face conflicting criteria that balance trade-offs associated with a proposed
solution. Practical applications of multiobjective optimization are pervasive and found in a diverse
set of domains, including, for example, production planning for supply chains (Dickersbach 2015),
radiotherapy optimization in healthcare (Yu et al. 2000), and aerodynamic design (Wang et al.
2011). These problems involve from a few up to potentially hundreds of objective functions, and
recent research has investigated using multiobjective optimization to solve single-objective integer
programming problems (Bodur et al. 2016), expanding even further on the importance of the study
of algorithms designed to tackle these computationally challenging problems.
A solution in multiobjective optimization typically consists of the set of points in the objective
space that belong to the Pareto frontier (Ehrgott 2006b) (i.e., points that improve upon at least one
∗

david.bergman@business.uconn.edu
bodur@mie.utoronto.ca
‡
carloscardonha@br.ibm.com
∗∗
acire@utsc.utoronto.ca
†

1

objective with respect to all feasible points). In practice, the Pareto frontier fully characterizes the
trade-offs among solutions and is post-processed in interactive decision support systems (Stewart
et al. 2008). This allows practitioners to prioritize objectives and operational aspects based on
their technical expertise, often on a case-by-case basis.
In this paper, we focus on the large class of multiobjective discrete optimization problems
(MODOs). There is a rich history of the study of tools and techniques for addressing MODOs
due to its relevance in industry. Most algorithms appearing in the literature leverage advances
in mathematical programming technology to solve this more general class of problems (see, e.g.,
in-depth surveys by Ehrgott 2006a and Zhou et al. 2011). Current state-of-the-art methodologies
rely on parametric single-objective reformulations of MODOs, employing commercial MILP solvers
as critical black-box components of their algorithms. The restriction to solutions whose variables
are all required to assume values from a discrete set may result in MODOs whose Pareto frontiers
cannot be fully recovered by such linear parametric models (Sayın and Kouvelis 2005). These
methods therefore require explicit enumeration techniques, which substantially limits the size of
problems that can be solved, both in terms of the number of variables and the number of objective
functions, especially when the Pareto frontier is large.
Contributions. This paper presents a fundamentally distinct approach for modeling and solving
a MODO. Our framework exploits the underlying discrete structure of a MODO by reformulating
the original problem as an optimization problem over a network that implicitly represent the Pareto
frontier. Each arc in the network is equipped with a multidimensional weight, and a point in the
objective space can be retrieved as the sum of the weights of the arcs in a given path of the network.
The full Pareto frontier can thus be enumerated by solving the classical multicriteria shortest path
problem (MSP) (Garroppo et al. 2010) over the network. Any network for which the MSP exactly
enumerates the Pareto frontier is referred to as a valid network.
Our first contribution is the design of two approaches for obtaining a valid network. One is
based on a recursive formulation of a MODO and the other is extracted from a decision diagram
representation of the problem (Bergman et al. 2016). Both techniques exploit construction procedures already available in the literature, some of which are inspired by early dynamic programming
(DP) approaches for multiobjective optimization (Carraway et al. 1990). We note that the proposed
approach is flexible in that it does not presuppose that the original MODO is modeled using DP,
mathematical programming, or any other restrictive modeling framework.
As our second contribution, we propose validity-preserving operations (VPOs) designed to reduce the size of a network while maintaining validity. A constructed network is, in general, exponentially large in the input size of the MODO. Even if the network itself is of manageable size,
computing the MSP may take a prohibitively long time. VPOs simplify the network without modifying the Pareto frontier, leading to significant reductions on the number of arcs and nodes and
hence computational time. We explore VPOs that are based solely on the network itself (e.g., removing arcs/nodes, merging nodes) and VPOs that explore domain-specific features of the problem
(e.g., using DP state-based information).
Finally, we present an extensive numerical study on four problem classes to compare our network model approach with two state-of-the-art MODO algorithms. We consider the knapsack, set
covering, and set partitioning problems, which are commonly used as benchmarks in the MODO
literature, in addition to a nonlinear MODO with absolute value functions in the objective. Our
experiments indicate that our approach solves MODO instances orders-of-magnitude faster than
the previous state-of-the-art algorithms, in addition to expanding on the size of problems that can
be solved, both in terms of the number of variables and the number of objective functions (up to
seven).
The benefits in adapting the approach developed in this paper are fourfold. First, network
2

representations do not require linear formulations, thereby broadening modeling expressiveness
over other popular techniques. Second, network models inherently exploit symmetry by combining
subpaths that correspond to common objective function evaluations, decreasing enumeration requirements. Third, whereas the design and execution of optimization algorithms were restricted by
highly-constrained memory environments in the past, modern-day computational platforms have
expansive memory resources, and the algorithms presented in this paper show how this new scenario
can be properly exploited in the design of modern algorithms for MODOs. Fourth, algorithms for
MODOs based on network models can immediately benefit from advances on the MSP, so continual
performance gains over the results reported in this work are expected in the future.
The remainder of this paper is organized as follows. §2 provides a literature review of MODO,
specifically as it relates to the present paper. §3 and §4 formally define MODOs and network
models, respectively. §5 describes network construction algorithms. §6 presents VPOs, and §7
presents the MSP algorithms we employ for enumerating the Pareto frontier. §8 describes the
results of an experimental evaluation on four problem classes. We conclude and describe future
work in §9. Proofs which do not appear in the main text are presented in the Appendix, §A.

2

Literature Review

There is an extensive literature on exact algorithms for generating the Pareto frontier of a MODO.
In general, these approaches can be divided into two main classes: those based on criterion-space
search, and those based on decision-space search.
Criterion-space search relies on scalarizations, most commonly based on a combination of
weighted sums and -constraints. Weighted-sum methods iteratively solve a single-objective optimization version of the problem where the single objective is defined by various positive-weight
combinations of the original objectives. For general MODOs, however, only a portion of the Pareto
frontier (those points referred to as supported efficient points) can be identified by this approach
alone. The remaining points (unsupported efficient points) can be found through the -constraint
method, introduced by Haimes et al. (1971), which optimizes one of the original objective functions
with the other constraints transformed into parametrized constraints. We refer to Ehrgott et al.
(2016) and Ehrgott (2006a) for technical details.
Kirlik and Sayın (2014) and Özlen et al. (2014) provide the state-of-the-art criterion-space
search algorithms for MODOs with an arbitrary number of objectives. Both algorithms are based
on variants of scalarization techniques. These variants build upon the early work by Klein and
Hannan (1982), who suggested an iterative approach to generate a subset of the Pareto frontier
while refining the search space by the addition of disjunctive constraints. Sylva and Crema (2004)
extended this work to an exact algorithm by reformulating the disjunctive conditions as big-M
constraints, which was further improved by Lokman and Köksalan (2013) and Bektas (2016).
Özlen and Azizoğlu (2009), in turn, developed an alternative approach called the augmented constraint method, which became one of the state-of-the-art methods after later enhancements by
Özlen et al. (2014). Another improvement to the augmented -constraint method is the recursive
methodology proposed by Laumanns et al. (2005, 2006), further refined by Kirlik and Sayın (2014)
into a computationally practical approach. We note that Boland et al. (2016b) developed an
extension of the so-called L-shape search method (specific to triobjective MODOs) that optimizes
a linear function over the set of nondominated points, that can also be used to enumerate the
Pareto frontier.
Other scalarization methods include the (lexicographic or augmented) weighted Tchebycheff
scalarization (Steuer and Choo 1983). The majority of the criterion-space search focuses on biob-

3

jective problems, where the special structure resulting from only having two objectives can be
exploited; see, e.g., Ralphs et al. (2006), Sayın and Kouvelis (2005), Boland et al. (2015), Parragh
and Tricoire (2015). Extensions of these ideas have also been proposed for triobjective MODOs,
which iteratively decompose the search space into smaller regions, and apply efficient ways to
explore and refine these regions (Dächert and Klamroth 2015, Boland et al. 2016a,c).
Decision-space search methods operate over the space defined by the original decision variables.
These techniques are typically based on branch-and-bound search developed for mixed-integer linear programs. The first of such algorithms was proposed by Mavrotas and Diakoulaki (1998) for
binary MODOs. The algorithm uses an artificial ideal point to define a bounding set, and discovers
nondominated points by solving (via a criterion-space search algorithm) the multiobjective linear
programs obtained when all binary variables are fixed. Mavrotas and Diakoulaki (2005) observed
that this branch-and-bound in fact generates a superset of the Pareto frontier, and proposed filtering algorithms to eliminate spurious points. Vincent et al. (2013) later showed that the previous
algorithm is still incomplete, and proposed a corrected and improved version for biobjective problems. Other branch-and-bound techniques have also been studied by Masin and Bukchin (2008),
Sourd and Spanjaard (2008), who suggested enhancements to the bounding aspects.
Stidsen et al. (2014) considered biobjective mixed-binary problems where continuous variables
appear in at most one objective function, and employed modified fathoming rules and a new branching paradigm they call Pareto branching. The balanced box method of Boland et al. (2015) can
also be used for this special case. This algorithm has recently been combined with ideas from
criterion-space search methods by Stidsen and Andersen (2017) to yield a hybrid algorithm. Gadegaard et al. (2016) extended Pareto branching and used it in their bound set-based branch-and-cut
algorithm.
The first decision-space algorithm for generic MODO based on branch-and-cut was developed
by Jozefowiez et al. (2012), where discrete sets are used for lower bounds. Recently, Adelgren
and Gupte (2017) developed a new branch-and-bound algorithm which employs multiobjective
extensions of many different aspects of branch-and-bound, such as (primal and dual) presolve,
preprocessing, node processing, and dual bounding via cutting planes.
Alternative methods that avoid scalarizations are based on DP and implicit enumerative methods for pure binary problems (e.g., Bitran 1977, Bitran and Rivera 1982). DP approaches are
typically focused on variants of the multiobjective knapsack problem, such as Villarreal and Karwan (1981) who employed lower and upper bound sets to eliminate dominated solutions. The
authors extended their work to general stagewise separable MODOs (Villarreal and Karwan 1982).
Klamroth and Wiecek (2000) presented distinct conceptual DP models for several variants of the
multiobjective knapsack problem, mostly based on extensions of the existing DP recursions for
the single-objective case. Bazgan et al. (2009) developed a new DP approach enhanced with complementary dominance relations for the 0-1 knapsack case. For the biobjective knapsack, Delort
and Spanjaard (2010) proposed a two-phase algorithm, and Rong and Figueira (2014) presented a
multiobjective DP algorithm.
Similar to the above DP approaches, the most relevant works to ours also considered the multiobjective knapsack problem. Captivo et al. (2003) proposed a transformation of biobjective 0–1
knapsack problem into a biobjective shortest path problem, which is solved via an enhanced version
of the labeling algorithm for MSPs. Figueira et al. (2010) developed a generic labeling algorithm for
the multiobjective integer knapsack problem applying existing reformulations from the literature.
More recently, state reduction techniques for the biobjective case have been proposed by Rong et al.
(2011) and Rong and Figueira (2013), with algorithmic enhancements by Figueira et al. (2013).

4

3

Multiobjective Discrete Optimization Problems

In this section we present the notation and formalism used throughout the text. For a ∈ N+ , we
let [a] := {1, 2, . . . , a} and, for any b ≥ a, b ∈ N, we let [a, b] := {a, a + 1, . . . , b}. Also, we represent
the vector of zeros and ones in an appropriate dimension by 0 and 1, respectively. The notation
B := {0, 1} indicates the Boolean set, while the operator (·)> denotes the transpose. For a given
vector v ∈ Rp , we denote the i-th component of v by vi or (v)i , for all i ∈ [p].
A multiobjective discrete optimization problem (MODO) is of the form
max {f (x) := (f1 (x), f2 (x), . . . , fK (x)) : x ∈ X } ,

(M)

where X is the feasible set and f : X → RK maps each solution x ∈ X into a corresponding Kdimensional objective vector (or image) y := f1 (x), f2 (x), . . . , fK (x) , with fk : X → R, k ∈ [K].
This paper focuses on binary optimization problems, i.e., we assume that X ⊆ Bn for some n ≥ 1,
although the techniques extend to MODOs with general integer variables. The objective functions
are not assumed to have any particular structure, except that they are well-defined in Bn .
For any two objective vectors y, y 0 ∈ RK , we say that y dominates y 0 , or simply y  y 0 , if (i)
yk ≥ yk0 for all k ∈ [K], and (ii) there exists at least one index k̃ for which yk̃ > yk̃0 . Conversely, y 0
is dominated by y, or y 0 ≺ y, if y  y 0 .
The image defined by the set of feasible solutions is denoted by Y := {f (x) : x ∈ X }. An
objective vector y ∗ ∈ Y is a nondominated point if there exists no other point y 0 ∈ Y for which
y 0  y ∗ . The set of all nondominated points of M is denoted by YN , also referred to as the Pareto
frontier. The typical goal of a MODO, and the focus of this paper, is to enumerate YN .
Example 1. We consider a set packing instance as a running example, with K = 3:
max

f1 (x) = 4x1 + 5x2 + 3x3 + 4x4 + 2x5 + 1x6 + 2x7 ,
f2 (x) = 8x1 + 7x2 + 1x3 + 5x4 + 3x5 + 3x6 + 8x7 ,

f3 (x) = 2x1 + 6x2 + 8x3 + 4x4 + 6x5 + 5x6 + 2x7

s.t. x1 + x2 + x3 ≤ 1, x2 + x3 + x4 ≤ 1, x4 + x5 ≤ 1,
x4 + x6 ≤ 1, x5 + x7 ≤ 1, x6 + x7 ≤ 1,
x ∈ B7 .
The nondominated set YN consists of the four points y 1 = (6, 7, 19), y 2 = (8, 13, 17), y 3 = (7, 14, 13),
and y 4 = (10, 21, 8). They are the images, respectively, of the feasible solution vectors x1 =
(0, 0, 1, 0, 1, 1, 0), x2 = (0, 1, 0, 0, 1, 1, 0), x3 = (1, 0, 0, 0, 1, 1, 0), and x4 = (1, 0, 0, 1, 0, 0, 1).

For any Ȳ ⊆ Y, let


ND Ȳ := y ∈ Ȳ : @ y 0 ∈ Ȳ with y 0  y
be an operator that returns the set of vectors within Ȳ that are not dominated by any other vector
in the same set. Note that YN = ND (Y). This operator has been studied in the context of relational
database systems, where it is known as the skyline operator (Borzsony et al. 2001). We refer to the
work by Gudala (2012) for a review of algorithms to compute ND (·) and their associated complexity
analysis. In particular, for K = 2, an efficient implementation of ND (S) for some S ⊆ RK has a
worst-case
time complexity
of O(|S| log (|S|)), and, for K > 2, it can be efficiently implemented in


K−2
O |S| · (log (|S|))
.

5

4

Network Models

This paper proposes the use of network models to enumerate YN for a MODO M. Formally, a
network model is a layered-acyclic multi-digraph N := (L, A) with node set L and arc set A. Such
a model is equipped with specific structure, properties, and attributesSassociated with M.
˙ Lj . Layers L1 and Ln+1
The node set L is partitioned into n + 1 non-empty layers L :=
j∈[n+1]

have cardinality one with L1 := {r} and Ln+1 := {t}. Nodes r and t are referred to as the root node
and the terminal node, respectively. The layer index of a node u ∈ L is ` (u), i.e., u ∈ L`(u) . Each
arc a := (r (a) , t (a)) ∈ A is directed from its arc-root r (a) ∈ Lj to its arc-terminal t (a) ∈ Lj+1 for
some j ∈ [n]. The layer of an arc is ` (a) := ` (r (a)). We denote by A+ (u) := {a ∈ A : r (a) = u}
the set of outgoing arcs from u and by A− (u) := {a ∈ A : t (a) = u} the set of incoming arcs to u.
We define P (u, v) as the set of arc-specified paths from node u to node v. The arguments
will be omitted when u = r and v = t, i.e., P := P (r, t). Each arc a has an arc-weight vector
w (a) ∈ RK . The arc-weight vectors provide the connection between arc-specifiedP
paths in N and
the objective space of M. Any path p = (aj1 , . . . , ajH ) has path-weight w (p) =
w (ajh ).
h∈[1,H]

The Pareto frontier of a network model N is defined as


[
PF (N ) := ND 
w(p) .
p∈P

A network model N is valid for a MODO M if PF (N ) = YN . In our structural results, we may
operate on distinct valid network models for the same M, in which case we append a subscript to
our notation so as to indicate the network. For example, rN will be used to represent the root node
of N . The subscript will be omitted when N is clear from context.
A network model can be interpreted as a data structure that supports the representation of a
MODO M as a MSP. By negating the arc-weight vectors (since we consider maximization), one
may apply any algorithm for solving MSPs to a network model N in order to find PF (N ) if N is
valid for M. Several MSP algorithms are available; see, e.g., the survey by Garroppo et al. (2010).
This paper provides a framework for building network models for MODOs, investigates techniques
for reducing a network model while preserving validity, and enhancements to MSP algorithms that
are specific to network models associated with MODOs arising in operations research applications.
Example 2. Figure 1 depicts a network model N for the MODO in Example 1. The arc-weight
vectors are shown (in black text) next to each arc (the additional details provided in the figure,
in red and blue, will be introduced and described in Example 6). There are 14 arc-directed paths
from r to t, and their path-weights are given by
[

w (p) = (8, 13, 17) , (6, 10, 11) , (7, 15, 8) , (6, 7, 19) , (4, 4, 13) , (5, 9, 10) , (3, 6, 11) ,
p∈P

(1, 3, 5) , (2, 8, 2) , (7, 14, 13) , (6, 13, 6) , (5, 11, 7) , (6, 16, 4) , (10, 21, 8)
and results in

PF (N ) = ND 


[

w(p) =



(8, 13, 17) , (6, 7, 19) , (7, 14, 13) , (10, 21, 8) .

p∈P

6

X (0, 0, 0)

r

x1

X (8, 13, 17)
(4, 8, 2)
X (7, 15, 8)
X (3, 6, 11)
X (6, 7, 19) X (4, 8, 2) u2
2 X (6, 13, 6)
× (5, 9, 10)
× (6, 13, 6)

(0, 0, 0)
u2
1

X (0, 0, 0)

x2

(5, 7, 6)

X (5, 6, 7) u3
1

x3

X (6, 7, 19)
X (5, 9, 10)
× (3, 6, 11)
X (0, 0, 0) u3
2
X (6, 13, 6)

(0, 0, 0)

(3, 1, 8)

X (5, 7, 6)
u4
1
X (3, 1, 8)

(0, 0, 0)
X (5, 7, 6)
X (3, 1, 8) u5
1
X (4, 8, 2)

x5

(2, 3, 6)
X (7, 10, 12)
X (5, 4, 14) u6
1
X (6, 11, 8)

x6

x7

X (4, 8, 2) u3
3

(0, 0, 0)

X (3, 6, 11)
X (2, 8, 2)

u7
1

X (3, 6, 11)
X (2, 8, 2) X (8, 13, 6) u5
2
× (1, 3, 5)

X (2, 8, 2)

(0, 0, 0)

(1, 3, 5)

X (0, 0, 0)

X (3, 6, 11)
× (2, 8, 2)
X (6, 13, 6)

(4, 5, 4)

(0, 0, 0)

X (5, 7, 6)
X (3, 1, 8) u6
2
X (4, 8, 2)

(1, 3, 5)

X (3, 6, 11)
X (6, 13, 6)

(0, 0, 0)

X (4, 8, 2)
u4
2
× (0, 0, 0)

(0, 0, 0)

X (1, 3, 5)

X (8, 13, 17)
X (6, 7, 19)
X (7, 14, 13)
× (6, 10, 11)
× (4, 4, 13)
× (5, 11, 7)

(0, 0, 0)

(0, 0, 0)

X (3, 6, 11)
X (2, 8, 2)

x4

X (8, 13, 17)
× (7, 15, 8)
X (6, 7, 19)
X (7, 14, 13)
X (10, 21, 8)

X (2, 8, 2)
X (1, 3, 5)

X (8, 13, 6) u6
3 X (2, 8, 2)

(0, 0, 0)
× (5, 7, 6)
X (3, 1, 8)
× (4, 8, 2)
X (8, 13, 6)

(0, 0, 0)
X (8, 13, 17)
X (6, 7, 19)
X (7, 14, 13)
t
× (5, 9, 10)
X (10, 21, 8)

(0, 0, 0)

u7
2 X (2, 8, 2)

(2, 8, 2)

X (0, 0, 0)

Figure 1: A valid network model for M in Example 1.

5

Network Model Construction

This section provides two frameworks for constructing valid network models for MODOs. The
first approach relies on recursive model of M. The second is a direct transformation from binary
decision diagrams, which is applicable when the objective functions are additively separable.

5.1

Recursive Formulations

Several classes of single-objective optimization problems admit recursive formulations, often written
as DP models. Bertsimas and Tsitsiklis (1997), for example, present guidelines on how to obtain
recursive models from linear programming formulations. These ideas were extended to the case
of multiple objectives in the early work by Villarreal and Karwan (1982), which we build upon
in this paper. In particular, while DP models are intrinsically associated with a state-transition
graph, we show that multiobjective recursive formulations are analogously associated with a valid
network model for the problem. Thus, we can leverage the reformulation and its associated network
structure to obtain more efficient algorithms for MODOs.
7

Formally, a multiobjective recursive formulation of a MODO M is written in terms of the
following elements which are depicted in Figure 2: (i) n+1 state variables s0 , s1 , . . . , sn ∈ S for some
state space S ⊆ Rm , where the initial state s0 is fixed a priori; (ii) n functions V1 , . . . , Vn : S → 2B
(where 2B denotes the powerset of B) that represent the variable-value assignments that can be
applied at a state (i.e., they provide state-dependent feasible decisions); (iii) n state transition
functions τ1 , . . . , τn : S × B → S, each mapping a pair (s, v) of a state s ∈ S and a variable-value
assignment v ∈ B to another state s0 ∈ S; and (iv) n reward functions δ1 , . . . , δn : S × B → RK
that map an analogous pair (s, v) to a reward vector in RK .
s0

x1 ∈ V1 (s0 )

s1 = τ1 (s0 , x1 )

δ1 (s0 , x1 )

x2 ∈ V2 (s1 )

s2 = τ2 (s1 , x2 )

xn ∈ Vn (sn−1 )

···

δ2 (s1 , x2 )

sn = τn (sn−1 , xn )

δn (sn−1 , xn )

Figure 2: Elements of a recursive formulation
Based on these four components, a multiobjective recursive problem is of the form
max
s,x

n
X

δj (sj−1 , xj )

(MR )

j=1

for all j ∈ [n],

s.t. sj = τj (sj−1 , xj ),
xj ∈ Vj (sj−1 ),

for all j ∈ [n].

In particular,
each objective function evaluation is represented by the k-th index of the total reward
P
tuple nj=1 δj (sj−1 , xj ) for k ∈ [K], i.e.,


n
X
fk (s, x) := 
δj (sj−1 , xj ) ,
j=1

k

while the feasible set of MR is

X := (s, x) ∈ S n+1 × Bn : sj = τj (sj−1 , xj ), xj ∈ Vj (sj−1 ) for all j ∈ [n] .
Assume, without loss of generality, that τn (s, v) = st for a fixed st and for all s ∈ S, v ∈ B.
That is, the final transition always leads to a common terminal state st , which can always be
accomplished by appropriately defining τn (·). If a MODO can be written in the form of MR , it
exposes a recursive structure that can be immediately leveraged for the construction of a network
model. More specifically, the network N = (L, A) is the state-transition graph defined by the
feasible set X , which is composed as follows (see, e.g., Bertsekas 2017 for algorithmic details):
• Node set L: The nodes of N are associated with nodes of the state-transition graph. For any
j ∈ [n + 1], there exists a node u ∈ Lj for every possible value of the state variable sj−1 in X .
More specifically, at layer j of N , there exists one node for each state in the set Projsj−1 (X ),
defined as the projection of the feasible set X into the space of the variable sj−1 . With a
slight abuse of notation, we let Lj = {uj1 , . . . , uj|Lj | } = Projsj−1 (X ). Note that since s0 and st
are fixed, layers L1 and Ln+1 are singletons; namely u11 = s0 and un+1
= st .
1
8

• Arc set A: The arcs of N are associated with arcs of the state-transition graph. In other
words, there exists an arc in N for every possible transition in the state-transition graph.
Again with a slight abuse of notation, we represent an arc as a triplet, rather than a pair,
(r (a) , t (a) , v) appending the variable-value information. Then, we have
A=

o
[
˙ n j j+1
(ui , ui0 , v) : uji ∈ Lj , uj+1
∈ Lj+1 , v ∈ Vj (uji ) with uj+1
= τj (uji , v) .
i0
i0
j∈[n]

• Arc weights: The reward functions provide the arc weights for N . If a = (uji , uj+1
i0 , v) ∈ A,
then,
w (a) = δj (uji , v).
Note that the weight of an arc does not depend on the state associated with its arc-terminal.
It follows from the definition of arc weights that a path weight corresponds to the objective
function in MR and, by construction of the state-transition graph, that there is a one-to-one
mapping between paths of N and solutions in X . Thus, N is a valid network model for MR .
Next, we provide a recursive formulation for our running example.
Example 3. The instance in Example 1 is a set packing instance, so it can be written as
o
n

max (c1 )> x, (c2 )> x, . . . , (cK )> x : Ax ≤ 1, x ∈ Bn
for an m × n matrix A with elements aij ∈ Bn , i ∈ [m], j ∈ [n], and cost vectors cj ∈ Rn , j ∈ [n].
A recursive formulation can be obtained as follows. A state variable s ∈ Bm is defined so that
for any i ∈ [m], si = 1 if and only if the i-th constraint of Ax ≤ 1 is satisfied as an equality or all
the decision variables appearing in the constraint have already been assigned to zero. Therefore,
we have the initial state s0 = 0, whereas any final transition will lead to the terminal state st = 1.
Variable-value assignments xj ∈ B represent the packing of element j. We are not allowed to set a
variable xj = 1 if any constraint which includes xj holds as an equality; i.e.,
Vj (s) := {v ∈ B : si + v ≤ 1 for all i ∈ [m] with aij = 1}.
The transition function ensures the consistency of the current state s and the next state s̃ when
the variable xj is set to the value v ∈ Vj (s). More specifically, let mi := argmax{j ∈ [n] : aij = 1}
be the maximum index of the variables with a nonzero coefficient in the i-th constraint. Then, we
have τj (s, v) = s̃ where

if j = mi
 1,
si + v, if j < mi and aij = 1
s̃i =
for all i ∈ [m],

si ,
otherwise
for any j ∈ [n]. Finally, the reward function for each j ∈ [n] is the contribution of xj to the
objective
function

 vector, which in this particular instance is state independent: δj (s, v) := v ×
c1j , c2j , . . . , cK
j .
For our particular instance in Example 1, this recursive formulation yields the network model
depicted in Figure 15 in the Appendix. We note that this network model matches the one given in
Figure 1 until layer six (inclusive). The network model in Figure 1 can be obtained from the one
in Figure 15 after applying Theorem 1 and Corollary 1 on layer seven.

9

5.2

Transformation from Binary Decision Diagrams

A binary decision diagram (BDD) is a network-based representation of a Boolean function (Bryant
1992), with a large list of applications in mathematical programming and operations research (see,
e.g., Behle 2007, Hadzic and Hooker 2007, Bergman et al. 2011). In such applications, BDDs are
used to compactly represent or approximate the set of feasible solutions to a discrete optimization
problem. We refer to the work by Bergman et al. (2016) for a survey of concepts and existing BDD
methodologies for optimization.
In our context, a BDD B is interpreted as a special type of network model where its paths
explicitly represent the feasible solutions of a discrete problem. Namely, B can be perceived as a
network model as defined in Section 4. However, B has no arc-weights, and it additionally associates
an arc-domain d (a) ∈ B with each arc a ∈ A. Each layer, except the last one, corresponds to a
unique decision variable mapped via the bijective function v : [n] → [n]; that is, layer l ∈ [n]
corresponds to variable xv(l) . Let PB be the set of paths from root to terminal in B. The arcdomains associate each path p = (a1 , . . . , an ) in PB with a solution vector x(p) ∈ Bn defined by
x(p)v(`(aj )) = d (aj ) for all j ∈ [n]. In particular, for a given feasible set X ⊆ Bn , B is said to be exact
for X if and only if x(p) ∈ X for all p ∈ PB and |PB | = |X |. That is, the paths of B encode X exactly.
Assume that, for a given MODO M, the objective function fk (·) is separable for all k ∈ [K];
i.e.,
X
fk (x) :=
fk,j (xj )
j∈[n]

for an appropriate choice of functions fk,j : B → R, j ∈ [n]. In that case, we can transform an
exact BDD into a network model N simply by defining suitable arc weights and removing the
arc-domain labels. More specifically, for every arc a ∈ A, we define its weight according to the
variable corresponding the layer of that arc, xv(`(a)) , and its assigned value, d (a), as

w (a) := f1,v(`(a)) (d (a)), . . . , fK,v(`(a)) (d (a)) .
The assignment of those weights to arcs implies that, for any path p = (a1 , . . . , an ) in the network,

X
w(p) =
f1,v(`(aj )) (d (aj )), . . . , fK,v(`(aj )) (d (aj )) = f (x(p)),
j∈[n]

thus PF (N ) is the Pareto frontier of M. This relation was first used in Bergman and Cire (2016).
There exists an extensive number of BDD construction algorithms for discrete optimization
problems with a (single) separable objective function (Bergman et al. 2016), some of which will
be employed for our numerical study in Section 8. Note that such construction algorithms can
be used in our case since we define multi-dimensional arc-weights after constructing the network
through a BDD, as explained above. They may either be constructed from a recursive formulation
of the problem, as in Section 5.1, or through methodologies that iteratively modify the diagram by
separating infeasible paths (Hadzic et al. 2008) or by intersecting distinct BDDs (Wegener 2000).
Note that in contrast with standard BDDs for optimization, our more general definition of a
network model is only concerned with the objective function image and does not include explicit
information about the solution vectors (see, e.g., Figure 1). Indeed, our network operators are
guaranteed only to preserve the Pareto frontier, and as such the diagram may be modified to
include infeasible paths or exclude feasible paths to allow for smaller network sizes (hence different
than an exact BDD).
As mentioned in Example 3, our running Example 1 is a set packing instance. §4 of Bergman
et al. (2016) provides a BDD representation for the set packing problem that can be readily transformed into a network model by applying the methodology described above.
10

6

Validity-preserving Operations

This section is devoted to the study of generic reduction algorithms on network models, enabling
MSP algorithms to run more efficiently. Formally, we study validity preserving operations (VPOs),
which are operations that transform a valid network model N into a new network model N 0 such
that PF (N ) = PF (N 0 ).

6.1

Weight Shifting and Node Merging

Because of the intrinsic connection between decision diagrams and network models, many operations
for single-objective BDDs can be extended for MODOs by exploiting the network structure. In this
section, we consider a classical operation denoted by reduction, which aims at eliminating nodes
while maintaining key properties of the graph.
Reduction was initially proposed by Bryant (1986) for merging nodes that share isomorphic subgraphs, thereby not affecting the representation of its underlying Boolean formula. Hooker (2013)
later extended this concept for BDDs encoding DPs with state-dependent rewards, showing that
further node merging can be achieved for a single-objective function by considering isomorphism
with respect to arc weights.
We now extend the concept of reduction proposed by Hooker (2013) to network models and
MODOs. We begin by defining a weight-shifting operation that, given a node u, transfers weights
from the arcs directed out of u (i.e., A+ (u)) to those directed into u (i.e., A− (u)).
For our results below, we assume that an initial valid network model N = (L, A) is available.
Proposition 1 (Weight shift). Let u ∈ L be a node in N such that u ∈
/ {r, t}. For any c̃ ∈ RK ,
the operation defined by setting
1. w (a) := w (a) − c̃ for all a ∈ A+ (u); and
2. w (a) := w (a) + c̃ for all a ∈ A− (u)
is a VPO.
Proof. Arc-weight w(p) remains unchanged for all p ∈ P and thus PF (N ) is unchanged.
Example 4. Consider the network model in Figure
 1. Let c̃ = (2, 8, 2). The network model
resulting from subtracting c̃ from arc-weight w u72 , t and adding c̃ to arcs w u62 , u72 and w u63 , u72
is a VPO. The last three layers of N are reproduced in Figure 3. The result of applying the weight
shift is depicted in Figure 4.

The main result of the section divulges a sufficient condition for when merging nodes is a VPO.
Theorem 1 (Node merge). Let u1 , u2 ∈ L, u1 6= u2 , be two nodes in N for which there exists a
one-to-one mapping between the arcs in A+ (u1 ) and A+ (u2 ) satisfying that, for every pair of arcs
(a1 , a2 ) such that a1 ∈ A+ (u1 ) maps to a2 ∈ A+ (u2 ), we have t (a1 ) = t (a2 ) and w(a1 ) = w(a2 ).
Then, the sequence of operations
1. Delete all arcs in A+ (u2 );
2. Redefine t (a) = u1 for all a ∈ A− (u2 ); and
3. Delete u2 ,
is a VPO.
11

u6
1

u6
2

u6
3

u6
1

u6
2

u6
3

u6
1

u6
2
(1, 3, 5)

(1, 3, 5)

(1, 3, 5)

(0, 0, 0)

u7
1

(0, 0, 0)

(0, 0, 0)

(1, 3, 5)

(1, 3, 5)

(2, 8, 2)

u7
2

u7
1

u7
2

(2, 8, 2)

(0, 0, 0)

(0, 0, 0)

t

Figure 3: Original arc-weights

(2, 8, 2)

t

Figure 4: Weight shift by
c̃ = (2, 8, 2) at node u72

u6
3
(2, 8, 2)

(1, 3, 5)

(2, 8, 2)

u7
1

(0, 0, 0)

t

Figure 5: Result of merge
operation on nodes u71 and u72

The proof of Theorem 1 is provided in Section A.1.
Example 5. Figure 5 depicts the result of merging nodes u71 with u72 in Figure 4. Note that
the resulting network model has fewer arcs and nodes, but the same number of paths. Such an
operation thereby decreases the size of the network model without altering PF (N ).

Proposition 1 and Theorem 1 define a strategy for simplifying a network model. In particular,
starting from the penultimate layer and moving upwards in the network, we process each node
in the layer in the following way. We construct the vector c̃(u) for each node u in the inspected
layer by taking the componentwise minimum arc-weight among arcs in A+ (u). Once c̃(u) has been
obtained, the arc-weights can then be shifted up as prescribed in Proposition 1. After applying
this operation to all nodes in a layer, the conditions of Theorem 1 and its VPO (when applicable)
are applied on the nodes of the same layer in the transformed network. Note that node-merge
operations can be performed in any order, in the sense that any sequence will lead to the same
reduced network.
Proposition 2. The componentwise minimum arc-weight shift and node merge VPO have worstcase time complexity of O(K|L|2 log(|L|)).
The proof of Proposition 2 is provided in Section A.2.

6.2

Arc removal

In this section, we investigate algorithms and structural results of VPOs that eliminate arcs of the
network model. For an arc a ∈ A, let N − a be the network model resulting from the removal of a
from N . By definition, removing arc a is a VPO if and only if PF (N ) = PF (N − a). The following
theorem shows that identifying when such condition holds in general is an NP-Hard problem.
Theorem 2. Given a valid network model N = (L, A) for a MODO M, let Ã ⊆ A be a subset of
arcs of N . Deciding whether there exists an arc a ∈ Ã such that PF (N ) = PF (N − a) is NP-hard.
Despite the hardness of determining whether an arc can be safely removed without changing
the Pareto frontier, we can still exploit strong sufficient conditions for designing arc-removal VPOs.
Additional notation is in order. Given a network model N = (L, A) and two nodes u, v ∈ L such
that ` (u) < ` (v), let N [u, v] be the network model containing only nodes and arcs in N that lie
on some path which starts at node u and ends at node v. We introduce the following concept:
Definition 1. A pair of nodes u, v ∈ L is isolating in N when, for every arc a ∈ AN \AN [u,v] ,
12

(i) t (a) ∈ LN [u,v] implies that t (a) = u; and
(ii) r (a) ∈ LN [u,v] implies that r (a) = v.
According to Definition 1, nodes u and v are isolating in N if N [u, v] contains all arcs from N
that are directed to nodes in LN [u, v]\{u} and all arcs that are directed out of nodes in LN [u, v]\{v}.
(For example, any pair of nodes are isolating in Figure 14(a).) Note that one can check whether
two nodes u and v are isolating in N in polynomial time, in the size N , with a breadth-first search.
Pairs of isolating vertices yield a sufficient condition for an arc-removal operation to be a VPO.
Theorem 3. Let u and v be isolating nodes in a network model N . For any a ∈ AN [u,v] , if
PF (N [u, v]) = PF (N [u, v] − a), then PF (N ) = PF (N − a). That is, if the removal of arc a is a
VPO in N [u, v], then it is also a VPO in N .
The proof of Theorem 3 is provided in Section A.4.
Theorem 3 shows that pairs of isolating nodes in an arbitrary network model N define subnetworks whose VPOs involving the removal of arcs are also VPOs for N . The simplest case reduces
to two arcs with the same endpoints, which yields the following immediate corollary.
Corollary 1. Let a1 and a2 be any two arcs of a network model N for which r (a1 ) = r (a2 ) and
t (a1 ) = t (a2 ). If w (a1 ) ≺ w (a2 ) or w (a1 ) = w (a2 ), then the removal of a1 is a VPO.
Theorem 3 can also be applied directly by choosing two nodes u, v such that N [u, v] is sufficiently
small so that all associated arcs can be removed efficiently. In particular, given δ := ` (v) − ` (u),
the number of paths in N [u, v] is bounded by O(2δ ), and hence for small δ the Pareto frontier (and
arc to be removed) can be identified quickly using, e.g., the procedures from Section 7, which we
discuss next. For our numerical evaluation, we fixed δ = 2. That is, (i) we find a pair (u, v) of
isolating nodes that are distant by at most two layers, (ii) obtain the network model N [u, v], (iii)
compute its Pareto frontier, and finally (iv) apply Theorem 3 to remove arcs that are VPOs in N .

7

On Finding the Pareto frontier

Given a valid network model, finding the Pareto frontier reduces to solving a MSP (by arc-weight
negation) in a layered-acyclic multi-digraph, for which an extensive literature exists (see, e.g.,
surveys by Tarapata 2007 and Garroppo et al. 2010).
In this section, we propose two algorithms (or compilation methods) for enumerating the Pareto
frontier based on our network model structure. The first is a direct modification of the unidirectional
recursion by Henig (1986), also applied, e.g., in Figueira et al. (2013) and Rong and Figueira 2014.
The second algorithm is an extension of Galand et al. (2013) and performs a bidirectional search
that combines the partial Pareto frontiers of each layer using a coupling operator. Both algorithms
assign a set (or a collection of sets) of K-dimensional vectors to nodes of the network. Each
K-dimensional vector is henceforth referred to as a label, as is done in the MSP literature.

7.1

Unidirectional Pareto frontier compilation

The unidirectional algorithms process one layer at a time, computing the partial nondominated
solutions at a node based on either the incoming arcs or the outgoing arcs. It is a direct application
of the recursion by Henig (1986) but using the underlying structure of the network model, similar
to the version presented by Rong and Figueira (2014).

13

The procedure works as follows. When processed from the root node r to the terminal node t,
the algorithm assigns a single set of top-down labels Z TD (u) to each node u. The label set of each
node is initialized as the empty set, except the root node, where Z TD (r) is initialized as {0}. For
each layer j from one to n, having constructed Z TD (u) for all u ∈ Lj , the labels are constructed for
the nodes in Lj+1 by considering
the arcs directed from nodes in Lj to Lj+1 , one by one. For each
S
such an arc a (i.e., a ∈ u∈Lj A+ (u)) and every label z ∈ Z TD (r (a)), the label z + w(a) is added
to Z TD (t (a)).
After all arcs directed out of nodes in Lj are processed, Z TD (u) is re-assigned to

ND Z TD (u) , to remove any labels that are dominated by other labels in the set. Note that one
can also do a simple check each time a label is added to see if it is dominated by another label
already existing for the node. At the culmination of the algorithm, Z TD (t) will be PF (N ).
One can also run the algorithm in the opposite direction, starting from t and flipping the
direction of the arcs. The terminal node is initialized as {0}, and the nodes are processed in the
opposite direction. We refer to the labels constructed in this direction as bottom-up labels Z BU (u).
The set of labels Z BU (r) coincides with Z TD (t) and is therefore equal to PF (N ).
Example 6. Consider Figure 1. The labels on the left of the nodes (shown in red) correspond to
the top-down labels (i.e., for every node u, they list Z TD (u)). A “X” is drawn next to labels that
remains in Z TD (u) after the application of the ND () operator, and a “×” is drawn otherwise. To
the right of each node u, Z BU (u) is listed (in blue), with symbols “X” and “×” indicating whether
the labels remain or are discarded after the application of the ND () operator, respectively.


7.2

Bidirectional Pareto frontier compilation

We now provide a compilation method that extends the work of Galand et al. (2013) for network
models. Namely, one may obtain the elements composing the Pareto frontier by constructing labels
in both directions simultaneously and coupling the top-down and bottom-up labels. Given two sets
of vectors Z1 , Z2 ⊆ RK , define the coupling of Z1 and Z2 as


CP (Z1 , Z2 ) := ND z : z = z 1 + z 2 , z 1 ∈ Z1 , z 2 ∈ Z2 .
The coupling function CP (Z1 , Z2 ) returns the nondominated set of vectors that result from every
pairwise sum of vectors from the two sets Z1 and Z2 .
Let us fix a layer j 0 and suppose we created the labels Z TD (u) for every node u ∈ Lj , j ≤ j 0 ,
and the labels Z BU (u) for every node u ∈ Lj , j ≥ j 0 . We define the operation of coupling on layer
j 0 as


[


CPL Lj 0 := ND 
CP Z TD (u) , Z BU (u)  .
u∈Lj 0

This yields the nondominated set that results from the coupling of the top-down and the bottom-up
labels on each node. Note that CPL (Lj ) = YN for any j ∈ [n + 1].
Example 7 shows that this approach can significantly reduce the number of operations required
to find the Pareto frontier of a network model. Since the nondominated
frontier of any set S

K−2
of K-dimensional vectors can be found in time O |S| · (log (|S|))
(Borzsony et al. 2001), the


coupling operation of sets Z1 and Z2 can be completed in time O |Z1 | · |Z2 | · (log (|Z1 | · |Z2 |))K−2 .
Example 7. Consider the network model in Figure 1. Suppose we fix L5 , composed of nodes
u51

and u52 , as the coupling layer. Only 14 top-down labels need to be created to find Z TD u51 and
14




Z TD u52 , and only 11 bottom-up labels need to be created in order to find Z BU u51 and Z BU u52 .
The coupling of these sets results in


CP Z TD u51 , Z BU u51
= {(8, 13, 17) , (6, 7, 19) , (7, 14, 13) , (7, 15, 8) , (6, 16, 4)}


CP Z TD u52 , Z BU u52
= {(10, 21, 8)} ,
and, finally,
CPL (L5 ) = {(8, 13, 17) , (6, 7, 19) , (7, 14, 13) , (10, 21, 8)} ,
as desired. Note that using either unidirectional approach requires the creation of a total of 36
labels, as opposed to the 25 required using coupling.

Determining the best layer to couple on (i.e., the one in which the number of labels that need
to be created is minimized) is a challenging optimization problem, even difficult to model concisely.
In particular, a unidirectional Pareto frontier compilation may require the creation of fewer labels
than the bidirectional variant. We therefore propose the following heuristic procedure. Starting
from r and t, we first create Z TD (u0 ) and Z BU (u00 ), respectively, for all u0 ∈ L2 and for all
u00 ∈ Ln . Then, having constructed top-down labels for each node up to Lj1 , 2 ≤ j1 , and bottomup labels for each node on or below Lj2 , j2 ∈ [j1 + 1, n], we pick among j1 and j2 the layer with
the
ofP
total labels in order to proceed with the extension operations. Namely, if
P fewer number
TD (u) | ≤
BU (u) |, extension of the top-down labels to L
|Z
j1 +1 is completed and
u∈Lj1
u∈Lj2 |Z
j1 is set to j1 + 1. Otherwise, extension of the bottom-up labels to Lj2 −1 is completed, and j2 is
set to j2 − 1. This procedure is repeated until j1 = j2 , upon which coupling on layer Lj1 is used to
calculate the Pareto frontier of the network model.

7.3

Label Removal Algorithms

Given a network model, the complexity of finding the Pareto frontier is largely determined by
the cardinality of Z TD (u) and Z BU (u). Thus, having a VPO for the reduction of |Z TD (u) | and
|Z BU (u) | is relevant for computational purposes.
The following proposition introduces such an operation. Given two nodes u, v in a same layer
Lj , the intuition behind the proposition is to identify when the subnetwork associated with N [u, t]
is dominated by N [v, t], in which case we can remove labels in u that are dominated by v.
Proposition 3 (Label filtering). Let u and v be two nodes in Lj for some j ∈ [n]. Suppose
ND (PF (N [u, t]) , PF (N [v, t])) = PF (N [v, t]) .
If there exists a pair of labels `u ∈ Z TD (u) and `v ∈ Z TD (v) for which `u ≺ `v (or `u = `v ), then
removing `u from Z TD (u) is a VPO. Similarly, if
ND (PF (N [r, u]) , PF (N [r, v])) = PF (N [r, v]) ,
and there exists `u ∈ Z BU (u) and `v ∈ Z BU (v) for which `u ≺ `v (or `u = `v ), then removing `u
from Z BU (u) is a VPO.
Proof. We provide a proof of the first case, as the proof of the other case follows by inverting the
network model. By the condition in the statement of the proof, for each path p in N[u,t] , it must
be that w(p) + `u is dominated by w(p0 ) + `v , for some path p0 in N[v,t] . Thus w(p) + `u does not
belong to PF (N ) for any such p. It follows that the removal of `u is a VPO.
Proposition 3 generalizes the concept of state-based dominance in DP (Bertsekas 2017) to network models. In particular, we can incorporate domain-specific information into a network model
to identify cases where the conditions of Proposition 3 are satisfied. We provide an example instantiation in Section 8 which results in significant speedups in enumerating the Pareto frontier.
15

8

Numerical Study

In this section, we provide a detailed numerical evaluation of the effectiveness of the proposed
algorithm on four different classes of problems. For each class, we discuss how the initial network
model is constructed, explain the best algorithmic configuration, and compare with existing stateof-the-art approaches for general MODOs. In particular, we consider the methodologies proposed
by Kirlik and Sayın (2014) and Özlen et al. (2014), hereafter denoted by algorithms K and O,
respectively. The source codes of these algorithms were kindly provided by the respective authors.
Experimental Design and Evaluation. All experiments ran on an Intel(R) Xeon(R) CPU
E5-2680 v2 at 2.80GHz. Each experiment was limited to one thread and subject to a time limit of
3,600 seconds and a memory limit of 16GB. The algorithms K and O depend on the resolution of
integer linear programs; we employed IBM ILOG CPLEX 12.7.1 with default settings for this task
(IBM ILOG 2017). The paired Wilcoxon signed-rank test was used to estimate p-values comparing
pairs of algorithms (Wilcoxon 1945) (this is a nonparametric test that will be used to compare if
population mean ranks of solution times differ between algorithms).
Direct comparisons between the algorithms for each class of problems are presented in cumulative distribution plots, which show the number of instances solved by each algorithm (y-axis)
within a given time limit (x-axis). We use the integral of the curve associated with an algorithm
to estimate its relative performance so that the relative quality of an algorithm depends on the
running time and on the number of instances solved. The order in which algorithms are listed in
the legends reflect this metric, with the best-performing algorithms appearing on the top.
Additionally, we present scatter plots to compare the best-performing network model algorithm
against the best previous state-of-the-art algorithm. These plots are presented in logarithmic
scale and represent the amount of time the algorithms require to solve each instance of the given
benchmark. We also employ a color code to indicate the size of the Pareto frontier of each instance.
The network model-based algorithms employing the bottom-up, top-down, and bidirectional
Pareto frontier compilation are represented by BU, TD, and Coup, respectively. For applications
where the domain-specific label filtering given in Proposition 3 has been applied, we denote the
extensions of BU, TD, and Coup by BU+, TD+, and Coup+, respectively.

8.1

Multiobjective 0-1 Knapsack Problem

Given n items, a capacity W > 0, and for each item i ∈ [n], a weight wi > 0 and K profits
p1i , p2i , . . . , pK
i > 0, the multiobjective 0-1 knapsack problem (MKP) is
n
o
max ((p1 )> x, (p2 )> x, . . . , (pK )> x) : w> x ≤ W, x ∈ Bn .
The MKP is a standard benchmark extensively used in the multiobjective literature.
Network model construction. The initial network is constructed via a multiobjective recursive
formulation as presented in Section 5.1. The recursive formulation is obtained by using a single
dimensional state variable s ∈ R+ , which corresponds to the total weight of the knapsack at a
certain stage. The root state is s0 = 0. We cannot set a variable xj = 1 if it weighs more than the
available capacity, i.e., Vj (s) := {v ∈ B : s + v · wj ≤ W }. The transition functions update the
total weight of the knapsack: τj (s, v) = s + v · wj for all j ∈ [n − 1] and τn (s, v) = W . Lastly, for

any j ∈ [n], the reward function is the profit vector of item j, i.e., δj (s, v) := v × p1j , p2j , . . . , pK
j ,

16

thus state-independent. Note that the transition function at the last layer ensures that there is a
unique terminal node t.
We can incorporate the label filtering of Proposition 3 by exploiting the classical DP state
dominance for knapsack problems. For any j ∈ [n], let uji , uji0 ∈ Lj be two possible states obtained
at layer j using the aforementioned recursive model. If uji ≥ uji0 , we have
 





ND PF N [uji , t] , PF N [uji0 , t] = PF N [uji0 , t] .
The equality above follows since each path (i.e., partial feasible solution) in N [uj , t] has a path
in N [uji0 , t] of same path-weight, given that we have more capacity available at uji0 than uji . That
is, we can assign the same variable values in a path starting at uji0 and incur the same objective
function contribution. Thus, we can remove labels at uji if they are dominated by a label at uji0 .
Data generation. We experimented on 450 knapsack instances with K ∈ {3, 4, . . . , 7} and n ∈
{20, 30, . . . , 100}. Instances were generated as in Kirlik and Sayın (2014): each profit pkj and weight
wj was randomly draw from the
Pninterval [1, 1000], for j ∈ [n] and k ∈ [K]. The capacity of the
knapsack was set to W := d0.5 i=1 wi e. Ten instances were generated for each (n, K) pair.
Computational evaluation. The detailed results are presented in Table 1 of the electronic
supplement. We provide an analysis of the results and summarize our findings.
350
300

Coup+ time (seconds)

Number solved

250
200
Coup+
Coup
TD
TD+
BU+
BU
O
K

150
100
50
0
0

900

1800

2700

140000

n=20
n=100

K=3
K=4
1000 K=5

120000
100000

100

80000

10

60000
1
40000
0.1
20000
0.01
0.01

3600

0.1

1

10

100

1000

0

O time (seconds)

Time (seconds)

Figure 6: Knapsack cumulative distribution plot

Figure 7: Knapsack scatter plot

The cumulative distribution plot for the knapsack instances is presented in Figure 6. The results
show a clear dominance of all network model algorithms over O and K. Overall, Coup+ delivered the
best results, solving 370 instances, whereas the configurations BU and BU+ were relatively weaker
and solved 280 and 279 and instances, respectively. O and K solved 154 and 149, respectively. The
figure also shows that the network model algorithms are considerably faster, as they solve more
instances within seconds than O and K in one hour.
Figure 7 shows a scatter plot comparing Coup+ with O for instances with up to five objective
functions. We removed from the plot instances with K = 6 and K = 7, since O and K have
considerably worse performance and the results do not provide much insight. Coup+ was at least as
efficient as O and K in every knapsack instance tested, and only three instances that could not be
solved by this network model configuration were solved by the others (one by TD and two by Coup).

17

The sizes of the Pareto frontier do not seem to affect the relative performance between Coup+ and O;
namely, Coup+ and the other network model algorithms perform better than the previous state of
the art in all cases, perhaps even more so in instances with smaller solution sets.
With respect to Pareto frontier compilation, the bidirectional strategy had the best results while
bottom-up was relatively poor, independently from the inclusion of the label filtering VPO. Filtering
affected network model algorithms in different ways, depending on the Pareto frontier compilation
strategy used. The solution time differences between BU and BU+ are statistically significant (pvalue of 10−9 ), although in practical terms they perform similarly; the number of solved instances
is almost the same (280 vs 279) and the average running time goes from 252 to 244 with filtering, a
gain of 3% on average. The inclusion of filtering decreased the quality of the top-down algorithm,
with 310 instances solved (as opposed to 323) and almost 40 seconds of additional computational
time, on average, to solve instances solved by both algorithms (from 184 to 221, with p-value of
10−13 ). Finally, Coup+ is significantly better than Coup (p-value of 10−40 ); more instances were
solved (370 in comparison with 362) in less time (average running time reduced from 318 to 228)
and less variability (standard deviation reduced from 710 to 525).

8.2

Multiobjective Set Covering and Set Partitioning Problems

We consider the multiobjective variants of the classical set covering and set partitioning problems.
Namely, let A ∈ Bm×n be a binary constraint coefficient matrix, and let c1 , . . . , cK be K cost vectors
in Rn . The multiobjective set covering problem (MSCP) is defined as
o
n
min ((c1 )> x, (c2 )> x, . . . , (cK )> x) : Ax ≥ 1, x ∈ Bn .
The multiobjective set packing problem (MSPP), in turn, is written as
o
n
max ((c1 )> x, (c2 )> x, . . . , (cK )> x) : Ax ≤ 1, x ∈ Bn .
Data generation. We generated 150 random instances based on the previous work by Stidsen
et al. (2014). Specifically, we considered n ∈ {100, 150, 200} variables, m = n/5 constraints, and
fixed 10 variables per constraint (i.e., for every k = 1, . . . , m, we pick 10 elements of the k-th
row of A to be 1 uniformly at random). We considered K = 3, . . . , 7 objectives. We sampled 10
instances per (n, K) pair. Each matrix was used both for the MSCP and the MSPP by selecting
the appropriate inequality direction in the ensuing constraint system.
Network model construction. The original networks are produced by employing the BDD
transformation discussed in Section 5.2. In particular, we used the set covering BDD by Bergman
et al. (2011) and the set packing BDD by Bergman et al. (2014). In our experiments, label filtering
for both applications did not impact performance, so we omit the corresponding results.
Computational evaluation. The detailed results for the MCSP and the MSPP are presented
in Tables 2 and 3, respectively.
For the MSCP, as shown in Figure 8, Coup delivered the best results among the network model
algorithms, solving 90 instances with an average running time of 81 seconds (and a standard
deviation of 188). TD solved more instances (91), but its average runtime was higher (122 with a
standard deviation of 263). Moreover, for the 88 instances solved by both, Coup had an average
running time of 64 (with a standard deviation of 127), against 104 (standard deviation of 229)
for TD.
18

9000
8000

1000
Coup time (seconds)

75
Number solved

10000

n=100
n=200

K=3
K=4
K=5

100

50

25

Coup
TD
BU
O
K

0
0

900

1800

2700

7000
6000

100

5000
4000
3000

10

2000
1000

1

3600

1

10

Time (seconds)

100

0

1000

O time (seconds)

Figure 8: Set covering cumulative distribution plot

Figure 9: Set covering scatter plot

Algorithms K and O solved 48 and 50 instances, respectively. Among these, there were instances
that the network model configurations could not enumerate the Pareto frontier (between 9 and 11,
depending on the configuration employed). These instances were relatively large, typically of size
n = 200, resulting in large network models that could not be solved within the given limits.
Figure 9 depicts the relative performance of O and Coup on the MSCP, in particular elucidating
the one configuration for which O significantly outperformed Coup (200 variables, K = 3). However,
the same plot also suggests that Coup is far more efficient for instances with relatively large Pareto
frontiers. This suggests that the performance of the network models are less sensitive to the size
of the Pareto frontier than K and O.
The results for the MSPP are presented in the cumulative distribution plot in Figure 10 and
in the scatter plot in Figure 11. Coup also delivered the best results for this problem class, slightly
K=3
K=4
1000 K=5

100

Coup time (seconds)

Number solved

75

50

25

Coup
TD
BU
K
O

0
0

900

1800

2700

25000

n=100
n=200

20000

100
15000
10
10000

1
0.1

5000

0.01

3600

0.01

Time (seconds)

0.1

1

10

100

1000

0

K time (seconds)

Figure 10: Set packing cumulative distribution plot

Figure 11: Set packing scatter plot

outperforming TD; Coup solved one more instances than TD (102 vs 101) and had smaller running
times (averages of 71 versus 123 and standard deviation of 170 versus 322). Coup and TD solved more
instances of the MSPP than the MSCP, whereas BU, K, and O had the inverse behavior. Algorithms
K and O again solved fewer instances (41 and 42) and, among these, between 9 and 10 instances
(depending on the configuration) were not solved by the network models. The size of the instances
19

played a role in the efficiency of the algorithms. Figure 11 shows that K outperforms Coup in some
instances, but in relatively fewer cases than in the MSCP. Moreover, the Pareto frontier sizes have
the same impact on the relative performance of the algorithms, as in the case of the MSCP, showing
the robustness of network models.

8.3

Multiobjective Cardinality-Constrained Absolute Value Problem

The multiobjective cardinality-constrained absolute value problem (MCCAVP) is defined as
n

o
min
(a1 )> x − b1 , (a2 )> x − b2 , . . . , (aK )> x − bK : 1> x ≤ C, x ∈ Bn ,
where a1 , . . . , aK ∈ Zn , b ∈ ZK , and C ∈ Z+ . The MCCAVP is a multiobjective variant of the
discrete L1 -norm minimization problems, classically applied in statistical data fitting and circuit
optimization (Jong 2012, Narula and Wellington 1982). For instance, in data fitting problems each
linear function represents a residual, and in the multiobjective case we wish to evaluate the Pareto
frontier of nondominated fitting configurations according to each residual.
The MCCAVP illustrates how the procedure generalizes to multiobjective nonlinear applications. Also note that, if any of the K objective functions is instead written as a linear function
raised to the power of α ∈ Z≥1 (as opposed to the absolute value as in the MCCAVP), the outer
function can be replaced by the absolute value (if α is even) or simply ignored (if α is odd) without affecting the Pareto frontier. The MCCAVP therefore provides a modeling framework for a
wide-range of objective functions.
Data generation. We experimented on 6,250 instances of the problem. Each was randomly generated with K ∈ {3, 4, 5, 6, 7} and n ∈ {5, 10, 15, 20, 25}. For all k ∈ [K], we drew the components
of ak uniformly at random from the set [−M, M ], where we considered M ∈ {50, 100, 150, 200, 250}.
We let bk = b1> ak /2c for all k ∈ [K]. For the cardinality constraints, we let C = bnδc for which
we considered δ ∈ {0.1, 0.2, 0.3, 0.4, 0.5}.
Network model construction. The initial network is constructed via a multiobjective recursive
formulation as presented in Section 5.1. The recursive formulation is obtained by using a (K + 1)dimensional state variable s := (θ, γ) ∈ RK × R, where θ1 , . . . , θK represent the partial evaluation
of each (ak )> x for all k ∈ [K], and γ is the number of variables that are set to one at that stage.
The root state is s0 = (b, 0). We cannot set a variable xj = 1 if it exceeds the available capacity,
i.e., Vj (s) := {v ∈ B : γ + v ≤ C}. The transition functions update the partial evaluation of (ak )> x
and the number of variables that are set to one; i.e.,
τj (s, v) = (θ1 + a1j · v, θ2 + a2j · v, . . . , θK + aK
j · v, γ + v).
for all j ∈ [n]. The reward function is the change in each objective function when transitioning
from state s = (θ, γ) to another, that is, the k-th component of δj (s, v) is given by


δj (s, v) := |θk + akj · v − bk | − |θk − bk |,
k

for k ∈ [k]. To verify its validity, fix k ∈ [K] and consider any feasible solution x = (v1 , . . . , vn )
and the associated state transitions s0 = (θ0 , γ 0 ), s1 = (θ1 , γ 1 ), . . . , sn = (θn , γ n ). Observe that, by

20

definition,

n 
X
δj (sj−1 , vj ) = | θk0 + ak1 · v1 −bk | − |θk0 − bk | + | θk1 + ak2 · v2 −bk | − |θk1 − bk |
| {z }
| {z }
k
j=1

θk1

θk2

+ | θk2 + ak3 · v3 −bk | − |θk2 − bk | + · · · + | θkn−1 + akn · vn −bk | − |θkn−1 − bk |
| {z }
{z
}
|
θkn

θk3

= |θkn − bk | − |θk0 − bk |
| {z }
0

=

n
X

akj vj − bk ,

j=1

which is the original objective for the k-th function.
Mixed-Integer Linear Programming Reformulation. Since the original formulation for the
MCCAVP is nonlinear, it cannot be directly input to K and O. We consider the following linear
reformulation of the MCCAVP for K and O:
min {y1 , . . . , yK }

x∈Bn

s.t. yk ≥ (ak )> x − bk ,
k >

yk ≥ −(a ) x + bk ,

k ∈ [K],
k ∈ [K],

1> x ≤ C.
Computational evaluation. The detailed results for the MCCAVP are presented in Table 4.
In particular, the 6,250 instances were all solved by each network-based configuration. For this
section, we restrict the analysis of the results to 450 of these instances, which have M = 250,
n ≥ 15, and C ≥ 30, as the others were solved within a few seconds.
K=3
K=4
1000 K=5

300

BU time (seconds)

Number solved

400

200
BU
TD
Coup
K
O

100

0
0

900

1800

2700

1800

n=15
n=25

1600
1400

100

1200
1000

10

800
1

600
400

0.1

200
0.01

3600

0.01

Time (seconds)

0.1

1

10

100

1000

0

K time (seconds)

Figure 12: Absolute value cumulative distribution plot Figure 13: Absolute value scatter plot
The cumulative distribution plot is presented in Figure 12. Whereas K and O solved less than
200 instances each, the network algorithms enumerated the Pareto frontier in less than 10 minutes
in each case. Over the complete benchmark set, BU delivered the best results. For the restricted
21

set of 450 instances, BU was also the best, although its results (average running time of 18 seconds
with standard deviation of 41) were not different from those obtained by TD (average of 24 and
standard deviation of 69) in a statistically significant way (p-value of 0.026).
Algorithm O solved more instances than K in the extended dataset (4,118 vs. 4,108), albeit with a
higher runtime (on average 30% larger). Alternatively, for the restricted family of instances, K solved
more instances (178 vs. 172) with a much shorter runtime (average of 258 and standard deviation
of 492 vs. average of 432 and standard deviation of 801), thus suggesting that K outperforms O for
harder instances. We therefore select K for further comparison.
Figure 13 shows a scatter plot comparing the performance of BU with K. Algorithm K outperforms BU in 5 instances by one order of magnitude (2 vs 30 seconds was the largest relative
difference). All instances solved by K were also successfully addressed by BU, with an average running time of 8 seconds and standard deviation of 17. Furthermore, similar to what was observed
for the MSCP and MSPP, there is a positive correlation between the size of the Pareto frontiers
and the relative superiority of network model procedures over K and O.

9

Conclusions and Insights

This paper presents a novel framework for solving multiobjective discrete optimization problems
(MODOs) through a reformulation into network models, enhanced by validity-preserving operations
that reduce the size of the network while preserving the Pareto frontier. The framework is applied
across a wide-range of problem classes, showing significant improvement over two previous stateof-the-art algorithms. The generality of the framework is established through application to a
nonlinear multiobjective optimization problem.
The experimental evaluation provides a strong indication that the proposed algorithm outperforms exisiting state-of-the-art MODO solvers in several multiobjective variants of classical
operations research problems. The results also indicate that, in comparison to the other algorithms
tested, the proposed methodology runtime is less sensitive to the Pareto frontier size, further elucidating the computational advantage of this approach.
Our methodology assumes memory availability exceeds memory requirements for constructing
and storing network models. Since multiobjective optimization problems become rapidly more
challenging as problem size grows, this was only a limiting factor for the largest of instances
generated, most of which were beyond the scope of other algorithms we tested against. As memory
availability in modern-day CPUs continues to grow, investigating algorithms specifically designed
to exploit this memory availability is of great interest, and this paper provides an aimed step in
this direction for MODO.

References
N. Adelgren and A. Gupte. Branch-and-bound for biobjective mixed integer programming, 2017. https:
//arxiv.org/pdf/1709.03668.pdf.
C. Bazgan, H. Hugot, and D. Vanderpooten. Solving efficiently the 0-1 multi-objective knapsack problem.
Comput. and Oper. Res., 36(1):260–279, 2009.
M. Behle. On threshold bdds and the optimal variable ordering problem. J. Comb. Optim., 16(2):107–118,
2007.
T. Bektas. Disjunctive programming for multiobjective discrete optimisation, 2016. http://www.
optimization-online.org/DB_FILE/2016/05/5442.pdf.
D. Bergman and Andre A. Cire. Multiobjective optimization by decision diagrams. In Michel Rueher,
editor, Principles and Practice of Constraint Programming: 22nd International Conference, CP 2016,

22

Toulouse, France, September 5-9, 2016, Proceedings, pages 86–95. Springer International Publishing,
Cham, 2016.
D. Bergman, W.-J. van Hoeve, and J.N. Hooker. Manipulating MDD Relaxations for Combinatorial Optimization. In Tobias Achterberg and J. Beck, editors, Integration of AI and OR Techniques in Constraint
Programming for Combinatorial Optimization Problems, volume 6697 of Lecture Notes in Computer
Science, chapter 5, pages 20–35. Springer Berlin / Heidelberg, Berlin, Heidelberg, 2011.
D. Bergman, A. A. Cire, W.-J. Jan van Hoeve, and T. H. Yunes. BDD-based heuristics for binary optimization. Journal of Heuristics, 20(2):211–234, 2014.
D. Bergman, A.A. Cire, W.J. van Hoeve, and J. Hooker. Decision diagrams for optimization. Springer, 2016.
D.P. Bertsekas. Dynamic Programming and Optimal Control, volume 1. Athena Scientific, 4th edition, 2017.
D. Bertsimas and J.N. Tsitsiklis. Introduction to Linear Optimization, volume 6. Athena Scientific, 1997.
G.R. Bitran. Linear multiple objective programs with zero-one variables. Math. Program., 13(1):121–139,
1977.
G.R. Bitran and J.M. Rivera. A combined approach to solve binary multicriteria problems. Naval Res. Logist.,
29(2):181–201, 1982.
M. Bodur, S. Ahmed, N. Boland, and G.L. Nemhauser. Decomposition of loosely coupled integer programs: A
multiobjective perspective, 2016. http://www.optimization-online.org/DB_FILE/2016/08/5599.
pdf.
N. Boland, H. Charkhgard, and M. Savelsbergh. A criterion space search algorithm for biobjective integer
programming: The balanced box method. INFORMS J. Comput., 27(4):735–754, 2015.
N. Boland, H. Charkhgard, and M. Savelsbergh. The L-shape search method for triobjective integer programming. Math. Program. Comput., 8(2):217–251, 2016a.
N. Boland, H. Charkhgard, and M. Savelsbergh. A new method for optimizing a linear function over the
efficient set of a multiobjective integer program. European J. Oper. Res., 2016b.
N. Boland, H. Charkhgard, and M. Savelsbergh. The quadrant shrinking method: A simple and efficient
algorithm for solving tri-objective integer programs. European J. Oper. Res., 2016c.
S. Borzsony, D. Kossmann, and K. Stocker. The skyline operator. In Data Engineering, 2001. Proceedings.
17th International Conference on, pages 421–430. IEEE, 2001.
R. E. Bryant. Graph-based algorithms for boolean function manipulation. IEEE Transactions on Computers,
C-35:677–691, 1986.
R. E. Bryant. Symbolic boolean manipulation with ordered binary decision diagrams. ACM Computing
Surveys, 24:293– 318, 1992.
M.E. Captivo, J. Clımaco, J. Figueira, E. Martins, and J.L. Santos. Solving bicriteria 0-1 knapsack problems
using a labeling algorithm. Comput. and Oper. Res., 30(12):1865–1886, 2003.
R.L. Carraway, T.L. Morin, and H. Moskowitz. Generalized dynamic programming for multicriteria optimization. European J. Oper. Res., 44(1):95–104, 1990.
K. Dächert and K. Klamroth. A linear bound on the number of scalarizations needed to solve discrete
tricriteria optimization problems. J. Global. Opt., pages 1–34, 2015.
C. Delort and O. Spanjaard. Using bound sets in multiobjective optimization: Application to the biobjective
binary knapsack problem. Lecture Notes in Computer Science, 6049:253–265, 2010.
J.T. Dickersbach. Supply Chain Management with APO: Structures, Modelling Approaches and Implementation of MySAP SCM 4.1. Springer, Berlin, 2 edition, 2015.
M. Ehrgott.
A discussion of scalarization techniques for multiple objective integer programming.
Ann. Oper. Res., 147(1):343–360, 2006a.
M. Ehrgott. Multicriteria optimization. Springer Science & Business Media, 2006b.
M. Ehrgott, X. Gandibleux, and A. Przybylski. Exact methods for multi-objective combinatorial optimisation. In Multiple Criteria Decision Analysis, pages 817–850. Springer, 2nd edition, 2016.
J.R. Figueira, G. Tavares, and M.M. Wiecek. Labeling algorithms for multiple objective integer knapsack
problems. Comput. and Oper. Res., 37(4):700–711, 2010.

23

J.R. Figueira, L. Paquete, M. Simões, and D. Vanderpooten. Algorithmic improvements on dynamic programming for the bi-objective {0, 1} knapsack problem. Comp. Opt. and Appl., 56(1):97–111, 2013.
S. Gadegaard, M. Ehrgott, and L. Nielsen. Bi-objective branch–and–cut algorithms: Applications to the
single source capacitated facility location problem, 2016. http://www.optimization-online.org/
DB_FILE/2016/04/5402.pdf.
L. Galand, A. Ismaili, P. Perny, and O. Spanjaard. Bidirectional preference-based search for state space
graph problems. In Proceedings of the Sixth Annual Symposium on Combinatorial Search, SOCS 2013,
Leavenworth, Washington, USA, July 11-13, 2013., 2013.
R. G. Garroppo, S. Giordano, and L. Tavanti. A survey on multi-constrained optimal path computation:
Exact and approximate algorithms. Computer Networks, 54(17):3081 – 3107, 2010.
S.G. Gudala. Skyline queries for multi-criteria decision support systems. Master’s thesis, Kansas State
University, 2012.
T. Hadzic and J. N. Hooker. Cost-bounded binary decision diagrams for 0-1 programming. In E. Loute
and L. Wolsey, editors, Proceedings of the International Workshop on Integration of Artificial Intelligence and Operations Research Techniques in Constraint Programming for Combinatorial Optimization
Problems (CPAIOR 2007), volume 4510 of Lecture Notes in Computer Science, pages 84–98. Springer,
2007.
T. Hadzic, J. N. Hooker, B. O’Sullivan, and P. Tiedemann. Approximate compilation of constraints into
multivalued decision diagrams. In P. J. Stuckey, editor, Principles and Practice of Constraint Programming (CP 2008), volume 5202 of Lecture Notes in Computer Science, pages 448–462. Springer,
2008.
Y.Y. Haimes, L.S. Lasdon, and D.A. Wismer. On a bicriterion formulation of the problems of integrated
system identification and system optimization. IEEE Trans. Syst. Man and Cybern., 1(1):296–297,
1971.
M. I. Henig. The shortest path problem with two objective functions. European J. Oper. Res., 25(2):281 –
291, 1986.
J.N. Hooker. Decision diagrams and dynamic programming. In Carla Gomes and Meinolf Sellmann, editors,
Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization
Problems: 10th International Conference, CPAIOR 2013, Yorktown Heights, NY, USA, May 18-22,
2013. Proceedings, pages 94–110. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
IBM ILOG. CPLEX Optimization Studio 12.7.1 user manual, 2017.
Y. Jong. An augmented smoothing method of l1 -norm minimization and its implementation by neural
network model, 2012. https://arxiv.org/abs/1207.1931.
N. Jozefowiez, G. Laporte, and F. Semet. A generic branch-and-cut algorithm for multiobjective optimization
problems: Application to the multilabel traveling salesman problem. INFORMS J. Comput., 24(4):
554–564, 2012.
G. Kirlik and S. Sayın. A new algorithm for generating all nondominated solutions of multiobjective discrete
optimization problems. European J. Oper. Res., 232(3):479–488, 2014.
K. Klamroth and M.M. Wiecek. Dynamic programming approaches to the multiple criteria knapsack problem. Naval Res. Logist., 47(1):57–76, 2000.
D. Klein and E. Hannan. An algorithm for the multiple objective integer linear programming problem.
European J. Oper. Res., 9(4):378–385, 1982.
M. Laumanns, L. Thiele, and E. Zitzler. An adaptive scheme to generate the pareto front based on the
epsilon-constraint method. In Dagstuhl Seminar Proceedings. Schloss Dagstuhl-Leibniz-Zentrum für
Informatik, 2005.
M. Laumanns, L. Thiele, and E. Zitzler. An efficient, adaptive parameter variation scheme for metaheuristics
based on the epsilon-constraint method. European J. Oper. Res., 169(3):932–942, 2006.
B. Lokman and M. Köksalan. Finding all nondominated points of multi-objective integer programs.
J. Global. Opt., 57(2):347–365, 2013.

24

M. Masin and Y. Bukchin. Diversity maximization approach for multiobjective optimization. Oper. Res., 56
(2):411–424, 2008.
G. Mavrotas and D. Diakoulaki. A branch and bound algorithm for mixed zero-one multiple objective linear
programming. European J. Oper. Res., 107(3):530–541, 1998.
G. Mavrotas and D. Diakoulaki. Multi-criteria branch and bound: A vector maximization algorithm for
mixed 0-1 multiple objective linear programming. Appl. Math. Comput., 171(1):53–71, 2005.
S.C. Narula and J.F. Wellington. The minimum sum of absolute errors regression: A state of the art survey.
International Statistical Review / Revue Internationale de Statistique, 50(3):317–326, 1982.
M. Özlen and M. Azizoğlu. Multi-objective integer programming: a general approach for generating all
non-dominated solutions. European J. Oper. Res., 199(1):25–35, 2009.
M. Özlen, B.A. Burton, and C.A.G. MacRae. Multi-objective integer programming: An improved recursive
algorithm. J. Optim. Theory Appl., 160(2):470–482, 2014.
S.N. Parragh and F. Tricoire. Branch-and-bound for bi-objective integer programming, 2015. http://www.
optimization-online.org/DB_FILE/2014/07/4444.pdf.
T.K. Ralphs, M.J. Saltzman, and M.M. Wiecek. An improved algorithm for solving biobjective integer
programs. Ann. Oper. Res., 147(1):43–70, 2006.
A. Rong and J.R. Figueira. A reduction dynamic programming algorithm for the bi-objective integer knapsack problem. European J. Oper. Res., 231(2):299–313, 2013.
A. Rong and J.R. Figueira. Dynamic programming algorithms for the bi-objective integer knapsack problem.
European J. Oper. Res., 236(1):85–99, 2014.
A. Rong, J.R. Figueira, and M.V. Pato. A two state reduction based dynamic programming algorithm for
the bi-objective 0–1 knapsack problem. Comp. & Math. with Appl., 62(8):2913–2930, 2011.
S. Sayın and P. Kouvelis. The multiobjective discrete optimization problem: A weighted min-max two-stage
optimization approach and a bicriteria algorithm. Management Sci., 51(10):1572–1581, 2005.
F. Sourd and O. Spanjaard. A multiobjective branch-and-bound framework: Application to the biobjective
spanning tree problem. INFORMS J. Comput., 20(3):471484, 2008.
R.E. Steuer and E.U. Choo. An interactive weighted tchebycheff procedure for multiple objective programming. Math. Program., 26(3):326–344, 1983.
T. Stewart, O. Bandte, H. Braun, N. Chakraborti, M. Ehrgott, M. Göbelt, Y. Jin, H. Nakayama, S. Poles,
and D. Di Stefano. Real-world applications of multiobjective optimization. Multiobjective Optimization,
pages 285–327, 2008.
T. Stidsen and K.A. Andersen. A hybrid approach for Bi-Objective Optimization, 2017. http://http:
//www.optimization-online.org/DB_FILE/2017/04/5940.pdf.
T. Stidsen, K.A. Andersen, and B. Dammann. A branch and bound algorithm for a class of biobjective
mixed integer programs. Management Sci., 60(4):1009–1032, 2014.
J. Sylva and A. Crema. A method for finding the set of non-dominated vectors for multiple objective integer
linear programs. European J. Oper. Res., 158(1):46–55, 2004.
Zbigniew Tarapata. Selected multicriteria shortest path problems: an analysis of complexity, models and
adaptation of standard algorithms. International Journal of Applied Mathematics and Computer Science, 17(2):269287, 2007. ISSN 1641-876X.
B. Villarreal and M.H. Karwan. Multicriteria integer programming: A (hybrid) dynamic programming
recursive approach. Math. Program., 21(1):204–223, 1981.
B. Villarreal and M.H. Karwan. Multicriteria dynamic programming with an application to the integer case.
J. Optim. Theory Appl., 38(1):43–69, 1982.
T. Vincent, F Seipp, S Ruzika, A Przybylski, and X Gandibleux. Multiple objective branch and bound
for mixed 0-1 linear programming: Corrections and improvements for biobjective case. Comput. and
Oper. Res., 40(1):498–509, 2013.
W. Wang, R. Mo, and Y. Zhang. Multi-objective aerodynamic optimization design method of compressor
rotor based on isight. Procedia Engineering, 15(Supplement C):3699 – 3703, 2011. CEIS 2011.

25

I. Wegener. Branching Programs and Binary Decision Diagrams: Theory and Applications. SIAM monographs on discrete mathematics and applications. Society for Industrial and Applied Mathematics,
2000.
F. Wilcoxon. Individual comparisons by ranking methods. Biometrics bulletin, 1(6):80–83, 1945.
Y. Yu, J.B. Zhang, G. Cheng, M.C. Schell, and Paul Okunieff. Multi-objective optimization in radiotherapy:
applications to stereotactic radiosurgery and prostate brachytherapy. Artificial Intelligence in Medicine,
19(1):39 – 51, 2000. Evolutionary Computation in Medicine.
A. Zhou, B.-Y. Qu, H. Li, S.-Z. Zhao, P. N. Suganthan, and Q. Zhang. Multiobjective evolutionary algorithms: A survey of the state of the art. Swarm and Evol. Comput., 1(1):32 – 49, 2011.

A

Proofs

In this section we present the proofs for the structural results presented in the main text.

A.1

Proof of Theorem 1

Let N 0 be the network model that results from the application of the three operations stated in
Theorem 1 on nodes u1 , u2 of N . We claim that there exists a one-to-one mapping between pathweights in PN and path-weights in PN 0 ; i.e., for each path in PN , there is a path in PN with same
path-weight and vice-versa.
First, every path p ∈ PN that does not cross u2 remains unchanged in PN 0 , thus the equivalence
follows in this case. Let us assume now that p = (a1 , . . . , ai , ai+1 , . . . , an ) is a path in PN such that
t (ai ) = r (ai+1 ) = u2 . After the sequence of operations specified in the statement of the theorem,
0 ) = w(a ),
arc ai is substituted for an arc a0i in N 0 such that r (a0i ) = r (ai ), t (a0i )= u1 , and w(a
i
i

whereas arc ai+1 is substituted for an arc a0i+1 in N 0 such that r a0i+1 = u1 , t a0i+1 = t (ai+1 ),
and w(a0i+1 ) = w(ai+1 ). These operations yield a one-to-one relationship between p and the path
p0 = (a1 , . . . , a0i , a0i+1 , . . . , an ) in N 0 which, by construction, is such that w(p0 ) = w(p). Finally, as
no additional paths are introduced, each path in PN 0 is also a path in PN , so the result follows.

A.2

Proof of Proposition 2

These weight-shift and node-merge operations can be implemented in a bottom-up iterative procedure as follows.
Starting from the penultimate layer, we construct the vector c̃(u) for each node u ∈ Ln by
taking the componentwise minimum arc-weight among those arcs in A+ (u). Once c̃(u) has been
obtained, the arc-weights can then be shifted up according to the VPO described in Proposition 1.
Each arc is inspected twice and each node is visited once; thus, this operation can be performed
over the complete network in time O(K(|L| + |A|)).
Checking whether nodes u1 , u2 in layer Li satisfy the merging conditions of Theorem 1 is an
operation that can be performed in time O(K|Li+1 |) if an adequate data structure is employed.
Let π be an (arbitrary) ordering associated with the nodes of layer Li+1 and let πj denote the j-th
node of Li+1 according to this ordering, with the first position starting from zero. Given π, one
can construct a K|Li+1 |-dimensional vector ξ(u) associated with each node u of layer Li such that
(ξ(u))Kj+k = w((u, πj ))k , i.e., entry Kj + k contains the value of the k-th coordinate of the arcweight associated with the arc connecting u to vertex πj in Li+1 ; some arbitrary value, such as
−∞, may be employed for entries associated with absent arcs.
Given vectors ξ(u1 ) and ξ(u2 ), the conditions of Theorem 1 can be verified in time O(K|Li+1 |). If
the verification above is performed for each pair of nodes composing layer Li , all merging operations
26

can be verified in time O(K|Li |2 |Li+1 |). As each pair of nodes in the network is compared at most
once, we have a total running time of O(K|L|3 ). By sorting nodes in Li using the lexicographic
order defined by vectors ξ(u), comparisons can be restricted to neighbors defined by the resulting
sequence. In this case, only O(|L| log(|L|)) comparisons are necessary and the running time of the
whole procedure is dominated by the sorting operation, which consumes time O(K|L|2 log(|L|)) for
the whole network.

A.3

Proof of Theorem 2

We consider a reduction from the following NP-Hard problem, written in terms of our notation and
objective sense:
Unconstrained Biobjective (Ehrgott
Proposition
Given c1 , c2 ∈ Zn and d1 , d2 ∈
Pn 2006b,
Pn 8.12):
n
1
2
Z, does there exist x ∈ B such that j=1 cj xj ≥ d1 and j=1 cj xj ≥ d2 ?
Let c1 , c2 , and d := (d1 , d2 ) define the parameters of an instance of the unconstrained biobjective
problem. Equivalently, the problem asks if the given vector d is dominated by a point in the image
space of the MODO defined as
n
o
max (c1 )> x, (c2 )> x : x ∈ Bn .
(UCB)
We construct a network Ñ where deciding if d ∈ YN (Ñ ) is equivalent to applying an arc-removal
VPO for an appropriate subset of the arcs of Ñ . To this end, let us first consider a valid network
model N = (L, A) for UCB. Since there are no constraints, we can construct a straightforward
valid network model with one node per layer Lj , j ∈ [n + 1], and two arcs a1 , a2 connecting the
node u1 from layer Lj to node u2 in layer
 Lj+1 , j ∈ [n]. In particular, w (a1 ) = 0 (representing the
assignment xj = 0) and w (a2 ) = c1j , c2j (representing the assignment xj = 1). The network N is
depicted in Figure 14(a). Note that it has size polynomial in n, consisting of n + 1 nodes and 2n
arcs.
rN

(c11 , c12 )

0

rN 0
0

d = w(a00 )

w(a ) = 0
rN

0
tN

(c1n , c2n )

0

0

0
tN 0

tN

(b) N 0

(a) N

Figure 14: A network model N of a binary unconstrained MODO (a), and an extended network
model Ñ of N (b), used in the proof of Theorem 2.
We now design a new network model Ñ by extending N as depicted in Figure 14(b). The root
node rÑ is the arc-root of two arcs a0 and a00 . Arc a0 has arc-weight w(a0 ) = 0 and connects rÑ
27

to the corresponding rN in the new network. Arc a00 has arc-weight w(a00 ) = d, and Ñ contains a
single path p00 = (rÑ , t (a00 ) , . . . , tÑ ) connecting rÑ to tÑ that traverses node t (a00 ); the arc-weight
of all arcs in p00 (except a00 ) is equal to zero. Finally, a single arc with arc-weight 0 connects tÑ to
the associated tN in the new node.
Let A0 be the subset of AÑ containing all arcs that compose path p00 . If d is dominated
0
by YN 
(Ñ ), any
i.e.,

 frontier,
 arc ain A may be removed from Ñ without changing its Pareto
PF Ñ = PF Ñ − a . Otherwise, by construction, d necessarily belongs to PF Ñ , thus d ∈
PF (N ). Hence, no arc is removed from A0 . Since all steps are polynomially bounded on n, the
result follows.

A.4

Proof of Theorem 3

Suppose by contradiction that there exists a pair of nodes u and v that are isolating in N , together
with an arc a whose removal is a VPO for N [u, v] but not for N . This implies in the existence of
a path p in PN from rN to tN containing a for which w(p) ∈ PF (N ) and such that there exists
no other path in PN with the same path-weight. Since a ∈ AN [u,v] , and u and v are isolating, p
must traverse u and v. Let p0 be the subpath of p directed from u to v. Additionally, let p1 be the
portion of p from rN to u, and let p2 be the portion of p from v to tN .
Since the removal of a from N [u, v] does not alter PF (N [u, v]), there must exist a path p00 from
u to v for which w(p0 ) ≺ w(p00 ). Let p̃ be the path in N constructed from concatenating p1 , p00 , and
p2 . Then, w(p̃)  w(p) (as w(p̃) 6= w(p) by assumption), which contradicts that w(p) ∈ PF (N ).

B

Recursive Model Example
X = {x ∈ B7 : x1 + x2 + x3 ≤ 1,
x2 + x3 + x4 ≤ 1,
x4 + x5 ≤ 1,
x4 + x6 ≤ 1,
x5 + x7 ≤ 1,
x6 + x7 ≤ 1}.

28

s0 = (0, 0, 0, 0, 0, 0)
x1 = 0

x1 = 1

u2
1 = (0, 0, 0, 0, 0, 0)

x2 = 1

u2
2 = (1, 0, 0, 0, 0, 0)

x2 = 0

u3
1 = (1, 1, 0, 0, 0, 0)

x2 = 0

u3
2 = (0, 0, 0, 0, 0, 0)

x3 = 0

x3 = 1

x3 = 0

u4
1 = (1, 1, 0, 0, 0, 0)

x4 = 0

x4 = 0

x4 = 1

u5
2 = (1, 1, 1, 1, 0, 0)

x5 = 1

x5 = 0

u6
1 = (1, 1, 1, 0, 1, 0)

u7
1 = (1, 1, 1, 1, 1, 1)

x5 = 0

u6
2 = (1, 1, 1, 0, 0, 0)

x6 = 0

x6 = 1

u7
2 = (1, 1, 1, 1, 1, 0)

x7 = 0

x3 = 0

u4
2 = (1, 0, 0, 0, 0, 0)

u5
1 = (1, 1, 0, 0, 0, 0)

x6 = 1

u3
3 = (1, 0, 0, 0, 0, 0)

x7 = 0

u7
3 = (1, 1, 1, 1, 0, 1)

x7 = 0

u6
3 = (1, 1, 1, 1, 0, 0)

x6 = 0

x6 = 0

u7
4 = (1, 1, 1, 1, 0, 0)

x7 = 0

x7 = 1

st = (1, 1, 1, 1, 1, 1)

Figure 15: A valid network model for M in Example 1 obtained by recursive formulation given
in Example 3.

C

Detailed Computational Results

In this section we provide a comprehensive reporting of the computational results. Each line in each
table presents the aggregated results of instances whose dimensions are described in the leftmost
column: n indicates the number of variables and K gives the number of objective functions. The
other columns are arranged in groups of three and present the results for the algorithm indicated on
the top. Column S gives the number of instances that were solved within the predefined memory
and time limits, and t̄ is the average runtime of the algorithm to solve these instances. Column
(M, T) report pairs (m, t), indicating that the algorithm could not solve m and t instances of that
particular configuration because it exceeded the predefined memory and time limits, respectively.

29

30

K

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

3
4
5
6
7

n

20
20
20
20
20

30
30
30
30
30

40
40
40
40
40

50
50
50
50
50

60
60
60
60
60

70
70
70
70
70

80
80
80
80
80

90
90
90
90
90

100
100
100
100
100

8
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
4
0
0
0

10
10
0
0
0

10
10
6
0
0

S

K

1631.27
-

778.15
-

401.49
-

271.77
-

107.96
-

36.97
-

22.07
1685.96
-

5.22
676.25
-

0.93
34.41
200.79
-

t̄

(0,2)
(0,10)
(2,8)
(2,8)
(0,10)

(0,0)
(0,10)
(3,7)
(0,10)
(0,10)

(0,0)
(0,10)
(3,7)
(0,10)
(0,10)

(0,0)
(0,10)
(1,9)
(0,10)
(1,9)

(0,0)
(0,10)
(1,9)
(1,9)
(0,10)

(0,0)
(0,10)
(0,10)
(1,9)
(0,10)

(0,0)
(0,6)
(1,9)
(0,10)
(0,10)

(0,0)
(0,0)
(1,9)
(0,10)
(0,10)

(0,0)
(0,0)
(0,4)
(2,8)
(0,10)

(M, T)

8
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
0
0
0
0

10
2
0
0
0

10
5
0
0
0

10
10
0
0
0

10
10
8
0
0

S

O

1530.55
-

850.51
-

489.45
-

340.95
-

143.68
-

52.2
2841.42
-

30.81
1988.99
-

6.98
452.18
-

1.26
61.45
766.91
-

t̄

(0,2)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,8)
(0,10)
(0,10)
(0,10)

(0,0)
(0,5)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,2)
(0,10)
(0,10)

(M, T)

10
0
0
0
0

10
1
0
0
0

10
4
0
0
0

10
10
2
0
0

10
10
8
6
1

10
10
10
9
9

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

TD

967.47
-

283.73
2421.44
-

95.62
846.36
-

28.85
1092.52
3419.18
-

9.04
169.13
1439.71
2722.71
2762.75

2.25
23.46
110.5
484.95
1737.94

0.67
3.22
11.66
37.16
89.58

0.14
0.18
0.3
0.72
1.75

0.02
0.03
0.03
0.03
0.04

t̄

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,9)
(0,10)
(0,10)
(0,10)

(0,0)
(0,6)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,8)
(0,10)
(0,10)

(0,0)
(0,0)
(0,2)
(0,4)
(0,9)

(0,0)
(0,0)
(0,0)
(0,1)
(0,1)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

7
0
0
0
0

10
0
0
0
0

10
1
0
0
0

10
4
0
0
0

10
10
2
0
0

10
10
10
5
2

10
10
10
10
9

10
10
10
10
10

10
10
10
10
10

S

BU

1684.11
-

748.69
-

356
2285.8
-

125.02
1605.46
-

25.51
792.27
2158.86
-

5.17
134.57
1122.29
1305.47
1654.82

1.69
20.57
46.41
251.22
281.21

0.22
0.53
1.42
4.78
6.24

0.02
0.03
0.03
0.04
0.07

t̄

(0,3)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,9)
(0,10)
(0,10)
(0,10)

(0,0)
(0,6)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,8)
(0,10)
(0,10)

(0,0)
(0,0)
(0,0)
(0,5)
(0,8)

(0,0)
(0,0)
(0,0)
(0,0)
(0,1)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

10
1
0
0
0

10
7
0
0
0

10
10
4
0
0

10
10
6
2
1

10
10
10
6
3

10
10
10
10
9

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

262.53
745.24
-

71.74
1597.38
-

29.46
997.93
2676.08
-

11.85
257.09
2052.85
1998.62
3338.91

4.6
51.03
601.92
1389.04
2454.84

1.71
11.04
64.76
507.23
1195.03

0.64
3.16
8.05
39.25
86.99

0.16
0.2
0.34
0.87
1.63

0.02
0.03
0.03
0.03
0.04

t̄

Coup

(0,0)
(0,9)
(0,10)
(0,10)
(0,10)

(0,0)
(0,3)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,6)
(0,10)
(0,10)

(0,0)
(0,0)
(0,4)
(0,8)
(0,9)

(0,0)
(0,0)
(0,0)
(0,4)
(0,7)

(0,0)
(0,0)
(0,0)
(0,0)
(0,1)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

Table 1: Knapsack results

10
0
0
0
0

10
1
0
0
0

10
4
0
0
0

10
9
0
0
0

10
10
6
2
0

10
10
10
9
6

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

TD+

1048.7
-

289.5
2720.16
-

98.89
912.5
-

28.43
903.07
-

9.07
191.15
1007.89
1900.28
-

2.17
28.59
136.59
724.16
1514.52

0.66
3.79
16.54
55.53
130.27

0.13
0.18
0.38
0.97
2.8

0.02
0.02
0.03
0.03
0.05

t̄

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

(0,0)
(0,9)
(0,10)
(0,10)
(0,10)

(0,0)
(0,6)
(0,10)
(0,10)
(0,10)

(0,0)
(0,1)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,4)
(0,8)
(0,10)

(0,0)
(0,0)
(0,0)
(0,1)
(0,4)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

7
0
0
0
0

10
0
0
0
0

10
1
0
0
0

10
4
0
0
0

10
10
1
0
0

10
10
10
5
2

10
10
10
10
9

10
10
10
10
10

10
10
10
10
10

S

BU+

1688.68
-

750.2
-

351.56
2396.79
-

126.12
1556.9
-

23.45
784.7
834.29
-

5.33
125.92
996.69
1249.05
1716.35

1.64
20.42
40.95
240.71
250.89

0.21
0.54
1.27
4.22
5.25

0.02
0.03
0.03
0.03
0.05

t̄

(0,3)
(0,10)
(0,10)
(0,10)
(3,7)

(0,0)
(0,10)
(0,10)
(1,9)
(0,10)

(0,0)
(0,9)
(0,10)
(1,9)
(1,9)

(0,0)
(0,6)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,9)
(0,10)
(0,10)

(0,0)
(0,0)
(0,0)
(0,5)
(0,8)

(0,0)
(0,0)
(0,0)
(0,0)
(0,1)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

10
2
0
0
0

10
8
0
0
0

10
10
4
0
0

10
10
8
2
2

10
10
10
7
5

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

212.64
1764.28
-

59.84
1499.22
-

23.84
1012.99
2225.35
-

9.08
163.71
1554.53
1165.76
2254.76

3.97
34.05
416.31
1124.37
2174.34

1.52
6.74
38.8
400.63
971.73

0.52
1.88
4.44
20.83
63.02

0.13
0.18
0.26
0.5
0.97

0.02
0.03
0.03
0.02
0.04

t̄

Coup+

(0,0)
(0,8)
(0,10)
(0,10)
(0,10)

(0,0)
(0,2)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,6)
(0,10)
(0,10)

(0,0)
(0,0)
(0,2)
(0,8)
(0,8)

(0,0)
(0,0)
(0,0)
(0,3)
(0,5)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

Table 2: Set covering results
K
n

K

S

100
100
100
100
100

3
4
5
6
7

10
10
1
0
0

150
150
150
150
150

3
4
5
6
7

200
200
200
200
200

3
4
5
6
7

t̄

O
(M, T)

S

8.33
326.36
105.35
-

(0,0)
(0,0)
(0,9)
(0,10)
(0,10)

10
10
1
0
0

10
7
0
0
0

34.6
1290.86
-

(0,0)
(0,3)
(0,10)
(0,10)
(0,10)

10
0
0
0
0

97.47
-

(0,0)
(0,10)
(2,8)
(0,10)
(0,10)

t̄

TD
(M, T)

S

t̄

10.81
427.21
553.49
-

(0,0)
(0,0)
(0,9)
(0,10)
(0,10)

10
10
10
10
10

10
9
0
0
0

44.22
1485.95
-

(0,0)
(0,1)
(0,10)
(0,10)
(0,10)

10
0
0
0
0

124.17
-

(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

BU
(M, T)

S

10.14
10.7
11.27
11.91
24.09

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

10
10
10
10
10

10
10
7
8
4

36.82
83.76
196.1
539.84
783.5

(0,0)
(0,0)
(0,3)
(0,2)
(0,6)

1
1
0
0
0

63.42
292.97
-

(0,9)
(0,9)
(0,10)
(0,10)
(0,10)

t̄

Coup
(M, T)

S

t̄

10.34
11.32
15.26
35.35
54.9

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

10
10
10
10
10

10.18
10.61
10.82
11.17
16.01

(M, T)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

10
8
5
5
3

59.11
197.08
652.01
803.55
1392.27

(0,0)
(0,2)
(0,5)
(0,5)
(0,7)

10
9
7
9
4

36.27
62.88
136.13
315.13
454.3

(0,0)
(0,1)
(0,3)
(0,1)
(1,5)

1
0
0
0
0

202.83
-

(0,9)
(0,10)
(0,10)
(0,10)
(0,10)

1
0
0
0
0

158.17
-

(0,9)
(0,10)
(0,10)
(0,10)
(0,10)

(M, T)

S

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,5)
(0,5)
(0,8)
(0,10)
(0,10)
(0,10)
(0,10)
(0,10)

10
10
10
10
10
10
10
10
10
10
1
1
0
0
0

Table 3: Set packing results
K
n

K

S

100
100
100
100
100
150
150
150
150
150
200
200
200
200
200

3
4
5
6
7
3
4
5
6
7
3
4
5
6
7

10
10
0
0
0
10
1
0
0
0
10
0
0
0
0

t̄
9.38
518.46
23.21
2916.26
76.47
-

O
(M, T)

S

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)
(0,0)
(0,9)
(2,8)
(1,9)
(0,10)
(0,0)
(0,10)
(0,10)
(1,9)
(0,10)

10
10
0
0
0
10
2
0
0
0
10
0
0
0
0

t̄
14.34
954.85
36.33
3057.22
111.18
-

TD
(M, T)

S

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)
(0,0)
(0,8)
(0,10)
(0,10)
(0,10)
(0,0)
(0,10)
(0,10)
(0,10)
(0,10)

10
10
10
10
10
10
10
10
10
10
1
0
0
0
0

t̄
0.09
0.15
0.28
0.77
3.6
12.11
27.97
87.74
203.97
888.75
121.92
-

31

BU
(M, T)

S

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,9)
(0,10)
(0,10)
(0,10)
(0,10)

10
10
10
10
10
10
10
5
5
2
0
0
0
0
0

t̄
0.12
0.23
0.56
0.86
3.73
17.94
75.58
155.74
318.31
196.25
-

Coup
t̄
0.08
0.13
0.24
0.57
1.28
12.41
23.16
65.68
125.62
475.53
137.78
336.82
-

(M, T)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,0)
(0,9)
(0,9)
(0,10)
(0,10)
(0,10)

32

K
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7
3
4
5
6
7

n
15
15
15
15
15
15
15
15
15
15
15
15
15
15
15
20
20
20
20
20
20
20
20
20
20
20
20
20
20
20
25
25
25
25
25
25
25
25
25
25
25
25
25
25
25

50
50
50
50
50

40
40
40
40
40

30
30
30
30
30

50
50
50
50
50

40
40
40
40
40

30
30
30
30
30

50
50
50
50
50

40
40
40
40
40

30
30
30
30
30

C

10
6
0
0
0

10
8
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
4
0
0

S

K

87.75
2144.1
-

86.15
2974
-

57.61
1209.4
-

42.55
933.77
-

32.85
553.77
-

7.86
89.82
-

2.76
31.58
-

2.05
23.02
-

1.21
8.31
1882.3
-

t̄

(0,0)
(0,4)
(0,10)
(0,10)
(0,10)

(0,0)
(0,2)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(1,9)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(2,8)
(0,10)

(0,0)
(0,0)
(1,9)
(0,10)
(1,9)

(0,0)
(0,0)
(1,9)
(2,8)
(1,9)

(0,0)
(0,0)
(1,9)
(1,9)
(0,10)

(0,0)
(0,0)
(0,10)
(1,9)
(0,10)

(0,0)
(0,0)
(0,6)
(1,9)
(0,10)

(M, T)

10
2
0
0
0

10
2
0
0
0

10
9
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
0
0
0

10
10
9
0
0

S

O

86.78
2135.3
-

89.88
3461.5
-

69.52
2037.7
-

55.92
1835.4
-

40.07
1046.7
-

9.84
222.09
-

3.71
87.62
-

2.85
65.87
-

1.71
25.9
2166.9
-

t̄

(0,0)
(0,8)
(0,10)
(0,10)
(0,10)

(0,0)
(0,8)
(0,10)
(0,10)
(0,10)

(0,0)
(0,1)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,10)
(0,10)
(0,10)

(0,0)
(0,0)
(0,1)
(0,10)
(0,10)

(M, T)

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

TD

49.21
54.91
67.79
137.12
427.1

24.76
26.58
32.49
60.98
148.45

3.12
3.35
3.78
5.69
10.54

1.37
1.47
1.84
3.29
7.08

0.63
0.7
0.8
1.28
2.62

0.14
0.15
0.17
0.24
0.4

0.01
0.01
0.02
0.03
0.07

0.01
0.01
0.01
0.02
0.04

0
0
0
0
0

t̄

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

Table 4: Absolute value results
BU

53.73
60.58
72.88
111.41
222.51

25.68
27.51
33.66
51.1
95.4

3.26
3.59
3.98
5.37
8.97

1.56
1.69
2.15
3.5
7.03

0.74
0.76
0.9
1.43
2.67

0.14
0.14
0.19
0.28
0.49

0.01
0.02
0.02
0.04
0.09

0.01
0.01
0.01
0.02
0.04

0
0
0
0
0

t̄

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

10
10
10
10
10

S

Coup

99.72
106.52
128.78
171.18
358.04

41.76
43.84
51.75
77.37
144.34

4.21
4.43
5.05
7.08
11.03

1.98
2.47
2.78
3.95
6.26

0.94
1.06
1.09
1.56
2.41

0.17
0.18
0.21
0.29
0.38

0.02
0.02
0.02
0.03
0.05

0.01
0.01
0.01
0.02
0.04

0
0
0
0
0

t̄

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(0,0)
(0,0)
(0,0)
(0,0)
(0,0)

(M, T)

