A LP Relaxation based Matheuristic for Multi-objective Integer

Programming

Duleabom An

, Sophie N. Parragh

1 a

, Markus Sinnl

1 b

and Fabien Tricoire

2 c

Institute of Production and Logistics Management, Johannes Kepler University Linz,

Altenberger Straße 69, 4040 Linz, Austria

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 signiﬁcantly 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 conﬂicting 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

conﬂicting objectives. To ﬁnd the optimal location

for the new store, the model maximises the number

of customers while minimising conﬂict 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

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-

https://orcid.org/0000-0002-7428-9770

https://orcid.org/0000-0003-1439-8702

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 efﬁcient (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 deﬁned 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-

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

 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 ﬁrst 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 brieﬂy 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-

ﬁcient solutions in the same area often share common

features. In each phase, their heuristic obtains feasi-

ble solutions by ﬁxing 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 ﬁrst 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 uniﬁed view using minimisation objectives.

y = min{C(x) : Ax ≥ b,x ∈ {0,1}

(MOIP)

where x

, j = 1,2, . . .,n, is the vector of decision vari-

ables and X := {Ax ≥ b, x ∈ {0,1}

} is the feasible

set. C is a p × n objective function matrix where c

(k = 1,2,..., p), is the k

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

(x) ≤ c

′

) for all

k ∈ {1,..., p} and c

(x) < c

′

) 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 (efﬁcient)

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

(x) + w

(x) + · · ·+ w

(x)

∑

i=1

= 1, w

∈ (0,1),

where w

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

(i.e. optimising a convex combination of all objective

functions), it is called a supported efﬁcient solution.

Otherwise, it is a non-supported efﬁcient 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 efﬁcient 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

≤ 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

denotes whether task l is as-

signed to the agent r (x

=1) or not (x

=0). When a

task is allocated to an agent, the corresponding non-

negative costs c

,...,c

are incurred. The MOAP

can be stated as follows:

min

∑

r=1

∑

l=1

j = 1,..., p

(1)

s.t.

∑

l=1

= 1 r = 1, . . .,n

(2)

∑

r=1

= 1 l = 1,...,n

(3)

∈ {0,1} r,l = 1, ...,n.

(4)

The objective of the MOAP is to ﬁnd 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-

ﬁcient 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

, 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

denotes whether item r is selected

for the knapsack (x

=1) or not (x

=0). Each item r

also has proﬁts v

...v

. Here, W, v

, and w

are non-

negative integer values. The MOKP model is stated

as follows:

max

∑

r=1

j = 1,..., p

(5)

s.t.

∑

r=1

≤ W

(6)

∈ {0,1} r = 1,...,n

(7)

Equation (5) denotes the objective functions maximis-

ing the p total proﬁts 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 ﬁrst 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

and Weißing (2017) and Inner approximation by Csir-

maz (2020) have been investigated to ﬁnd LB sets

ﬁrst. 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 ﬁnds 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 ﬁltered 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

X: an empty list

3 IR ← RoundingDown(L);

4 i = 0;

5 for i< iteration limit do

6 candX ←PathRelinking(IR,S

,IGPair);

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

) and a guiding

solution (S

) 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 ﬁrst 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 speciﬁc, 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 ﬁrst

row (Neighbourhood) of Table 1. When the P matrix

(8) represents the proﬁts of the three objectives, the

sets of total proﬁts that correspond to the ﬁrst neigh-

bourhood are [7,6,8], [5,4,6], and [0,0,0]. Since the

dominance relation among the solution points is clear,

the ﬁrst 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

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 ﬁnds diverse solutions.

Depending on the speciﬁc problem, the following

components of PR can be designed differently in the

algorithm.

1) Building an initial solution set

2) Selecting S

and S

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

and S

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

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

-S

pairs during the PR process. Since the

performance of PR can vary depending on the choice

of the S

-S

pair, we investigate different ways, and

call SelectionRuleI. The ﬁrst approach is to select two

random solutions from IR. For the second and third

methods, S

is chosen at random in IR ﬁrst, then the

similarity between S

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

by ﬁnding the most simi-

lar solution to S

. The third approach ﬁnds the most

different solution to S

for S

. When S

and S

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

and S

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

and S

has already been chosen.

To decide the number of intermediate paths, we need

to ﬁnd the positions of the differing variable values in

and S

. They are stored in ∆items as indices (lines

6-7). Once ∆items is deﬁned, we generate the neigh-

bourhood by the following rules: For each index in

∆items, the variable values in S

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

(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

. 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

(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 ﬁnding 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

(line 14). If S

is feasible and not in-

cluded in IR, it is stored in the archive candX (line

16). The infeasible S

is not archived, though, it is still

used for the next iterations as it might help to ﬁnd ad-

ditional feasible solutions in yet unexplored parts of

the search region. IR is updated by adding S

(line

15). To prevent the case in which the same pair of ini-

tial and guiding solution is chosen, the current [S

]

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 ﬁgures out which path

shows the biggest improvement compared to the cur-

rent solution S

. Algorithm 3 shows a precise descrip-

tion of the operation.

Let ND be a set of non-dominated intermediate

points and ob jS

be the objective values of S

. 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

, 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

(line 13).

ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems

Algorithm 2: The framework of path relinking.

Input: IR, S

, S

, IGPair

1 IGPair: an archive of S

-S

pairs

2 candX ←

3 IGPair ←

4 Select S

from IR following

SelectionRuleI;

5 while S

̸= S

and [S

] /∈ 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

is chosen by SelectionRuleII

12 else

13 S

← randomly choose one solution

from the neighbourhood;

14 Feasibility check of S

;

15 if S

is feasible and S

/∈ IR then

16 candX ← S

;

17 IR ← IR ∪ S

;

18 IGPair ← [S

];

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 ﬁrst 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 ﬁnds the most similar solution

among the remaining solution set for the guiding

solution.

• PRdif: This version ﬁnds 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

, 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

[ 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

column

11 for i=1,...,|ND| do

12 ND degree[i] ← sum(rank table i

column)

13 S

← ND with the highest degree

Output: S

• 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 classiﬁed 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

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

be a set of k

objective values of the true

PF and y ∈ R

be an arbitrary non-dominated point

obtained from a heuristic algorithm. Then, the nor-

malised values of the obtained point are:

− min(Y

)

max(Y

) − min(Y

)

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 ﬁgures are average results

over 10 test instances. The ﬁgures 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 ﬁnds 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 (n≥15), 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 signiﬁ-

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

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 ﬁnds the most

solutions among all PR variants while its CPU time

increases ﬁvefold. We also observe that it ﬁnds 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 ﬁnds 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

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 inﬂuences that

of PR variants. If very small IR sets are provided, the

choice of the S

-S

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 ﬁnd high-quality solutions in most in-

stances. After embedding PR with the ﬁrst 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 ﬁnd 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 efﬁ-

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

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 ﬂeet 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.

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 ﬁnding 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 ﬁnding 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 ﬁgures 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

Table 6: Comparing algorithms on MOLP instances with p=3.

CPUtime(sec) #solved LP

Problem #item Ben Inner

OK Ben Inner

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