A LP relaxation based matheuristic for multi-objective integer programming
Duleabom An, Sophie N. Parragh, Markus Sinnl, Fabien Tricoire
aa r X i v : . [ m a t h . O C ] F e b A LP relaxation based matheuristic for multi-objective integerprogramming
Duleabom An , Sophie N. Parragh ,Markus Sinnl and Fabien Tricoire 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 { duleabom.an, sophie.parragh, markus.sinnl } @jku.at, [email protected] Keywords: Three-objective binary integer programming, matheuristics, path relinking, multi-objective knapsack problemAbstract: Motivated by the success of matheuristics in the single-objective domain, we propose a very simple linearprogramming-based matheuristic for three-objective binary integer programming. To tackle the problem, weobtain lower bound sets by means of the vector linear programming solver
Bensolve . Then, simple heuristicapproaches, such as rounding and path relinking, are applied to this lower bound set to obtain high-qualityapproximations of the optimal set of trade-off solutions. The proposed algorithm is compared to a recentlysuggested 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 benchmarkinstances for the three-objective knapsack problem.
Many real-world optimisation problems involvemultiple conflicting objectives, concerning, e.g.,costs, environmental impact or service level, andcan be formulated as integer linear programs.Kolli and Evans (1999), for example, deal with facil-ity location problems for a franchise company. Whena new franchise is launched, the franchisee and thefranchisor have conflicting objectives. To find the op-timal location for the new store, the model maximisesthe number of customers while minimising conflictamong existing franchises. Sawik et al. (2016) inves-tigate vehicle routing problems occurring in the logis-tics of the food industry. The introduced model min-imises both the total travel distance and CO emis-sions. A home care routing and scheduling problem isinvestigated by Braekers et al. (2016). Given a groupof nurses and a number of care tasks at the patients’home locations, the model assigns tasks to each nursewhile minimising operating cost and client inconve-nience with respect to the timing of the visit and thenurse performing the task. Kovacs et al. (2015) study a https://orcid.org/0000-0002-7428-9770 b https://orcid.org/0000-0003-1439-8702 c https://orcid.org/0000-0002-3700-5134 a multi-objective consistent vehicle routing problem.To provide consistent service in the routing indus-try, they maximise driver consistency and arrival timeconsistency 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 thesingle objective domain, we present a very sim-ple matheuristic for three-objective integer program-ming. To the best of our knowledge, only one othermatheuristic for three-objective integer program-ming has been developed so far (Pal and Charkhgard,2019b). The algorithm we propose relies on a lowerbound set which is obtained from the linear pro-gramming (LP) relaxation, using the vector linearprogramming 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 obtainhigh-quality approximations of the true
Pareto fron-tier ( PF ).The contribution of this paper is twofold: We show that, in the case of the three-objectiveassignment problem, since the extreme solutionswhich Bensolve produces are integer, it alreadyprovides high-quality approximations of the truePareto frontier, even without additional ingredi-ents.- We propose the first LP relaxation-basedmatheuristic algorithm for three-objective integerprogramming, 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 backgroundfor MOIP. The proposed algorithms are described inSection 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.
Over the past years a number of genericexact methods for solving MOIP have been pro-posed (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 inthe single objective domain, only comparably fewcontributions on matheuristc methods exist.A heuristic and a matheuristic approach for binaryinteger linear programming are proposed by Soylu(2015). The methods are variants of variable neigh-bourhood search and local branching. Both algo-rithms collect segments of the PF during the searchand combine them at the end. However, the pro-posed algorithms are not compared with a bench-mark method and struggle with generating a high-quality approximation of the true PF in quite a fewproblem classes. To deal with bi-objective binary IP,Leitner et al. (2016) suggest an exact method and amatheuristic framework. In each phase, their heuris-tic obtains feasible solutions by fixing a large num-ber of variables and reducing the associated feasi-ble region, respectively. Based on the feasibilitypump (FP) idea introduced by Fischetti et al. (2005),Pal and Charkhgard (2019a) suggested a FP basedheuristic for bi-objective IP and extended the methodto higher dimensions in (Pal and Charkhgard, 2019b).The authors suggest a two-stage heuristic algorithmfor MO mixed integer linear programming. The au-thors focus on finding diverse non-dominated points in the first stage by using the weighted sum and localsearch method. In the second stage, they rely on localbranching. To the best of our knowledge, this is theonly matheuristic method that deals with MOIP withmore than two objectives. Therefore, we use it as abenchmark and refer to it as FPBH. The problem we consider is a MOIP, with binaryinteger decision variables and three objectives. In thefollowing, we state the MOIP in its general form.We present a unified view using minimisation ob-jectives (any maximisation objective function can beconverted into a minimisation one by multiplying itby -1). y = min { C ( x ) : Ax ≥ b , x ∈ { , } n } , (MOIP)where x j , j = , , . . . , n , is the vector of decision vari-ables and X : = { Ax ≥ b , x ∈ { , } n } is the feasibleset. C is a p × n objective function matrix where c k ,( k = , , . . . , p ), is the k th row of C . In our case, p = y is a set of points in the objective space (cri-terion space), each of which corresponds to at leastone solution vector x ∈ X . A is an m × n constraintmatrix and b is the vector of right-hand-side valuesfor these constraints. In MO optimisation, the quality of a solution is de-termined by
Pareto dominance. Suppose there aretwo solutions x and x ′ of Problem (MOIP). Then, x dominates x ′ if and only if c k ( x ) ≤ c k ( x ′ ) for all k ∈ { , . . . , 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 ofall non-dominated points is called the Pareto front .If a solution x ∗ can be found by a convex combina-tion of all objective functions, it is called a supported efficient solution. Ohterwise, it is a non-supported ef-ficient solution . The notion of bound set was introduced byEhrgott and Gandibleux (2007). In the context ofroblem (MOIP), a lower bound (LB) set and an up-per bound (UB) set provide information about thevariable range that efficient solutions of the MOIP canattain. One common way to obtain a LB set for a min-imisation problem is to solve the LP relaxation of theoriginal problem. For constructing the LP relaxationof Problem (MOIP), the integrality conditions are re-moved, i.e. 0 ≤ x j ≤ j = , . . . , n and Bensolve canthen 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 thesesolutions as LB set.
The multi-objective assignment problem (MOAP)and the multi-objective knapsack problem (MOKP)are commonly used benchmark problems in MO. Weuse standard benchmark instances of these two prob-lems for our computational experiments.
The well-known assignment problem is a type oftransportation problem: a certain number of tasks l ∈ { , . . . , n } and agents r ∈ { , . . . , n } are given.The decision variable x rl denotes whether task l isassigned to the agent r ( x rl =1) or not ( x rl =0). Whena task is allocated to an agent, the correspondingnon-negative costs c rl , . . . , c prl are incurred. TheMOAP can be stated as follows:min n ∑ r = n ∑ l = c jrl x rl j = , . . . , p (1)s.t. n ∑ l = x rl = r = , . . . , n (2) n ∑ r = x rl = l = , . . . , n (3) x rl ∈ { , } r , l = , . . . , 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 eachagent is assigned to only one task. Equation (3) limitseach task to be assigned to one agent only.It is well-known that the constraint matrix of APis totally unimodular. According to the property, wecan conclude that every solution of an LP relaxationis naturally an integer vector. Although it is restrictedto supported efficient solutions as Bensolve is based on Benson’s outer approximation, it already providesa 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 suitablebenchmark for MOIP heuristics. In the multi-objective knapsack problem, a set ofitems is given, each with a certain weight w r , andwe must select a subset of these items such that thetotal weight does not exceed a given capacity W .The decision variable x r denotes whether item r isselected for the knapsack ( x r =1) or not ( x r =0). Eachitem r also has profits v r . . . v pr . Here, W , v r , and w r are non-negative integer values. The MOKP model isstated as follows:max n ∑ r = v jr x r j = , . . . , p (5)s.t. n ∑ r = w r x r ≤ W (6) x r ∈ { , } r = , . . . , n (7)Equation (5) denotes the objective functions maximis-ing the p total profits of selected items. Equation (6) isthe capacity constraint. The total weight of the itemsplaced in the knapsack cannot exceed the given ca-pacity.In this paper, we convert a maximisation objectivefunction into a minimisation one by multiplying itby − This section provides the overview of the pro-posed algorithm and its main ingredients.
The proposed matheuristic algorithm follows a two-stage approach. At the first stage, we obtain LBsets and adapt them to feasible integer roundedsets. A variant of the weighted sum methodby ¨Ozpeynirci and K¨oksalan (2010),
Bensolve byL¨ohne and Weißing (2017) and
Inner approximation by Csirmaz (2020) have been investigated to find LBsets first. Among the three methods,
Bensolve showscompetitive performance in the experiment and pro-duces all the required bound set information we need.he 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 variableto produce a feasible solution. These are referred toas integer rounded (IR) sets. At the second stage, PR is employed to improve the solution quality. Untilreaching 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 limitof iterations of PR to the number of solutions of theIR set times 50 . Since the dominance relation is notchecked in ˜ X during the search, dominated solutionsare filtered after the algorithm terminates, and the setof integer feasible solutions X is returned. Algorithm1 describes the general framework of the proposed al-gorithm. Algorithm 1:
The LP relaxation-basedmatheuristic framework
Input: L : points describing an LB set candX : an archive of newly found feasible IPsolutions by PR ˜ X : an empty list IR ← RoundingDown ( L ); i = for i < iteration limit do candX ← PathRelinking ( IR , S i , S g , IGPair ); Update candX i = i + ˜ X ← DominanceCheck ( candX ); Output: ˜ X 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). Themain idea of it is that there should be common char-acteristics among high-quality solutions. The methodproduces new solutions by exploring solution ”paths”between pairs of known solutions. To generate a newsolution, PR chooses two solutions from a set of ini-tial solutions; an initiating solution ( S i ) and a guidingsolution ( S g ) represent the starting and ending pointsof the path, respectively. Then it explores the trajecto-ries that link the pair of solutions in a neighbourhoodspace. In each step, a certain number of intermedi-ate paths, called a neighbourhood, are created. Forthe next step, the new initiating solution is selected in there. As the search proceeds, more attributes of theguiding 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 toMO problems such as the travelling salesman prob-lem (Jaszkiewicz, 2005), the dial-a-ride problem(Parragh et al., 2009), and a school bus routing prob-lem (de Souza Lima et al., 2017). It also produceshigh-quality solutions on large MOKP instanceswithin a reasonable timeframe (Beausoleil et al.,2008; Mart´ı et al., 2015). For these reasons, we com-bine PR with the first stage approach.For illustration purposes, we provide a short ex-ample describing the PR process on a MOKP withthree objectives. Suppose we have a four-item MOKPin which the initiating solution is (0 0 1 0) and theguiding 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 followingrules. To create one path, one variable value in theinitiating solution is switched. To be specific, if theitem is placed in the knapsack and its value is 1, itis taken out of the knapsack and the value changes to0. If the item is not chosen and the value equals 0,it is placed in the knapsack and the value changes to1. This procedure applies to all the different variablevalues. The outcome of the process is seen in the firstrow (Neighbourhood) of Table 1. When the P matrix(8) represents the profits of the three objectives, thesets of total profits that correspond to the first neigh-bourhood are [7,6,8], [5,4,6], and [0,0,0]. Since thedominance relation among the solution points is clear,the first path (1 0 1 0) is selected as the new initiatingsolution (marked with * in Table 1). This process re-peats until the initiating solution becomes identical tothe guiding solution. Table 1 illustrates the transform-ing process of PR in this example. P = (8)Initiating Guiding Neighbourhood0 0 1 Table 1:
Path relinking procedure
In our PR process, an infeasible intermediate so-lution is not stored in the solution archive, but stillcan be a new initiating solution. This is for letting thelgorithm explore a broader search region so that itpossibly finds diverse solutions.Depending on the specific problem, the followingcomponents of PR can be designed differently in thealgorithm.1) Building an initial solution set2) Selecting S i and S g
3) Generating intermediate paths (a neighbourhood)4) Choosing the next initiating solutionIn 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. Since the performanceof 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 tworandom solutions from IR . For the second andthird methods, S i is chosen at random in IR first,then the similarity between S i and all the other so-lutions in IR is calculated. The similarity is mea-sured by counting the number of variables whichhave equal values. Thus, higher values imply agreater similarity . The second method selects S g by finding the most similar solution to S i . Thethird approach finds the most different solution to S i for S g . When S i and S g are similar, fewer inter-mediate paths are likely to be created. However,when they are different, the PR process can takelonger as it explores relatively diverse intermedi-ate paths.3) A new neighbourhood is generated as many timesas 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 forthe proposed algorithm is given in Algorithm 2.The algorithm maintains two archives candX and IGPair . Each archive stores newly found solutionsand used S i - S g pairs during the PR process. Under theiteration limit, the whole process repeats. Firstly, S i and S g are chosen by SelectionRuleI (line 7). A PR iteration ends when one of the following two termina-tion conditions is met (line 8);- 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 needto find the positions of the differing variable values in S i and S g . They are stored in ∆ items as indices (lines9-10). 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 11). Let us suppose items 1 and 5 have a dif-ferent value, ∆ items = { , } . If item 1 is selected in S i (variable value=1), then it is taken out of it. If itis out of the knapsack (variable value=0), then it isadded 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 neighbouris not systematically selected. Rather, it is selectedwith 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 bestneighbour, we analyse the dominance relation amongintermediate solutions and select the best solution as S i (line 14). For this, SelectionRuleII is introduced. Ifone non-dominated solution exists in the neighbour-hood, it becomes the next initiating solution. If thereare mutually non-dominating solutions, we check the improved ratio of each solution point. The approachfor finding the most-improved point is explained indetail in Section 4.2.1. For the case where we do notwant to select the best neighbour, the analysis is notneeded and simply one random solution is selectedfrom the neighbourhood (lines 15-16). Before mov-ing on to the next iteration, the algorithm checks thefeasibility of S i (line 18). If S i is feasible and not in-cluded in IR , it is stored in the archive candX (line19). The infeasible S i is not archived, though, it is stillused for the next iterations as it might help to find ad-ditional feasible solutions in yet unexplored parts ofthe search region. A newly found feasible solution isadded to IR (line 20). To prevent the case in which thesame pair of initial and guiding solution is chosen, thecurrent [ S i , S g ] pair is archived after every PR iteration(line 21). After post-processing that filters dominatedsolutions (line 23), The algorithm returns the set ofinteger feasible solutions ˜ X . ImprovedND operation
The
ImprovedND operation figures out which pathshows 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 intermediatepoints and ob jS i be the objective values of S i . Torecord 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 (lines7-10). The rankings of each non-dominated point areadded up (lines 11-12) and stored in ND degree . Thenon-dominated path with the highest degree is set to S i (line 13). lgorithm 2: LP relaxation-basedmatheuristic
Input: IR candX : an archive of newly found feasible IPsolutions IGPair : an archive of S i - S g pairs ˜ X : an empty list candX ← /0 ; IGPair ← /0 ; i =0; for i < iteration limit do Select S i , S g from IR following SelectionRuleI ; while S i = S g and [S i , S g ] / ∈ IGPair do ∆ items ← index set of differentvariable values; n ← ∆ items ; Create n neighbourhood; The best move analysis: if rand() < then S i is chosen by SelectionRuleII else S i ← randomly choose onesolution from theneighbourhood; Feasibility check of S i ; if S i is feasible and S i / ∈ IR then candX ← S i ; Update IR by adding S i ; IGPair ← [ S i , S g ] ; i = i + ˜ X ← DominanceCheck ( candX ); Output: ˜ X We evaluate the performance of the differentversions of our heuristic method, focusing on thecomparison with FPBH (Pal and Charkhgard, 2019b).The variant with rounding down is named RD . Thefollowing PR variants are tested. The first three PR variants randomly choose the next initiating solutionduring the iterations.• PRrand : This version randomly selects both theinitiating and guiding solutions from the initial so-lution set.• PR sim: This version randomly chooses the initi-ating solution, then finds the most similar solutionamong the remaining solution set for the guidingsolution.• PR dif: This version finds the most different solu- Algorithm 3:
ImprovedND
Input: ob jS i , ND ratio table : | ND | × p matrix rank table : | ND | × p matrix ND degree : list of size | ND | for i=1,..., | ND | do for j=1,...,p do ratio table [ i ][ j ] = ND [ i ][ j ] ob jS i [ j ] for i=1,..., | ND | do for j=1,...,p do rank table [ i ][ j ] = rank of ND [ i ][ j ] in j th column for i=1,..., | ND | do ND degree [ i ] ← sum( rank table i th column) S i ← ND with the highest degree Output: S i tion to the already chosen initiating solution forthe guiding solution.The three other variants consider the improvementamong mutually non-dominating intermediate pathsto choose the next initiating solution.• PI : This version selects the most improved inter-mediate solution as the next initiating solution.• PI sim: This variant uses ImprovedND within PR sim.• PI dif: This variant uses ImprovedND in PR dif.The proposed algorithms use Bensolve byL¨ohne and Weißing (2017) to obtain the bound sets.The heuristic integration ( rounding down and PR ),is implemented in Julia. For the benchmark algo-rithm, we used the Julia implementation of FPBH(with the default setting) which is publicly availableat https://github.com/aritrasep/FPBH.jl . Allexperiments of matheuristics are carried out on In-tel® Core™ i5-8250U CPU running at 1.60GHz with16GB RAM. The exact MOIP solver proposed byKirlik and Sayın (2014) (KS) is also used in the ex-periment to obtain the true PF for comparison pur-pose. KS is run on Quad-core X5570 Xeon [email protected] with 48GB RAM. The KS results are forreference only (i.e. not for benchmarking).The test instances we use are the same ones onwhich FPBH is tested. The instances were generatedby Kirlik and Sayın (2014) and are publicly availableat http://home.ku.edu.tr/˜moolibrary/ . Eachproblem class has 100 instances divided into 10 sub-lasses, each of which contains 10 instances. MOAPinstances 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 itemswhich varies from 10 to 100 in increments of 10. One widely used indicator to measure the quality ofa 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 isthe “worst possible” point in the objective space. Inthis study, all the HV values are calculated with nor-malised objective values.Let Y kN be a set of k th objective values of the true PF and y ∈ R p be an arbitrary non-dominated pointobtained from a heuristic algorithm. Then, the nor-malised values of the obtained point are: y k − min ( Y kN ) max ( Y kN ) − min ( Y kN ) k = , . . . , p . As all non-dominated points are normalised, theirvalues exist in [0,1]. Therefore, the reference point is(1,1,1). Higher HV values indicate a better approx-imation. We used the publicly available HV com-puting program provided by Fonseca et al. (2006) at http://lopez-ibanez.eu/hypervolume to obtain HV values.
We report the following results of each algorithm:the number for solutions ( | Y | ), CPU time (sec), HVvalue, and HV as a percentage of the HV indicatorvalue for the exact non-dominated set as provided byKirlik and Sayın (2014). All the figures are averageresults over 10 test instances. The figures of PR vari-ants and FPBH are averaged over 10 runs for eachinstance because they have random components. Theresults of the experiments on MOAP and MOKP arereported in Tables 2-5. Bensolve already finds integer feasible solutionsfor all the MOAP instances. In addition, it showsbetter performance than FPBH in all the subclasses.Therefore, we do not apply our matheuristic to MOAPinstances.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 greateras the size of the instances grows. For example, forthe instances with more than 15 tasks (n ≥ Bensolve is more than doublethat of FPBH. Further, solutions are found in signifi-cantly less time. Hence, not only the quantity but alsothe quality of the solutions of
Bensolve are better thanFPBH. It reaches more than 99% of the maximumHV throughout all the problem sets. Furthermore, theHV value increases as the problem size gets larger.On the other hand, the highest HV value of FPBH is97.49% in the smallest problem class (n=5), and it de-creases as the problem size increases. This suggeststhat MOAP is not a suitable benchmark problem forMOIP heuristics.In the case of the KP, FPBH and the PR variantsare competitive in terms of the number of solutionsfor instances with fewer than or equal to 70 items (n ≤ n ≥
80, FPBH generates more solutions than PR variants. PR variants take less computation timethan FPBH in all instances. Notably, PRrand takesless than a quarter of the time than FPBH does. Skip-ping
SelectionRules and the best move analysis, butinstead relying on randomness helps to reduce CPUtime. Except for PI , embedding the ImprovedND operation into the PR heuristics does not bring anynoticeable difference in the number of solutions orrun time. PI finds the most solutions among all PR variants while its CPU time increases fivefold. Wealso observe that it finds better solutions while spend-ing longer running time in Table 5. RD always pro-duces the fewest solutions among the heuristic meth-ods. However, it takes considerably less CPU time.For instance, it takes less than 1 second regardless ofthe problem size as it is seen in Table 4. AlthoughFPBH finds more solutions for larger instances, every PR variant achieves a higher HV than FPBH, exceptfor the smallest problem class with n=10. In particu-lar, PI generates the highest quality solutions through-out the experiments. RD also outperforms FPBH forthe instances with more than 30 items.We observe that the RD and PR variants do notshow 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 issmall, a small fractional value has a relatively big im-pact on each objective. When a fractional value isrounded down, the solution quality deteriorates com-parably more on smaller instances. In addition, thenumber of initially provided LB set solutions is lim-ited in smaller instances. For these reasons, RD has alarge HV gap. The quality of RD also influences thatof PR variants. If very small IR sets are provided, the | Y | CPUtime(sec) HV HV(%)KS* FPBH
Bensolve
KS FPBH
Bensolve
KS FPBH
Bensolve
FPBH
Bensolve
10 176.8 21.0
15 674.9 40.4
20 1860.5 62.5
25 3567.8 90.6
30 6181.3 140.0
35 8972.3 163.2
40 14679.7 242.9
45 17702.2 238.6
50 24916.8 337.5
Table 2: Comparing algorithm performance on MOAP for p =3, * indicates optimal values, best heuristic values are in bold. PR variantsn KS* FPBH RD PRrand PR sim PR dif PI PI sim PI dif10 9.8 5.1 4.3 5.7 4.8 4.8 Table 3: Comparing | Y | of algorithms on MOKP for p =3, * indicates optimal values, best heuristic values are in bold. PR variantsn KS* FPBH RD PRrand PR sim PR dif PI PI sim PI dif10 0.140 0.023 Table 4: Comparing CPU time (sec) of algorithms on MOKP for p =3, * indicates optimal values, best heuristic values are inbold. choice of the S i - S g pair is restricted. This causes alimited number of new paths to be generated. In this study, we propose a LP relaxation-basedmatheuristic for three-objective binary IP. The pro- posed algorithm relies on a high-performing vec-tor LP solver,
Bensolve , which provides bound sets, rounding down and PR . In the computational study,we show that simple rounding down can already findhigh-quality solutions in most instances. After em-bedding PR with the first stage, the proposed heuris-tic generates more solutions, which show higher qual-ity than that of FPBH in most problem classes. The R variantsn KS* FPBH RD PRrand PR sim PR dif PI PI sim PI dif10 6.58
90 7.20 7.01 (97.4) 7.08 (98.3) )100 7.19 6.99 (97.2) 7.08 (98.5) 7.09 (98.6) 7.09 (98.6) 7.09 (98.6)
Table 5: Comparing HV(%) of algorithms on MOKP for p =3, * indicates optimal values, best heuristic values are in bold. number of solutions and CPU times of the PR vari-ants are similar to each other. Notably, the PRrand algorithm takes less than a quarter of the computa-tion time FPBH does. The biggest advantage of theproposed algorithm is that it can find high-quality ap-proximations fast, which also shows its effectiveness.For future work, we plan to extend the algorithm todeal with general MOIPs.
ACKNOWLEDGEMENTS
This research was funded in whole, or in part, bythe Austrian Science Fund (FWF) [P 31366-NBL].For the purpose of open access, the author has ap-plied a CC BY public copyright licence to any AuthorAccepted Manuscript version arising from this sub-mission.
REFERENCES
Beausoleil, R. P., Baldoquin, G., and Montejo, R. A. (2008).Multi-start and path relinking methods to deal withmultiobjective 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 clientinconvenience.
European Journal of Operational Re-search , 248(2):428–443.Csirmaz, L. (2020). Inner approximation algorithm forsolving 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. , 15(4):359–386.Ehrgott, M. and Gandibleux, X. (2007). Bound setsfor 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 , 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 Meetingof EWGT (Euro Working Group Transportation) .Kirlik, G. and Sayın, S. (2014). A new algorithm for gen-erating all nondominated solutions of multiobjectivediscrete optimization problems.
European Journal ofOperational 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). Themulti-objective generalized consistent vehicle routingproblem.
European Journal of Operational Research ,247(2):441–458.Leitner, M., Ljubi´c, I., Sinnl, M., and Werner, A. (2016). Ilpheuristics and a new exact method for bi-objective 0/1ilps: Application to fttx-network design.
Computers& Operations Research , 72:128–146.L¨ohne, A. and Weißing, B. (2017). The vector linearprogram solver bensolve–notes on theoretical back-round.
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 versionof the augmented ε -constraint method (augmecon2)for finding the exact pareto set in multi-objective in-teger programming problems. Applied Mathematicsand Computation , 219(18):9652–9669.¨Ozpeynirci, ¨O. and K¨oksalan, M. (2010). An exactalgorithm 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 pumpand local search based heuristic for bi-objective pureinteger linear programming.
INFORMS Journal onComputing , 31(1):115–133.Pal, A. and Charkhgard, H. (2019b). Fpbh: A feasibilitypump 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 approachfor 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 optimizationmodels for green vrp approaches.
ICIL 2016 , page247 `A254.Soylu, B. (2015). Heuristic approaches for biobjectivemixed 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 mathematicalinteger programming problems. European Journal ofOperational 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) byL¨ohne and Weißing (2017) and
Inner solver (In-ner) by Csirmaz (2020) are implemented in C andpublicly available at and https://github.com/lcsirmaz/inner ,respectively. The algorithm suggested by¨Ozpeynirci and K¨oksalan (2010) ( ¨OK) is imple-mented in Julia. GLPK is used as LP solver. Thetime limit of the experiment is 3600 seconds. Allthe experiments are carried out on a Quad-coreX5570 Xeon CPUs @2.93GHz with 48GB RAM. The figures are the average results of 10 test in-stances over 10 runs. As benchmark instances,we used the multi-objective assignment problem(MOAP), the multi-objective knapsack problem(MOKP), and multi-objective general integer linearprogramming problems (MOILP) which are allgenerated by Kirlik and Sayın (2014) and avail-able at http://home.ku.edu.tr/˜moolibrary/ .Each problem class is divided into subclasses. Thesubclasses 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 subclassesare 10/30/50/70/100. Each subclass has 10 instances;thus, there are 50 instances per class in total. Wereport the CPU time (sec) and the number of LPs.n/a. indicates that the algorithm did not terminatewithin the time limit. ∗ number indicates the number ofinstances solved out of 10 instances.Throughout the experiment, Ben and Inner arehighly competitive. In terms of the number of solvedLPs, the Inner solver comprehensively outperformsthe other two methods. Overall, the
Inner solver performs the best. However, it does not provide allthe bound set information we need for our heuristic.Thus, we use
Bensolve to obtain the initial bound sets.PUtime(sec) n/a.50 n/a.MOKP 10 0.003 n/a.70 n/a. 978.6 n/a.100 0.16 n/a. 1962.0 n/a.MOILP 10 0.003 . ∗ . ∗ . ∗ . ∗
50 0.02 . ∗ . ∗
70 0.04 . ∗ . ∗
100 0.07 . ∗ . ∗ Table 6: Comparing algorithms on MOLP instances, p =3=3