A LP Relaxation based Matheuristic for Multi-objective Integer
Programming
Duleabom An
1
, Sophie N. Parragh
1 a
, Markus Sinnl
1 b
and Fabien Tricoire
2 c
1
Institute of Production and Logistics Management, Johannes Kepler University Linz,
Altenberger Straße 69, 4040 Linz, Austria
2
Institute for Transport and Logistics Management, Vienna University of Economics and Business,
Welthandelsplatz 1, 1020 Vienna, Austria
Keywords:
Three-objective Binary Integer Programming, Matheuristics, Path Relinking, Multi-objective Knapsack
Problem.
Abstract:
Motivated by the success of matheuristics in the single-objective domain, we propose a very simple linear
programming-based matheuristic for three-objective binary integer programming. To tackle the problem, we
obtain lower bound sets by means of the vector linear programming solver Bensolve. Then, simple heuristic
approaches, such as rounding and path relinking, are applied to this lower bound set to obtain high-quality
approximations of the optimal set of trade-off solutions. The proposed algorithm is compared to a recently
suggested algorithm which is, to the best of our knowledge, the only existing matheuristic method for three-
objective integer programming. Computational experiments show that our method produces a better approxi-
mation of the true Pareto front using significantly less time than the benchmark method on standard benchmark
instances for the three-objective knapsack problem.
1 INTRODUCTION
Many real-world optimisation problems involve mul-
tiple conflicting objectives, concerning, e.g., costs,
environmental impact or service level, and can be for-
mulated as integer linear programs. Kolli and Evans
(1999), for example, deal with facility location prob-
lems for a franchise company. When a new franchise
is launched, the franchisee and the franchisor have
conflicting objectives. To find the optimal location
for the new store, the model maximises the number
of customers while minimising conflict among exist-
ing franchises. Sawik et al. (2016) investigate vehi-
cle routing problems occurring in the logistics of the
food industry. The introduced model minimises both
the total travel distance and CO
2
emissions. A home
care routing and scheduling problem is investigated
by Braekers et al. (2016). Given a group of nurses and
a number of care tasks at the patients’ home locations,
the model assigns tasks to each nurse while minimis-
ing operating cost and client inconvenience with re-
spect to the timing of the visit and the nurse per-
a
https://orcid.org/0000-0002-7428-9770
b
https://orcid.org/0000-0003-1439-8702
c
https://orcid.org/0000-0002-3700-5134
forming the task. Kovacs et al. (2015) study a multi-
objective consistent vehicle routing problem. To pro-
vide consistent service in the routing industry, they
maximise driver consistency and arrival time consis-
tency while minimising total routing costs.
The main goal in multi-objective (MO) optimi-
sation is to generate the set of optimal trade-off so-
lutions, known as efficient (or Pareto optimal) solu-
tions. Our study focuses on binary integer program-
ming (IP) problems with three objectives.
Motivated by the success of matheuristics in the
single objective domain, we present a very sim-
ple matheuristic for three-objective integer program-
ming. To the best of our knowledge, only one other
matheuristic for three-objective integer programming
has been developed so far (Pal and Charkhgard,
2019b). The algorithm we propose relies on a lower
bound set which is obtained from the linear pro-
gramming (LP) relaxation, using the vector linear
programming solver Bensolve (L
¨
ohne and Weißing,
2017). The lower bound set is defined by its ex-
treme points (and edges). Starting from those solu-
tions which give rise to the extreme points, we apply
rounding in combination with path relinking to obtain
high-quality approximations of the true Pareto fron-
88
An, D., Parragh, S., Sinnl, M. and Tricoire, F.
A LP Relaxation based Matheuristic for Multi-objective Integer Programming.
DOI: 10.5220/0010347000880098
In Proceedings of the 10th International Conference on Operations Research and Enterprise Systems (ICORES 2021), pages 88-98
ISBN: 978-989-758-485-5; ISSN: 2184-4372
Copyright
c
2023 by SCITEPRESS Science and Technology Publications, Lda. Under CC license (CC BY-NC-ND 4.0)
tier (PF).
The contribution of this paper is twofold:
We show that, in the case of the three-objective
assignment problem, since the extreme solutions
which Bensolve produces are integer, it already
provides high-quality approximations of the true
Pareto frontier, even without additional ingredi-
ents.
We propose the first LP relaxation-based
matheuristic algorithm for three-objective integer
programming, combining Bensolve with rounding
and path relinking (PR)
The remainder of the paper is structured as fol-
lows. Section 2 briefly reviews existing work in multi-
objective integer programming (MOIP) matheuristics.
Section 3 provides basic concepts and background
for MOIP. The proposed algorithms are described in
Section 4 and empirically evaluated on the multi-
objective assignment problem (MOAP) and multi-
objective knapsack problem (MOKP) in Section 5.
Finally, the paper concludes with a summary and sug-
gestions for future work in Section 6.
2 RELATED WORK
Over the past years a number of generic exact meth-
ods for solving MOIP have been proposed (see, e.g.,
Mavrotas and Florios (2013), Zhang and Reimann
(2014), Kirlik and Sayın (2014), Boland et al. (2017)).
In spite of their popularity in the single objec-
tive domain, only comparably few contributions on
matheuristc methods exist.
A heuristic and a matheuristic approach for bi-
objective mixed binary integer linear programming
are proposed by Soylu (2015). The methods are
variants of variable neighbourhood search and local
branching. Both algorithms collect segments of the
PF during the search and combine them at the end. To
deal with bi-objective binary IP, Leitner et al. (2016)
suggest an exact method and a matheuristic frame-
work. The general idea is inspired by the fact that ef-
ficient solutions in the same area often share common
features. In each phase, their heuristic obtains feasi-
ble solutions by fixing a large number of variables and
reducing the associated feasible region, respectively.
Based on the feasibility pump (FP) idea introduced
by Fischetti et al. (2005), Pal and Charkhgard (2019a)
suggested a FP based heuristic for bi-objective IP and
extended the method to higher dimensions in (Pal and
Charkhgard, 2019b). The proposed method works in
two stages. Both stages employ a FP based heuristic.
In the first stage, where the purpose is to generate well
spread solutions, fractional solutions are computed by
means of a weighted sum method. The purpose of the
second stage is to generate additional solutions. It re-
lies on local branching in combination with the devel-
oped FP based heuristic. To the best of our knowl-
edge, this is the only matheuristic method that deals
with MOIP with more than two objectives. Therefore,
we use it as a benchmark and refer to it as FPBH.
3 BASIC CONCEPTS AND
BACKGROUND
The problem we consider is a MOIP, with binary in-
teger decision variables and three objectives. In the
following, we state the MOIP in its general form. We
present a unified view using minimisation objectives.
y = min{C(x) : Ax b,x {0,1}
n
},
(MOIP)
where x
j
, j = 1,2, . . .,n, is the vector of decision vari-
ables and X := {Ax b, x {0,1}
n
} is the feasible
set. C is a p × n objective function matrix where c
k
,
(k = 1,2,..., p), is the k
th
row of C. In our case, p = 3.
Further, y is a set of points in the objective space (cri-
terion space), each of which corresponds to at least
one solution vector x X. A is an m × n constraint
matrix and b is the vector of right-hand-side values
for these constraints.
3.1 Pareto Dominance
In MO optimisation, the quality of a solution is de-
termined by Pareto dominance. Suppose there are
two solutions x and x
of Problem (MOIP). Then,
x dominates x
if and only if c
k
(x) c
k
(x
) for all
k {1,..., p} and c
k
(x) < c
k
(x
) for at least one k.
If there does not exist any x
that dominates x, then x
is Pareto optimal. If x
is a Pareto optimal (efficient)
solution, then y
is a non-dominated point. The set of
all non-dominated points is called the Pareto front.
3.2 Weighted Sum Method and
Supported Solutions
The weighted sum method is a commonly used ap-
proach in solving MO optimisation problems. It com-
bines multiple objectives into one as follows:
F(x) = w
1
f
1
(x) + w
2
f
2
(x) + · · ·+ w
p
f
p
(x)
p
i=1
w
i
= 1, w
i
(0,1),
where w
i
is the weight of objective function i. If a
solution x
can be found by the weighted sum method
A LP Relaxation based Matheuristic for Multi-objective Integer Programming
89
(i.e. optimising a convex combination of all objective
functions), it is called a supported efficient solution.
Otherwise, it is a non-supported efficient solution.
3.3 LP Relaxation and Bound Set
The notion of bound set was introduced by Ehrgott
and Gandibleux (2007). In the context of Problem
(MOIP), a lower bound (LB) set and an upper bound
(UB) set provide information about the variable range
that efficient solutions of the MOIP can attain. One
common way to obtain a LB set for a minimisation
problem is to solve the LP relaxation of the original
problem. For constructing the LP relaxation of Prob-
lem (MOIP), the integrality conditions are removed,
i.e. 0 x
j
1, j = 1,. . . ,n and Bensolve can then
be used to compute the LB set. In our heuristics, we
will use the solutions associated with the LB set. In
a slight abuse of notation, we will refer to these solu-
tions as LB set.
3.4 Benchmark Problems
The multi-objective assignment problem (MOAP)
and the multi-objective knapsack problem (MOKP)
are commonly used benchmark problems in MO. We
use standard benchmark instances of these two prob-
lems for our computational experiments.
3.4.1 Multi-objective Assignment Problem
The well-known assignment problem is a type of
transportation problem: a certain number of tasks
l {1, . . .,n} and agents r {1,..., n} are given.
The decision variable x
rl
denotes whether task l is as-
signed to the agent r (x
rl
=1) or not (x
rl
=0). When a
task is allocated to an agent, the corresponding non-
negative costs c
1
rl
,...,c
p
rl
are incurred. The MOAP
can be stated as follows:
min
n
r=1
n
l=1
c
j
rl
x
rl
j = 1,..., p
(1)
s.t.
n
l=1
x
rl
= 1 r = 1, . . .,n
(2)
n
r=1
x
rl
= 1 l = 1,...,n
(3)
x
rl
{0,1} r,l = 1, ...,n.
(4)
The objective of the MOAP is to find an optimal as-
signment of all tasks to agents while minimising p
cost functions (1). Equation (2) ensures that each
agent is assigned to only one task. Equation (3) limits
each task to be assigned to one agent only.
It is well-known that the constraint matrix of AP
is totally unimodular, then every vertex of the LP re-
laxation is an integer vector. Thus, by solving the LP
relaxation we can naturally obtain integer solutions of
the AP.
Although Bensolve only computes supported ef-
ficient solutions, these solutions already provide a
high-quality approximation of PF. The correspond-
ing computational results are given in Table 2 in Sec-
tion 5. In conclusion, the MOAP may not be a suitable
benchmark for MOIP heuristics.
3.4.2 Multi-objective Knapsack Problem
In the multi-objective knapsack problem, a set of
items is given, each with a certain weight w
r
, and we
must select a subset of these items such that the total
weight does not exceed a given capacity W . The de-
cision variable x
r
denotes whether item r is selected
for the knapsack (x
r
=1) or not (x
r
=0). Each item r
also has profits v
1
r
...v
p
r
. Here, W, v
r
, and w
r
are non-
negative integer values. The MOKP model is stated
as follows:
max
n
r=1
v
j
r
x
r
j = 1,..., p
(5)
s.t.
n
r=1
w
r
x
r
W
(6)
x
r
{0,1} r = 1,...,n
(7)
Equation (5) denotes the objective functions maximis-
ing the p total profits of selected items. Equation (6) is
the capacity constraint. The total weight of the items
placed in the knapsack cannot exceed the given ca-
pacity.
In this paper, we convert a maximisation objective
function into a minimisation one by multiplying it
by 1.
4 LP RELAXATION-BASED
MATHEURISTIC
This section provides the overview of the proposed
algorithm and its main ingredients.
4.1 Algorithmic Framework
The proposed matheuristic algorithm follows a two-
stage approach. At the first stage, we obtain LB
sets and adapt them to feasible integer rounded
sets. A variant of the weighted sum method by
¨
Ozpeynirci and K
¨
oksalan (2010), Bensolve by L
¨
ohne
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
90
and Weißing (2017) and Inner approximation by Csir-
maz (2020) have been investigated to find LB sets
first. Among the three methods, Bensolve shows
competitive performance in the experiment and pro-
duces all the required bound set information we need.
The performance comparison can be found in the Ap-
pendix. In the context of MOKP, once LB sets (L)
are obtained, we round down the fractional variable
to produce a feasible solution. These are referred to
as integer rounded (IR) sets. At the second stage, PR
is employed to improve the solution quality. Until
reaching the iteration limit, the PR process repeats.
If PR finds a new solution, it is stored in the archive
candX . Based on initial experiments, we set the limit
of iterations of PR to the number of solutions of the
IR set times 50. Since the dominance relation is not
checked in
˜
X during the search, dominated solutions
are filtered after the algorithm terminates, and the set
of integer feasible solutions X is returned. Algorithm
1 describes the general framework of the proposed al-
gorithm.
Algorithm 1: The LP relaxation-based matheuristic
framework.
Input: L: points describing an LB set
1 candX: an archive of newly found feasible IP
solutions by PR
2
˜
X: an empty list
3 IR RoundingDown(L);
4 i = 0;
5 for i< iteration limit do
6 candX PathRelinking(IR,S
i
,S
g
,IGPair);
7
˜
X
˜
X candX;
8 i = i + 1;
9 X DominanceCheck(IR
˜
X);
Output: X
4.2 Path Relinking
Although IR sets obtained from rounding are feasi-
ble, they are not necessarily of high-quality. There-
fore, we employ PR to increase the number of solu-
tions and improve the quality of the approximate PF.
PR was originally introduced by Glover (1997). The
main idea of it is that there should be common char-
acteristics among high-quality solutions. The method
produces new solutions by exploring solution ”paths”
between pairs of known solutions. To generate a new
solution, PR chooses two solutions from a set of ini-
tial solutions; an initiating solution (S
i
) and a guiding
solution (S
g
) represent the starting and ending points
of the path, respectively. Then it explores the trajecto-
ries that link the pair of solutions in a neighbourhood
space. In each step, a certain number of intermedi-
ate paths, called a neighbourhood, are created. For
the next step, the new initiating solution is selected in
there. As the search proceeds, more attributes of the
guiding solutions are gradually passed on to interme-
diate solutions. The search continues until the initiat-
ing solution becomes identical to the guiding solution.
The method has been applied successfully to MO
problems such as the travelling salesman problem
(Jaszkiewicz, 2005), the dial-a-ride problem (Par-
ragh et al., 2009), and a school bus routing problem
(de Souza Lima et al., 2017). It also produces high-
quality solutions on large MOKP instances within a
reasonable timeframe (Beausoleil et al., 2008; Mart
´
ı
et al., 2015). For these reasons, we combine PR with
the first stage approach.
For illustration purposes, we provide a short ex-
ample describing the PR process on a MOKP with
three objectives. Suppose we have a four-item MOKP
in which the initiating solution is (0 0 1 0) and the
guiding solution is (1 1 0 0). In this case, the only se-
lected item in common is the fourth item. Therefore,
three intermediate paths are created by the following
rules. To create one path, one variable value in the
initiating solution is switched. To be specific, if the
item is placed in the knapsack and its value is 1, it
is taken out of the knapsack and the value changes to
0. If the item is not chosen and the value equals 0,
it is placed in the knapsack and the value changes to
1. This procedure applies to all the different variable
values. The outcome of the process is seen in the first
row (Neighbourhood) of Table 1. When the P matrix
(8) represents the profits of the three objectives, the
sets of total profits that correspond to the first neigh-
bourhood are [7,6,8], [5,4,6], and [0,0,0]. Since the
dominance relation among the solution points is clear,
the first path (1 0 1 0) is selected as the new initiating
solution (marked with * in Table 1). This process re-
peats until the initiating solution becomes identical to
the guiding solution. Table 1 illustrates the transform-
ing process of PR in this example.
P =
4 2 3 6
5 3 1 8
6 4 2 7
(8)
Table 1: Path relinking procedure.
Initiating Guiding Neighbourhood
0 0 1 0 1 1 0 0 1 0 1 0* 0 1 1 0 0 0 0 0
1 0 1 0 1 1 0 0 1 1 1 0* 1 0 0 0
1 1 1 0 1 1 0 0 1 1 0 0*
1 1 0 0 1 1 0 0
A LP Relaxation based Matheuristic for Multi-objective Integer Programming
91
In our PR process, an infeasible intermediate so-
lution is not stored in the solution archive, but still
can be a new initiating solution. This is for letting the
algorithm explore a broader search region so that it
possibly finds diverse solutions.
Depending on the specific problem, the following
components of PR can be designed differently in the
algorithm.
1) Building an initial solution set
2) Selecting S
i
and S
g
3) Generating intermediate paths (a neighbourhood)
4) Choosing the next initiating solution
In this study:
1) We use IR sets as the initial solution set.
2) S
i
and S
g
are chosen either randomly or based on
similarity between them.
3) A new neighbourhood is generated as many times
as the number of different variable values.
4) The next S
i
is determined either by dominance re-
lation based analysis or randomly.
The description of the entire PR implementation
for the proposed algorithm is given in Algorithm 2.
The algorithm maintains two archives candX and
IGPair. Each archive stores newly found solutions
and used S
i
-S
g
pairs during the PR process. Since the
performance of PR can vary depending on the choice
of the S
i
-S
g
pair, we investigate different ways, and
call SelectionRuleI. The first approach is to select two
random solutions from IR. For the second and third
methods, S
i
is chosen at random in IR first, then the
similarity between S
i
and all the other solutions in IR
is calculated. The similarity is measured by count-
ing the number of variables which have equal values.
Thus, higher values imply a greater similarity. The
second method selects S
g
by finding the most simi-
lar solution to S
i
. The third approach finds the most
different solution to S
i
for S
g
. When S
i
and S
g
are sim-
ilar, fewer intermediate paths are likely to be created.
However, when they are different, the PR process can
take longer as it explores relatively diverse intermedi-
ate paths.
Once S
i
and S
g
are chosen by SelectionRuleI (line 4),
PR continues until one of the following two termina-
tion conditions is met (line 5):
The initial and guiding solutions are identical.
The pair of S
i
and S
g
has already been chosen.
To decide the number of intermediate paths, we need
to find the positions of the differing variable values in
S
i
and S
g
. They are stored in items as indices (lines
6-7). Once items is defined, we generate the neigh-
bourhood by the following rules: For each index in
items, the variable values in S
i
change one by one
(line 8). Let us suppose items 1 and 5 have a differ-
ent value, items = {1, 5}. If item 1 is selected in
S
i
(variable value=1), then it is taken out of it. If it
is out of the knapsack (variable value=0), then it is
added to it. The same process is conducted on item 5.
Once intermediate paths are built, we select the next
S
i
. In order to avoid local optima, the best neighbour
is not systematically selected. Rather, it is selected
with a certain probability, set to 0.7 after initial exper-
iments. Otherwise, another neighbour is selected ran-
domly. In the case where we do want to select the best
neighbour, we analyse the dominance relation among
intermediate solutions and select the best solution as
S
i
(line 11). For this, SelectionRuleII is introduced. If
one non-dominated solution exists in the neighbour-
hood, it becomes the next initiating solution. If there
are mutually non-dominating solutions, we check the
improved ratio of each solution point. The approach
for finding the most-improved point is explained in
detail in Section 4.2.1. For the case where we do not
want to select the best neighbour, the analysis is not
needed and simply one random solution is selected
from the neighbourhood (lines 12-13). Before mov-
ing on to the next iteration, the algorithm checks the
feasibility of S
i
(line 14). If S
i
is feasible and not in-
cluded in IR, it is stored in the archive candX (line
16). The infeasible S
i
is not archived, though, it is still
used for the next iterations as it might help to find ad-
ditional feasible solutions in yet unexplored parts of
the search region. IR is updated by adding S
i
(line
15). To prevent the case in which the same pair of ini-
tial and guiding solution is chosen, the current [S
i
,S
g
]
pair is archived after every PR iteration (line 16). The
algorithm returns the set of integer feasible solutions
candX .
4.2.1 ImprovedND Operation
The ImprovedND operation figures out which path
shows the biggest improvement compared to the cur-
rent solution S
i
. Algorithm 3 shows a precise descrip-
tion of the operation.
Let ND be a set of non-dominated intermediate
points and ob jS
i
be the objective values of S
i
. To
record the improvement of each non-dominated point,
we create two matrices and one list. The ratio table
stores the ratio of a non-dominated point to the cur-
rent point ob jS
i
, for each objective (lines 4-6). Af-
terwards, we assign a rank to each column of the ra-
tio table and enter rankings into the rank table (lines
7-10). The rankings of each non-dominated point are
added up (lines 11-12) and stored in ND degree. The
non-dominated path with the highest degree is set to
S
i
(line 13).
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
92
Algorithm 2: The framework of path relinking.
Input: IR, S
i
, S
g
, IGPair
1 IGPair: an archive of S
i
-S
g
pairs
2 candX
/
0;
3 IGPair
/
0;
4 Select S
i
,S
g
from IR following
SelectionRuleI;
5 while S
i
̸= S
g
and [S
i
,S
g
] / IGPair do
6 items index set of different variable
values;
7 n #indices in items;
8 Create n neighbourhood;
9 The best move analysis:
10 if rand() < 0.7 then
11 S
i
is chosen by SelectionRuleII
12 else
13 S
i
randomly choose one solution
from the neighbourhood;
14 Feasibility check of S
i
;
15 if S
i
is feasible and S
i
/ IR then
16 candX S
i
;
17 IR IR S
i
;
18 IGPair [S
i
,S
g
];
Output: candX
5 COMPUTATIONAL
EXPERIMENTS
We evaluate the performance of the different versions
of our heuristic method, focusing on the comparison
with FPBH (Pal and Charkhgard, 2019b). The vari-
ant with rounding down is named RD. The following
PR variants are tested. The first three PR variants ran-
domly choose the next initiating solution during the
iterations.
PRrand: This version randomly selects both the
initiating and guiding solutions from the initial so-
lution set.
PRsim: This version randomly chooses the initi-
ating solution, then finds the most similar solution
among the remaining solution set for the guiding
solution.
PRdif: This version finds the most different solu-
tion to the already chosen initiating solution for
the guiding solution.
The three other variants consider the improvement
among mutually non-dominating intermediate paths
to choose the next initiating solution.
PI: This version selects the most improved inter-
mediate solution as the next initiating solution.
Algorithm 3: ImprovedND.
Input: ob jS
i
, ND
1 ratio table: |ND| × p matrix
2 rank table: |ND| × p matrix
3 ND degree : list of size |ND|
4 for i=1,...,|ND| do
5 for j=1,...,p do
6 ratio table[i][ j] =
ND[i][ j]
ob jS
i
[ j]
7 for i=1,...,|ND| do
8 for j=1,...,p do
9 rank table[i][ j] = rank of ND[i][ j]
10 in j
th
column
11 for i=1,...,|ND| do
12 ND degree[i] sum(rank table i
th
column)
13 S
i
ND with the highest degree
Output: S
i
PIsim: This variant uses ImprovedND within
PRsim.
PIdif: This variant uses ImprovedND in PRdif.
The proposed algorithms use Bensolve by L
¨
ohne
and Weißing (2017) to obtain the bound sets. The
heuristic integration (rounding down and PR), is im-
plemented in Julia. For the benchmark algorithm, we
used the Julia implementation of FPBH (with the de-
fault setting) which is publicly available at https://
github.com/aritrasep/FPBH.jl. All experiments
of matheuristics are carried out on Intel® Core™ i5-
8250U CPU running at 1.60GHz with 16GB RAM.
The exact MOIP solver proposed by Kirlik and Sayın
(2014) (KS) is also used in the experiment to obtain
the true PF for comparison purpose. KS is run on
Quad-core X5570 Xeon CPUs @2.93GHz with 48GB
RAM. The KS results are for reference only (i.e. not
for benchmarking).
The test instances we use are the same ones on
which FPBH is tested. The instances were generated
by Kirlik and Sayın (2014) and are publicly available
at http://home.ku.edu.tr/
˜
moolibrary/. Each
problem class has 100 instances divided into 10 sub-
classes, each of which contains 10 instances. MOAP
instances are formed in the number of tasks (to be as-
signed) which varies from 5 to 50 in increments of 5.
MOKP instances are classified by the number of items
which varies from 10 to 100 in increments of 10.
A LP Relaxation based Matheuristic for Multi-objective Integer Programming
93
5.1 Performance Measure:
Hypervolume Indicator
One widely used indicator to measure the quality of
a solution set in MO optimisation is the hypervolume
(HV) indicator. HV measures the volume of the dom-
inated space of all the solutions contained in a solu-
tion set. To calculate the dominated space, a refer-
ence point must be used. Usually, a reference point is
the “worst possible” point in the objective space. In
this study, all the HV values are calculated with nor-
malised objective values.
Let Y
k
N
be a set of k
th
objective values of the true
PF and y R
p
be an arbitrary non-dominated point
obtained from a heuristic algorithm. Then, the nor-
malised values of the obtained point are:
y
k
min(Y
k
N
)
max(Y
k
N
) min(Y
k
N
)
k = 1,..., p.
As all non-dominated points are normalised, their
values exist in [0,1]. Therefore, the reference
point is (1,1,1). Higher HV values indicate a bet-
ter approximation. We used the publicly avail-
able HV computing program provided by Fon-
seca et al. (2006) at http://lopez-ibanez.eu/
hypervolume#intro to obtain HV values.
5.2 Results and Discussion
We report the following results of each algorithm: the
number for solutions (|Y|), CPU time (sec), HV value,
and HV as a percentage of the HV indicator value
for the exact non-dominated set as provided by Kirlik
and Sayın (2014). All the figures are average results
over 10 test instances. The figures of PR variants and
FPBH are averaged over 10 runs for each instance be-
cause they have random components. The results of
the experiments on MOAP and MOKP are reported in
Tables 2-5.
Bensolve already finds integer feasible solutions
for all the MOAP instances. In addition, it shows
better performance than FPBH in all the subclasses.
Therefore, we do not apply our matheuristic to MOAP
instances.
In general, both the number of solutions and com-
putation time increase as the size of instances be-
comes larger. The difference between FPBH and Ben-
solve is clearly noticeable in Table 2, which shows
that Bensolve outperforms FPBH in all the MOAP in-
stances. Furthermore, the difference becomes greater
as the size of the instances grows. For example, for
the instances with more than 15 tasks (n15), the
number of solutions of Bensolve is more than double
Table 2: Comparing algorithm performance on MOAP for
p=3, * indicates optimal values, best heuristic values are in
bold.
n |Y | CPUtime(sec) HV HV(%)
KS* FPBH Bensolve KS FPBH Bensolve KS FPBH Bensolve FPBH Bensolve
5 14.1 6.7 7.5 0.11 0.05 0.001 6.77 6.60 6.71 97.49 99.11
10 176.8 21.0 39.0 10.71 0.28 0.019 7.23 6.92 7.18 95.71 99.31
15 674.9 40.4 83.1 92.51 1.17 0.068 7.33 6.96 7.29 94.95 99.45
20 1860.5 62.5 161.3 359.07 2.04 0.222 7.40 6.99 7.37 94.46 99.59
25 3567.8 90.6 253.1 872.19 4.99 0.596 7.46 7.05 7.44 94.50 99.73
30 6181.3 140.0 379.4 1859.74 15.59 1.157 7.48 7.04 7.46 94.12 99.73
35 8972.3 163.2 501.4 3285.57 26.39 2.082 7.50 7.06 7.48 94.13 99.73
40 14679.7 242.9 699.1 6425.98 57.95 3.824 7.53 7.09 7.51 94.16 99.73
45 17702.2 238.6 838.0 9239.01 82.16 6.097 7.56 7.10 7.54 93.92 99.74
50 24916.8 337.5 1034.8 15814.82 119.52 9.739 7.58 7.12 7.56 93.93 99.74
that of FPBH. Further, solutions are found in signifi-
cantly less time. Hence, not only the quantity but also
the quality of the solutions of Bensolve are better than
FPBH. It reaches more than 99% of the maximum
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
94
Table 3: Comparing |Y| of algorithms on MOKP for p=3, * indicates optimal values, best heuristic values are in bold.
PR variants
n KS* FPBH RD PRrand PRsim PRdif PI PIsim PIdif
10 9.8 5.1 4.3 5.7 4.8 4.8 6.3 4.6 4.6
20 38.0 18.1 10.7 22.4 20.7 19.0 26.1 20.4 18.6
30 115.8 43.2 20.7 47.9 46.9 45.6 58.6 46.3 45.0
40 311.2 95.7 33.8 95.8 96.4 91.3 122.0 97.8 93.2
50 444.2 111.8 41.7 119.1 118.1 114.9 143.5 118.4 112.6
60 917.1 195.1 71.5 209.1 209.8 203.6 266.2 207.8 208.4
70 1643.4 348.2 90.2 263.0 264.4 259.8 353.4 271.3 264.9
80 2295.8 439.0 113.1 305.1 301.1 310.1 399.9 313.7 309.2
90 3107.8 501.9 130.6 322.7 319.4 327.0 412.7 338.7 343.3
100 5849.0 919.2 176.7 442.8 453.0 439.8 581.3 469.4 471.9
Table 4: Comparing CPU time (sec) of algorithms on MOKP for p=3, * indicates optimal values, best heuristic values are in
bold.
PR variants
n KS* FPBH RD PRrand PRsim PRdif PI PIsim PIdif
10 0.140 0.023 0.001 0.002 0.004 0.003 0.006 0.002 0.002
20 1.030 0.080 0.004 0.020 0.056 0.044 0.067 0.051 0.036
30 5.540 0.324 0.008 0.085 0.314 0.281 0.307 0.312 0.268
40 23.23 1.071 0.013 0.224 0.949 0.896 0.824 0.885 0.763
50 40.07 1.941 0.019 0.379 1.752 1.629 1.517 1.684 1.416
60 116.0 5.332 0.041 1.105 5.490 5.095 4.872 5.445 4.686
70 283.5 12.68 0.054 2.118 10.69 9.818 10.72 10.64 9.118
80 440.0 20.77 0.079 3.533 18.18 16.21 17.49 17.98 15.31
90 833.9 42.17 0.102 6.121 28.89 26.50 30.04 28.18 24.41
100 2478.4 82.54 0.129 11.23 59.78 57.50 66.97 60.57 53.24
HV throughout all the problem sets. Furthermore, the
HV value increases as the problem size gets larger.
On the other hand, the highest HV value of FPBH is
97.49% in the smallest problem class (n=5), and it de-
creases as the problem size increases. This suggests
that MOAP is not a suitable benchmark problem for
MOIP heuristics.
In the case of the MOKP, FPBH and the PR vari-
ants are competitive in terms of the number of so-
lutions for instances with fewer than or equal to 70
items (n 70). For n 80, FPBH generates more
solutions than PR variants. PR variants take less
computation time than FPBH in all instances. No-
tably, PRrand takes less than a quarter of the time
than FPBH does. Skipping SelectionRules and the
best move analysis, but instead relying on randomness
helps to reduce CPU time. Except for PI, embed-
ding the ImprovedND operation into the PR heuris-
tics does not bring any noticeable difference in the
number of solutions or run time. PI finds the most
solutions among all PR variants while its CPU time
increases fivefold. We also observe that it finds better
solutions while spending longer running time in Table
5. RD always produces the fewest solutions among
the heuristic methods. However, it takes consider-
ably less CPU time. For instance, it takes less than
1 second regardless of the problem size as it is seen
in Table 4. Although FPBH finds more solutions for
larger instances, every PR variant achieves a higher
HV than FPBH, except for the smallest problem class
with n=10. In particular, PI generates the highest
quality solutions throughout the experiments. RD also
outperforms FPBH for the instances with more than
30 items.
We observe that the RD and PR variants do not
show better performance than FPBH on the small-
est instances. The reason for this may be the struc-
ture of the test instances. When a problem size is
small, a small fractional value has a relatively big im-
pact on each objective. When a fractional value is
rounded down, the solution quality deteriorates com-
parably more on smaller instances. In addition, the
number of initially provided LB set solutions is lim-
ited in smaller instances. For these reasons, RD has a
A LP Relaxation based Matheuristic for Multi-objective Integer Programming
95
Table 5: Comparing HV(%) of algorithms on MOKP for
p=3, * indicates optimal values, best heuristic values are in
bold.
PR variants
n KS* FPBH RD PRrand PRsim PRdif PI PIsim PIdif
10 6.58 6.28 (95.4) 5.91 (89.8) 6.02 (91.5) 5.96 (90.6) 5.95 (90.4) 6.03 (91.6) 5.94 (90.3) 5.93 (90.1)
20 6.97 6.70 (96.1) 6.61 (94.8) 6.75 (96.8) 6.72 (96.4) 6.71 (96.3) 6.77 (97.1) 6.71 (96.3) 6.70 (96.1)
30 7.21 6.92 (96.0) 6.96 (96.5) 7.05 (97.8) 7.04 (97.6) 7.04 (97.6) 7.06 (97.9) 7.04 (97.6) 7.04 (97.6)
40 7.14 6.86 (96.1) 6.95 (97.3) 7.01 (98.2) 7.01 (98.2) 7.01 (98.2) 7.03 (98.5) 7.01 (98.2) 7.01 (98.2)
50 7.23 7.00 (96.8) 7.05 (97.5) 7.09 (98.1) 7.09 (98.1) 7.09 (98.1) 7.10 (98.2) 7.09 (98.1) 7.09 (98.1)
60 7.18 6.95 (96.8) 7.03 (97.9) 7.06 (98.3) 7.06 (98.3) 7.06 (98.3) 7.07 (98.5) 7.06 (98.3) 7.06 (98.3)
70 7.19 6.98 (97.1) 7.05 (98.1) 7.08 (98.5) 7.08 (98.5) 7.08 (98.5) 7.09 (98.6) 7.08 (98.5) 7.08 (98.5)
80 7.22 7.05 (97.6) 7.10 (98.3) 7.12 (98.6) 7.11 (98.5) 7.12 (98.6) 7.12 (98.6) 7.11 (98.5) 7.12 (98.6)
90 7.20 7.01 (97.4) 7.08 (98.3) 7.10 (98.6) 7.10 (98.6) 7.10 (98.6) 7.10 (98.7) 7.10 (98.6) 7.10 (98.6)
100 7.19 6.99 (97.2) 7.08 (98.5) 7.09 (98.6) 7.09 (98.6) 7.09 (98.6) 7.10 (98.7) 7.09 (98.6) 7.09 (98.6)
large HV gap. The quality of RD also influences that
of PR variants. If very small IR sets are provided, the
choice of the S
i
-S
g
pair is restricted. This causes a
limited number of new paths to be generated.
6 CONCLUSION
In this study, we propose a LP relaxation-based
matheuristic for three-objective binary integer pro-
gramming. The proposed algorithm relies on a high-
performing vector LP solver, Bensolve, which pro-
vides bound sets, and two simple heuristics, round-
ing down and path relinking (PR). In the compu-
tational study, we show that simple rounding down
can already find high-quality solutions in most in-
stances. After embedding PR with the first stage, the
proposed heuristic generates more solutions, which
show higher quality than that of FPBH in most prob-
lem classes. The number of solutions and CPU times
of the PR variants are similar to each other. Notably,
the PRrand algorithm takes less than a quarter of the
computation time FPBH does. The biggest advan-
tage of the proposed algorithm is that it can find high-
quality approximations fast, which also shows its ef-
fectiveness.
For future work, we plan to extend the algorithm
to deal with general MOIPs. Further, it could be tai-
lored to real-world applications such as supply chain
network design. As the size of real-world problems be
much larger, approaches to further reduce the com-
putation time will be investigated. In addition, the
proposed method can be embedded into an interactive
algorithm in order to facilitate decision making.
ACKNOWLEDGEMENTS
This research was funded in whole, or in part, by the
Austrian Science Fund (FWF) [P 31366-NBL]. For
the purpose of open access, the author has applied a
CC BY public copyright licence to any Author Ac-
cepted Manuscript version arising from this submis-
sion.
REFERENCES
Beausoleil, R. P., Baldoquin, G., and Montejo, R. A. (2008).
Multi-start and path relinking methods to deal with
multiobjective knapsack problems. Annals of Oper-
ations Research, 157(1):105–133.
Boland, N., Charkhgard, H., and Savelsbergh, M. (2017).
The quadrant shrinking method: A simple and effi-
cient algorithm for solving tri-objective integer pro-
grams. European Journal of Operational Research,
260(3):873–885.
Braekers, K., Hartl, R. F., Parragh, S. N., and Tricoire, F.
(2016). A bi-objective home care scheduling prob-
lem: Analyzing the trade-off between costs and client
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
96
inconvenience. European Journal of Operational Re-
search, 248(2):428–443.
Csirmaz, L. (2020). Inner approximation algorithm for
solving linear multiobjective optimization problems.
Optimization, pages 1–25.
de Souza Lima, F. M., Pereira, D. S., da Conceic¸
˜
ao, S. V.,
and de Camargo, R. S. (2017). A multi-objective ca-
pacitated rural school bus routing problem with het-
erogeneous fleet and mixed loads. 4OR, 15(4):359–
386.
Ehrgott, M. and Gandibleux, X. (2007). Bound sets
for biobjective combinatorial optimization problems.
Computers & Operations Research, 34(9):2674–
2694.
Fischetti, M., Glover, F., and Lodi, A. (2005). The feasi-
bility pump. Mathematical Programming, 104(1):91–
104.
Fonseca, C. M., Paquete, L., and L
´
opez-Ib
´
anez, M. (2006).
An improved dimension-sweep algorithm for the hy-
pervolume indicator. In 2006 IEEE international con-
ference on evolutionary computation, pages 1157–
1163. IEEE.
Glover, F. (1997). Tabu search and adaptive memory pro-
gramming—advances, applications and challenges. In
Interfaces in computer science and operations re-
search, pages 1–75. Springer.
Jaszkiewicz, J. (2005). Path relinking for multiple objec-
tive combinatorial optimization: Tsp case study. In
The 16th Mini-EURO Conference and 10th Meeting
of EWGT (Euro Working Group Transportation).
Kirlik, G. and Sayın, S. (2014). A new algorithm for gen-
erating all nondominated solutions of multiobjective
discrete optimization problems. European Journal of
Operational Research, 232(3):479–488.
Kolli, S. and Evans, G. W. (1999). A multiple objec-
tive integer programming approach for planning fran-
chise expansion. Computers & Industrial Engineer-
ing, 37(3):543–561.
Kovacs, A. A., Parragh, S. N., and Hartl, R. F. (2015). The
multi-objective generalized consistent vehicle routing
problem. European Journal of Operational Research,
247(2):441–458.
Leitner, M., Ljubi
´
c, I., Sinnl, M., and Werner, A. (2016). Ilp
heuristics and a new exact method for bi-objective 0/1
ilps: Application to fttx-network design. Computers
& Operations Research, 72:128–146.
L
¨
ohne, A. and Weißing, B. (2017). The vector linear
program solver bensolve–notes on theoretical back-
ground. European Journal of Operational Research,
260(3):807–813.
Mart
´
ı, R., Campos, V., Resende, M. G., and Duarte, A.
(2015). Multiobjective grasp with path relinking. Eu-
ropean Journal of Operational Research, 240(1):54–
71.
Mavrotas, G. and Florios, K. (2013). An improved version
of the augmented ε-constraint method (augmecon2)
for finding the exact pareto set in multi-objective in-
teger programming problems. Applied Mathematics
and Computation, 219(18):9652–9669.
¨
Ozpeynirci,
¨
O. and K
¨
oksalan, M. (2010). An exact
algorithm for finding extreme supported nondomi-
nated points of multiobjective mixed integer pro-
grams. Management Science, 56(12):2302–2315.
Pal, A. and Charkhgard, H. (2019a). A feasibility pump
and local search based heuristic for bi-objective pure
integer linear programming. INFORMS Journal on
Computing, 31(1):115–133.
Pal, A. and Charkhgard, H. (2019b). Fpbh: A feasibility
pump based heuristic for multi-objective mixed inte-
ger linear programming. Computers & Operations Re-
search, 112:104760.
Parragh, S. N., Doerner, K. F., Hartl, R. F., and Gandibleux,
X. (2009). A heuristic two-phase solution approach
for the multi-objective dial-a-ride problem. Networks:
An International Journal, 54(4):227–242.
Sawik, B., Faulin, J., Perez-Bernabeu, E., and Serrano-
Hernandez, A. (2016). Bi-objective optimization
models for green vrp approaches. ICIL 2016, page
247
`
A254.
Soylu, B. (2015). Heuristic approaches for biobjective
mixed 0–1 integer linear programming problems. Eu-
ropean Journal of Operational Research, 245(3):690–
703.
Zhang, W. and Reimann, M. (2014). A simple augmented
εconstraint method for multi-objective mathematical
integer programming problems. European Journal of
Operational Research, 234(1):15–24.
APPENDIX
Table 6 shows the comparison of three state-of-
the-art algorithms that can deal with multi-objective
linear programming. Bensolve (Ben) by L
¨
ohne
and Weißing (2017) and Inner solver (Inner) by
Csirmaz (2020) are implemented in C and pub-
licly available at http://www.bensolve.org/ and
https://github.com/lcsirmaz/inner, respectively. The
algorithm suggested by
¨
Ozpeynirci and K
¨
oksalan
(2010) (
¨
OK) is implemented in Julia. GLPK is used
as LP solver. The time limit of the experiment
is 3600 seconds. All the experiments are carried
out on a Quad-core X5570 Xeon CPUs @2.93GHz
with 48GB RAM. The figures are the average re-
sults of 10 test instances over 10 runs. As bench-
mark instances, we used the multi-objective assign-
ment problem (MOAP), the multi-objective knapsack
problem (MOKP), and multi-objective general inte-
ger linear programming problems (MOILP) which
are all generated by Kirlik and Sayın (2014) and
available at http://home.ku.edu.tr/ moolibrary/. Each
problem class is divided into subclasses. The sub-
classes are categorised by the number of items.
For the MOAP, it is categorised by 5/10/15/30/50,
whereas for the MOKP and MOILP, the subclasses
are 10/30/50/70/100. Each subclass has 10 instances;
A LP Relaxation based Matheuristic for Multi-objective Integer Programming
97
Table 6: Comparing algorithms on MOLP instances with p=3.
CPUtime(sec) #solved LP
Problem #item Ben Inner
¨
OK Ben Inner
¨
OK
MOAP
5 0.004 0.003 2.46 30.2 23.0 42.0
10 0.03 0.02 30.20 117.2 110.1 1916.8
15 0.10 0.07 844.04 237.6 230.5 12845.7
30 1.57 1.59 n/a. 1015 1008 n/a.
50 13.51 19.09 n/a. 2705 2698 n/a.
MOKP
10 0.003 0.002 4.27 46.0 39.0 300.7
30 0.02 0.01 320.98 223.1 216.0 7749.0
50 0.03 0.02 320.98 464.6 454.7 n/a.
70 0.07 0.07 n/a. 978.6 920.0 n/a.
100 0.16 0.14 n/a. 1962.0 1397.0 n/a.
MOILP
10 0.003 0.001 3.71
7
28.1 18.3 12.97
7
30 0.01 0.05 449.9
8
79.5 70.1 3481.8
8
50 0.02 0.01 1295.4
2
134.71 124.3 3764.5
2
70 0.04 0.02 32.1
2
178.2 167.0 857.0
2
100 0.07 0.03 3082.2
1
237.4 223.53 12835.0
1
thus, there are 50 instances per class in total. We re-
port the CPU time (sec) and the number of LPs. n/a.
indicates that the algorithm did not terminate within
the time limit.
number
indicates the number of in-
stances solved out of 10 instances.
Throughout the experiment, Ben and Inner are
highly competitive. In terms of the number of solved
LPs, the Inner solver comprehensively outperforms
the other two methods. Overall, the Inner solver
performs the best. However, it does not provide all
the bound set information we need for our heuristic.
Thus, we use Bensolve to obtain the initial bound sets.
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
98