Short rational generating functions for multiobjective linear integer programming
Abstract.
This paper presents algorithms for solving multiobjective integer programming problems. The algorithm uses Barvinok’s rational functions of the polytope that defines the feasible region and provides as output the entire set of nondominated solutions for the problem. Theoretical complexity results on the algorithm are provided in the paper. Specifically, we prove that encoding the entire set of nondominated solutions of the problem is polynomially doable, when the dimension of the decision space is fixed. In addition, we provide polynomial delay algorithms for enumerating this set. An implementation of the algorithm shows that it is useful for solving multiobjective integer linear programs.
Key words and phrases:
Multiple objective optimization, integer programming, generating functions2000 Mathematics Subject Classification:
90C92, 90C10, 05A152.cm2.cm2.cm2.cm
Introduction
Short rational functions were used by Barvinok [2] as a tool to develop an algorithm for counting the number of integer points inside convex polytopes, based in the previous geometrical papers by Brion [6],Khovanskii and Puhlikov [22], and Lawrence [24]. The main idea is encoding those integral points in a rational function in as many variables as the dimension of the space where the body lives. Let be a given convex polyhedron, the integral points may be expressed in a formal sum with , where . Barvinok’s aimed objective was representing that formal sum of monomials in the multivariate polynomial ring , as a “short” sum of rational functions in the same variables. Actually, Barvinok presented a polynomialtime algorithm when the dimension, , is fixed, to compute those functions. A clear example is the polytope : the long expression of the generating function is , and it is easy to see that its representation as sum of rational functions is the well known formula .
Brion proved in 1988 [6], that for computing the short generating function of the formal sum associated to a polyhedron, it is enough to do it for tangent cones at each vertex of . Barvinok applied this function to count the number of integral points inside a polyhedron , that is, , that is not possible to compute using the original expression, but it may be obtained using tools from complex analysis over the rational function .
The above approach, apart from counting lattice points, has been used to develop some algorithms to solve, exactly, integer programming. Actually, De Loera et al [9] and Woods and Yoshida [35] presented different methods to solve this family of problems using Barvinok’s rational function of the polytope defined by the constraints of the given problem.
The goal of this paper is to present new methods for solving multiobjective integer programming problems. In contrast to usual integer programming problems, in multiobjective problems there are at least two (and in this case the problem is called biobjective) or more objective functions to be optimized.
The importance of multiobjective optimization is not only due to its theoretical implications but also to its many applications. Witnesses of that are the large number of realworld decision problems that appear in the literature formulated as multiobjective programs. Examples of them are flowshop scheduling (see [18]), analysis in Finance (see [13], Chapter 20), railway network infrastructure capacity (see [11]), vehicle routing problems (see [20, 29]) or trip organization (see [31]).
Multiobjective programs are formulated as optimization (we restrict ourselves without lost of generality to the maximization case) problems over feasible regions with at least two objective functions. Usually, it is not possible to maximize all the objective functions simultaneously since objective functions induce a partial order over the vectors in the feasible region, so a different notion of solution is needed. A feasible vector is said to be a nondominated (or Pareto optimal) solution if no other feasible vector has componentwise larger objective values. The evaluation through the objectives of a nondominated solution is called efficient solution.
This paper studies multiobjective integer linear programs (MOILP). Thus, we assume that there are at least two objective functions involved, the constraints that define the feasible region are linear, and the feasible vectors are integers.
Even if we assume that the objective functions are also linear, there are nowadays relatively few exact methods to solve general multiobjective integer and linear problems (see [13]). Some of them, as branch and bound with bound sets, which belong to the class of implicit enumeration methods, combine optimality of the returned solutions with adaptability to a wide range of problems (see for example [36, 37, 26, 27] for details). Other methods, as Dynamic Programming, are general methods for solving, not very efficiently, general families of optimization problems (see [21, 7]). A different approach, as the TwoPhase method (see [33]), looks for supported solutions (those that can be found as solutions of a singleobjective problem over the same feasible region but with objective function a linear combination of the original objectives) in a first stage and nonsupported solutions are found in a second phase using the supported ones. The Twophase method combines usual singlecriteria methods with specific multiobjective techniques.
Apart from those generic methods, there are specific algorithms for solving some combinatorial biobjective problems: biobjective knapsacks ([34]), biobjective minimum spanning tree problems ([30]) or biobjective assignment problems ([28]), as well as heuristics and metaheuristics algorithms that decrease the CPU time for computing the nondominated solutions for specific biobjective problems.
Nowadays, new approaches for solving multiobjective problems, using tools from Algebraic Geometry and Computational Algebra, have been proposed in the literature aiming to provide new insights into the combinatorial structure of the problems. This new research line seems to be prolific in a near future. An example of that is presented in [5] where a notion of partial Gröbner basis is given that allows to build a test family (analogous to the test set concept but for solving multiobjective problems) to solve general multiobjective linear integer programming problems.
Another witness of this trend is the recent work by Deloera et al. [10]. In this paper, the authors present several algorithms for multiobjective integer linear programs using generating functions. Nevertheless, their approach differs from ours in that their requires, in addition, to fix the dimension of the objective space to prove polynomiality of their algorithms, and their proofs are totally different. Moreover, no actual implementation of the algorithms is shown in that paper although it addresses an interesting shortest distance problems respect to a prespecified Pareto point.
In this paper, we also use rational generating function of polytopes for solving multiobjective integer linear programs.
In Section 1, the main results on Barvinok’s rational functions, which we use in our approach, are presented. Section 2 presents the multiobjective integer problem and the notion of dominance in order to clarify which kind of solutions we are looking for. The two following sections analyze different algorithms for solving general multiobjective problems. In Section 3, fixing the dimension of the decision space, a polynomial time algorithm that encodes the set of nondominated solutions of the problem as a short sum of rational functions is detailed. Next, a digging algorithm that computes the entire set of nondominated solutions using the multivariate Laurent expansion for the Barvinok’s function of the polytope defined by the constraints of the problem is given in Section 4. In that section, a polynomial delay algorithm for solving multiobjective problems is also presented (fixing only the dimension of the decision space).
Section 5 shows the results of a computational experiment and its analysis. Here, we solve biobjective knapsack problems, report on the performance of the algorithms and draw some conclusions on their results and their implications.
1. Barvinok’s rational functions
In this section, we recall some results on short rational functions for polytopes, that we use in our development. For details the interested reader is referred to [2, 3, 4].
Let be a rational polytope in . The main idea of Barvinok’s Theory was encoding the integer points inside a rational polytope in a “long” sum of monomials:
where .
The following results, due to Barvinok, allow us to reencode, in polynomialtime for fixed dimension, these integer points in a “short” sum of rational functions.
Theorem 1.1 (Theorem 5.4 in [2]).
Assume , the dimension, is fixed. Given a rational polyhedron , the generating function can be computed in polynomial time in the form
where is a polynomialsize indexing set, and where and for all and .
As a corollary of this result, Barvinok gave an algorithm for counting the number of integer points in . It is clear from the original expression of that this number is , but is a pole for the rational function, so, the number of integer points in the polyhedron is . This limit can be computed using residue calculation tools from elementary complex analysis.
Another useful result due to Barvinok and Wood [4], states that computing the short rational function of the intersection of two polytopes, given the respective short rational function for each polytope, is doable in polynomial time.
Theorem 1.2 (Theorem 3.6 in [4]).
Let , be polytopes in and . Let and be their short rational functions with at most binomials in each denominator. Then there exists a polynomial time algorithm that computes
with , where the are rational numbers and , are nonzero integral vectors for and .
In the proof of the above theorem, the Hadamard product of a pair of power series is used. Given and , the Hadamard product is the power series
The following Lemma is instrumental to prove Theorem 1.2.
Lemma 1.1 (Lemma 3.4 in [4]).
Let us fix . Then there exists a polynomial time algorithm, which, given functions and such that
(1) 
where and such that there exists with for all , computes a function in the form
with , and such that possesses the Laurent expansion in a neighborhood of and .
For proving Theorem 1.2, it is enough to assure that for given polytopes , their rational functions satisfy conditions (1). It is not difficult to ensure that the conditions are verified after some changes are done in the expressions for the short rational functions (for further details, the interested reader is referred to [4]).
Actually, with this result a general theorem can be proved ensuring that for a pair of polytopes, , there exists a polynomial time algorithm to compute, given the rational functions for and , the short rational function of any boolean combination of and .
Finally, we recall that one can find, in polynomial time, rational functions for polytopes that are images of polytopes with known rational function.
Lemma 1.2 (Theorem 1.7 in [2]).
Let us fix . There exists a number and a polynomial time algorithm, which, given a rational polytope and a linear transformation such that , computes the function for , in the form
where , and for all .
2. Multiobjective combinatorial optimization problems
In this section we present the problem to be solved as well as the new concept of solutions motivated by the nature of the problem.
A multiobjective integer linear program (MOILP) can be formulated as:
(2)  
with integers and non negative. Without loss of generality, we will consider the above problem in its standard form, i.e., the coefficient of the objective functions are nonnegative and the constraints are in equation form. In addition, we will assume that the constraints define a polytope (bounded) in . Therefore, from now on we deal with .
It is clear that Problem (2) is not a standard optimization problem since the objective function is a coordinate vector, thus inducing a partial order among its feasible solutions. Hence, solving the above problem requires an alternative concept of solution, namely the set of nondominated (or Paretooptimal) points.
A vector is said to be a nondominated (or Pareto optimal) solution of if there is no other feasible vector such that
with at least one strict inequality for some . If is a nondominated solution, the vector is called efficient. Note that is a subset of (decision space) and is a subset of (objectives space).
We will say that a dominated point, , is dominated by if for all .^{1}^{1}1We are denoting by the binary relation ”more than or equal to” and where it is assumed that at least one of the inequalities in the list is strict. We denote by the set of all nondominated solutions for (2) and by the image under the objective functions of , that is, .
From the objective function , we obtain a linear partial order on as follows:
Notice that since , the above relation is not complete. Hence, there may exist noncomparable vectors. We will use this partial order, induced by the objective functions of problem (2) as the input for the multiobjective integer programming algorithm developed in this paper.
Sometimes, the same efficient value is the image of several nondominated solutions. At this point, different problems can be tackled. We say that two nondominated solutions, and are equivalent if . Then, the solutions for are one of the following:

Complete set: A subset such that for all there is with .

Minimal complete set: A complete set with no equivalent solutions.

Maximal complete set: All equivalent solutions.
Through this paper, we are looking for the entire set of nondominated solutions, equivalently the maximal complete set for .
3. A short rational function expression of the entire set of nondominated solutions
We present in this section an algorithm for solving using Barvinok’s rational functions technique.
Theorem 3.1.
Let , , , , and assume that the number of variables is fixed. Suppose is a rational convex polytope in . Then, we can encode, in polynomial time, the entire set of nondominated solutions for in a short sum of rational functions.
Proof.
Using Barvinok’s algorithm, compute the following generating function in variables:
(3) 
where . is clearly a rational polytope. For fixed , the degrees, , in the monomial of represent the solutions dominated by .
Now, for any function , let be the projections of onto the  and variables, respectively. Thus encodes all dominated feasible integral vectors (because the degree vectors of the variables dominate them, by construction), and it can be computed from in polynomial time by Lemma 1.2.
Let be the set of extreme points of the polytope and choose an integer (we can find such an integer via linear programming). For this positive integer, , let be the rational function for the polytope , its expression is:
Define as above, the projection of onto the second set of variables as a function of the variables and the short generating function of . They are computed in polynomial time by Lemma 1.2 and Theorem 1.1 respectively. Compute the following difference:
This is the sum over all monomials where is a nondominated solution, since we are deleting, from the total sum of feasible solutions, the set of dominated ones.
This construction gives us a short rational function associated with the sum over all monomials with degrees being the nondominated solutions for . As a consequence, we can compute the number of nondominated solutions for the problem. The complexity of the entire construction being polynomial since we only use polynomial time operations among four short rational functions of polytopes (these operations are the computation of the short rational expressions for , and ). ∎
Remark 3.1.
To prove the above result one can use a different approach to compute the nondominated solutions assuming that there exists a polynomially bounded (for fixed dimension) feasible lower bound set, , for , i.e., a set of feasible solutions such that every nondominated solution is either one element in , or it dominates at least one the elements in .
First, compute the following operations with generating functions:
This is the sum over all monomials where , is a nondominated solution and is dominated by . In , each nondominated solution, , appears as many times as the number of feasible solutions that it dominates.
Next, compute a feasible lower bound set (see [16, 15]), . This way the set of nondominated solutions is encoded using the following construction:
Let be the following short sum of rational functions
Taking into account that for each , the element is common factor for and it is the unique factor where the variables appear, we can define , to be the sum of rational functions that encodes the nondominated solutions that dominate , . Therefore, the entire set of nondominated solutions for is encoded in the short sum of rational functions .
4. Digging algorithm for the set of nondominated solutions of MOILP
Section 3 proves that encoding the entire set of efficient solutions of MOILP can be done in polynomial time for fixed dimension. This is a compact representation of the solution concept. Nevertheless, one may be interested in an explicit description of this list of points. This task could be performed, by expanding the short rational expression which is ensured by Theorem 3.1, but it would require the implementation of all operations used in the proof. As far we know, they have never been efficiently implemented.
An alternative algorithm for enumerating the nondominated solutions of a multiobjective integer programming problem, which uses rational generating functions, is the digging algorithm. This algorithm is an extension of a heuristic proposed by Lasserre [23] for the singleobjective case.
Let and be as in Problem (2), and assume that is a polytope. Then, by Theorem 1.1, we can compute a rational expression for in the form
in polynomial time for fixed dimension, . Each addend in the above sum will be referred to as , .
If we make the substitution , in the monomial description we have , where are the rows in . It is clear that for enumerating the entire set of nondominated solutions, it would suffice to look for the set of leader terms, in the variables, in the partial order induced by , , of the multipolynomial . After the above changes we have:
(4) 
where Now, we can assume, wlog, that is negative or zero. If it were zero, then we could assume that is negative. Otherwise, we would repeat the argument until the first non zero element is found (it is assured that this element exists, otherwise the factor would not appear in the expression of the short rational function). Indeed, if the first non zero element were positive, we would make the change:
and the sign of the degree would be negative.
With these assumptions, the multivariate Laurent series expansion for each rational function, , in is
The following result allows us to develop a finite algorithm for solving using Barvinok’s rational generating functions.
Lemma 4.1.
Obtaining the entire set of nondominated solutions for a MOILP requires only an explicit finite, polynomially bounded (in fixed dimension) number of terms of the long sum in the Laurent expansion of .
Proof.
Let , and define , and . and are welldefined because , defined above, is non empty and bounded since, by construction, for each there exists such that .
Then, it is enough to search for the nondominated solutions in the finite sum
Let (resp. ) be the greatest (resp. smallest) value that appears in the nonzero absolute values of the entries in , , . Set . First, . Then, by applying Cramer’s rule one can see that is bounded above by . Thus, the explicit number of terms in the expansion of , namely , is polynomial, when the dimension, is fixed. ∎
The digging algorithm looks for the leader terms in the variables, with respect to the partial order induced by . At each rational function (addends in the above sum (4)) multiplications are done in lexicographical order in their respective bounded hypercubes. If the degree of a specific multiplication is not dominated by one of the previous factors, it is kept in a list; otherwise the algorithm continues augmenting lexicographically the lambdas. To simplify the search at each addend, the following consideration can be taken into account: if is dominated, then any term of the form , being componentwise larger than , is dominated as well.
The above process is done on each rational function that appears in the representation of . As an output we get a set of leader terms (for each rational function), that are the candidates to be nondominated solutions. Terms that appear with opposite signs will be cancelled. Removing terms in the list of candidates (to be nondominated solutions) implies consideration of those terms that were dominated by the cancelled ones. These terms are included in the current list of candidates and the process continues until no more terms are added.
At the end, some dominated elements may appear in the union of the final list. Deleting them in a simple cleaning process gives the list that contains only the entire set of nondominated solutions for the multiobjective problem.
Algorithm 1 details the pseudocode of the digging algorithm.
Recall that , where is the greatest value that appears in the nonzero absolute values of the entries in , , and is the least value among these values.
Taking into account Lemma 4.1 and the fact that Algorithm 1 never cycles, we have the following statement.
Theorem 4.1.
It is well known that enumerating the nondominated solutions of MOILP is NPhard and Phard ([12, 14]). Thus, one cannot expect to have very efficient algorithms for solving the general problem (when the dimension is part of the input).
In the following, we concentrate on a different concept of complexity that has been already used in the literature for slightly different problems. Computing maximal independent sets on graphs is known to be Phard ([17]), nevertheless there exist algorithms for obtaining these sets which ensure that the number of operations necessary to obtain two consecutive solutions of the problem is bounded by a polynomial in the problem input size (see e.g. [32]). These algorithms are called polynomial delay. Formally, an algorithm is said polynomial delay if the delay, which is the maximum computation time between two consecutive outputs, is bounded by a polynomial in the input size ([1, 19]).
In our case, a polynomial delay algorithm, in fixed dimension, for solving a multiobjective linear integer program means that once the first nondominated solution is computed, either in polynomial time a next nondominated solution is found or the termination of the algorithm is given as an output.
Next, we present a polynomial delay algorithm, in fixed dimension,
for solving multiobjective integer linear programming problems.
This algorithm combines the theoretical construction of Theorem
3.1 and a digging process in the Laurent
expansion of the short rational functions of the polytope
associated with the constraints of the problem.
The algorithm proceeds as follows.
Let be the short rational function that encodes the nondominated solutions (by Theorem 3.1, the complexity of computing is polynomial in fixed dimension). Make the changes , for , in . Denote by each of the rational functions of after the above changes. Next, the Laurent expansion over each rational function, , is done in the following way: (1) Check if contains nondominated solutions computing the Hadamard product of with . If does not contain nondominated solutions, discard it and set (termination); (2) if encodes nondominated solutions, look for an arbitrary nondominated solution (expanding ); (3) once the first nondominated solution, , is found, check if there exist more nondominated solutions encoded in the same rational function computing . If there are more solutions encoded in , look for them in . Repeat this process until no new nondominated solutions can be found in .
The process above describes the pseudocode written in Algorithm 2.
Theorem 4.2.
Assume is a constant. Algorithm 2 provides a polynomial delay procedure to obtain the entire set of nondominated solutions of .
Proof.
Let be the rational function that encodes the nondominated solutions of . Theorem 3.1 ensured that is a sum of short rational functions that can be computed in polynomial time.
Algorithm 2 digs separately on each one of the rational functions , ,that define . (Recall that ).
Fix . First, the algorithm checks whether encodes some nondominated solutions. This test is doable in polynomial time by Theorem 1.2. If the answer is positive, an arbitrary nondominated solution is found among those encoded in . This is done using digging and the Intersection Lemma. Specifically, the algorithm expands on the hyperbox and checks whether each term is nondominated. The expansion is polynomial, for fixed , since the number of terms is polynomially bounded by Lemma 4.1. The test is performed using the Hadamard product of each term with .
The process is clearly a polynomial delay algorithm. We use digging separately on each rational function that encodes nondominated points. Thus, the time necessary to find a new nondominated solution from the last one is bounded by the application of digging on a particular which, as argued above, is polynomially bounded. ∎
Instead of the above algorithm one can use a binary search procedure to solve multiobjective problem using short generating functions. In the worst case, digging algorithm may need to expand every nonnegative term to obtain the set of nondominated solutions. Therefore, as it is stated in Theorem 4.1, the number of steps to solve the problem can be polynomially bounded on . With a binary search approach, the number of steps to obtain consecutive solutions of our problem decreases to a number polynomially bounded on . A binary search approach was already used in [10]. Here, the novelty is that our analysis does not require to fix the dimension of the objective space whereas in [10] it was required.
The process is as follows. Let be defined as above. By construction . We proceed by dividing the hypercube into hypercubes of smaller dimensions, and recursively repeating the division process over those hypercubes containing at least one nondominated solution (until only one solution is included in each element of the partition), whereas those hypercubes that at a given stage of the process do not contain nondominated solutions are discarded for any further consideration.
The division process is done by bisecting each dimension. Testing for nondominated solutions on a given hypercube (at any stage of the process) is always done using the same tool based on Theorem 3.1. That result allows us to construct, in polynomial time in fixed dimension, the function that encodes all nondominated solutions. Moreover, it is easy to see that the short rational function that encodes the integer points in the hypercube , with , , is:
Thus, the Hadamard product, encodes the subset of nondominated solutions that lie in ; and hence by Barvinok’s theory we can also count, in polynomial time, the number of integer points encoded by (Lemma 3.4 in [4]).
The elements in our search space (hypercubes) are organized on a search tree and we use a depth first search strategy. Each node is a hypercube containing nondominated solutions. Descendants of a given node are hypercubes obtained bisecting the edges on the previous one (parent). It is clear that the maximum depth of the tree is . The above construction ensures that, provided that the set of nondominated solutions is nonempty, finding a first nondominated solution can be done testing at most nodes in the search tree. Since testing a node is polynomial, in fixed dimension, this operation is polynomial. Moreover, finding a new nondominated solution from a given one is also polynomial. Indeed, it consists of backtracking at most nodes until we find a branch containing nondominated points and then we have to explore, at most, nodes; or detecting that none of the branches contain solutions.
An illustrative example of this procedure is shown in Figure 1 where can be seen how the initial hypercube, , is divided successively in subhypercubes, until an isolated nondominated solution is located in one of them.
The finiteness of this procedure is assured since the number of times that the hypercube can be divided in subhypercubes is bounded by .
The pseudocode for this procedure is shown in Algorithm 3.
Theorem 4.3.
Assume is a constant. Algorithm 3 provides a polynomial delay (polynomially bounded on ) procedure to obtain the entire set of nondominated solutions of .
Remark 4.1.
The application of the above algorithm to the single criterion case provides an alternative proof of polynomiality for the problem of finding an optimal solution of integer linear problems, in fixed dimension.
Assume that the number of objectives, , is , and that there exists a unique optimal value for the problem. Applying Theorem 3.1 ensures that the optimal solution of the problem is found in polynomial time, if the dimension is fixed.
Remark 4.2 (Optimization over the set of nondominated solutions).
In practice, a decision maker expects to be helped by the solutions of the multiobjective problem. In many cases, the set of nondominated solutions is too large to make easily the decision, so it is necessary to optimize (using a new criterion) over the set of nondominated solutions.
With our approach, we are able to compute, in polynomial time for fixed dimension, a “short sum of rational functions”representation, , of the set of nondominated solutions of . This representation allows us to reoptimize with a linear objective, , based in the algorithms for solving singleobjective integer programming problems using Barvinok’s functions (see e.g. [35]) or the algorithm proposed in Remark 4.1. The above discussion proves that solving the problem of optimizing a linear function over the efficient region of a multiobjective problem is doable in polynomial time, for fixed dimension.
5. Computational Experiments
For illustrative propposes, a series of computational experiments have been performed in order to evaluate the behavior of a simple implementation of the digging algorithm (Algorithm 1). Computations of short rational functions have been done with Latte v1.2 [8] and Algorithm 1 has been coded in MAPLE 10 and executed in a PC with an Intel Pentium 4 processor at 2.66Gz and 1 GB of RAM. The implementation has been done in a symbolic programming language, available upon request, in order to make the access easy for the interested readers.
The performance of the algorithm was tested on randomly generated instances for biobjective (two objectives) knapsack problems. Problems from 4 to 8 variables were considered, and for each group, the coefficients of the constraint were randomly generated in . The coefficients of the two objective matrices range in and the coefficients of the right hand side were randomized in . Thus, the problems solved are in the form:
(5) 
The computational tests have been done on this way for each number of variables: (1) Generate 5 constraint vectors and right hand sides and compute the shorts rational functions for each of them; (2) Generate a random biobjective matrix and run digging algorithm for them to obtain the set of nondominated solutions.
Table 1 contains a summary of the average
results obtained for the considered knapsack multiobjective
problems. The second and third columns show the average CPU times
for each stage in the Algorithm: srf
is the CPU time for
computing the short rational function expression for the polytope
with LattE and modigging
the CPU time for running the
multiobjective digging algorithm for the problem. The total
average CPU times are summarized in the total column.
Columns latpoints and nosrf represent the number
of lattice points in the polytope and the number of short rational
functions, respectively. The average number of efficient solutions
that appear for the problem is presented under effic. The
problems have been named as knapN
where N
is the
number of variables of the biobjective knapsack problem.
problem  srf  latpoints  nosrf  modigging  effic  total 

knap4 
0.018  12.25  25.75  4.863  4.5  4.881 
knap5 
0.038  31  62.5  487.640  9.25  487.678 
knap6 
0.098  217.666  124.25  2364.391  7.666  2364.489 
knap7 
0.216  325  203  2869.268  20  2869.484 
knap8 
0.412  3478  342  10245.533  46  10245.933 
As can be seen in Table 1, the computation times are clearly divided into two steps (srf and modigging), being the most expensive the application of the digging algorithm (Algorithm 1). In all cases more than 99% of the total time is spent expanding the short rational function using “digging algorithm”.
The CPU times and sizes in the two steps are highly sensitive to the number of variables. It is clear that one cannot expect fast algorithm for solving MOILP, since all these problems are NPhard and #Phard. Nevertheless, this approach gives exact tools for solving any MOILP problem, independently of the combinatorial nature of the problem.
Finally, from our computational experiments, we have detected that an easy, promising heuristic algorithm could be obtained truncating the expansion at each rational function. That algorithm would accelerate the computational times at the price of obtaining only heuristics nondominated points.
Acknowledgements
This research has been partially supported by Junta de Andalucía grant number P06–FQM–01366 and by the Spanish Ministry of Science and Education grant number MTM2007–67433–C02–01. The authors acknowledge the useful comments received from J. De Loera and M. Köppe on an earlier version of this paper.
References
 [1] Arimura, H., Uno, T. (2005). A Polynomial Space and Polynomial Delay Algorithm for Enumeration of Maximal Motifs in a Sequence. ISAAC 2005: 724–737
 [2] Barvinok, A. A polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed, Mathematics of Operations Research , 19 (1994), 769–779.
 [3] Barvinok, A. and Pommersheim, J.E. An algorithmic theory of lattice points in polyhedra, in: New Perspectives in Algebraic Combinatorics (Berkeley, CA, 19961997), 91147, Math. Sci. Res. Inst. Publ. 38, Cambridge Univ. Press, Cambridge, 1999.
 [4] Barvinok, A. and Woods, K. Short rational generating functions for lattice point problems, Journal of the American Mathematical Society, 16 (2003), 957–979.
 [5] Blanco, V. and Puerto, J. (2007). Partial Gröbner bases for multiobjective combinatorial optimization. Submitted. arXiv:0709.1660
 [6] Brion, M. Points entiers dans les polyèdres convexes. Annales scientifiques de l’Ècole Normale Supèrieure Sér. 4, 21 no. 4 (1988), p. 653–663.
 [7] Daellenbach, H.G., C.A. De Kluyver (1980) . Note on multiple objective dynamic programming. Journal of the Operational Research Society 31 591–594.
 [8] De Loera, J.A., Haws, D., Hemmecke, R., Huggins, P., Tauzer, J., Yoshida, R. A User’s Guide for LattE v1.1. 2003, software package LattE is available at http://www.math.ucdavis.edu/latte/
 [9] De Loera, J.A, Haws, D., Hemmecke, R., Huggins, P., Sturmfels, B., and Yoshida, R. (2004). Short rational functions for toric algebra and applications. Journal of Symbolic Computation, Vol. 38, 2 , 2004, 959–973.
 [10] De Loera, J.A., Hemmecke, R., Köppe, M. (2008). Pareto Optima of Multicriteria Integer Linear Programs. To appear in: INFORMS Journal on Computing, 2008
 [11] Delorme, X. , Gandibleux, X. and Degoutin, F. (2003). Resolution approché du probleme de set packing biobjectifs. In Proceedings de l’ecole d’Automne de Recherche Operationnelle de Tours (EARO), 74–80.
 [12] Ehrgott, M. (2000). Approximation algorithms for combinatorial multicriteria optimization problems. International Transactions in Operational Research 7:5–31.
 [13] Ehrgott, M. and Gandibleux, X. (editors) (2002). Multiple Criteria Optimization. State of the Art Annotated Bibliographic Surveys. Boston, Kluwer.
 [14] Ehrgott, M. and Gandibleux, X. (2004). Approximative solution methods for multiobjective combinatorial optimization. TOP 12(1):1–88.
 [15] Ehrgott, M. and Gandibleux, X. (2007). Bound sets for biobjective combinatorial optimization problems. Comput. Oper. Res. 34, 9 (Sep. 2007), 26742694.
 [16] Fernández, E. and Puerto, J. (2003). The multiobjective solution of the uncapacitated plant location problem. European Journal of Operational Research. vol. 45, n.3 509529.
 [17] Garey, M. R. and Johnson, D. S. (1979). Computers and Intractability: a Guide to the Theory of Np.Completeness. W. H. Freeman & Co.
 [18] Ishibuchi, H. and Murata, T. (1998), A multiobjective genetic local search algorithm and its application to flowshop scheduling, IEEE Trans. Syst., Man, Cybern. C 28, 392–403.
 [19] Johnson, D. S. and Papadimitriou, C. H. (1988). On generating all maximal independent sets. Inf. Process. Lett. 27, 3 (Mar. 1988), 119–123.
 [20] Jozefowiez, N. , Semet, F. and Talbi, EG. (2004). A multiobjective evolutionary algorithm for the covering tour problem, Chapter 11 in ”Applications of multiobjective evolutionary algorithms”, C. A. Coello and G. B. Lamont (editors), p 247267, World Scientific.
 [21] Villarreal, B. and Karwan, M.H. (1982). Multicriteria Dynamic Programming with an Application to the Integer Case. Journal of Optimization Theory and Applications. Vol. 31. pp 4369.
 [22] Khovanskii, A.G. and Pukhlikov, A.V. , The RiemannRoch theorem for integrals and sums of quasipolynomials on virtual polytopes, (Russian) Algebra i Analiz 4 (1992), no. 4, 188–216; translation in St. Petersburg Mathematical Journal, 4 (1993), no. 4, 789–812.
 [23] Lasserre, J.B. Integer programming, Barvinok’s counting algorithm and Gomory relaxations. Operations Research Letters, 32, 2003, 133–137.
 [24] Lawrence,J. , Rationalfunctionvalued valuations on polyhedra, in: Discrete and Computational Geometry (New Brunswick, NJ, 1989/1990), 199–208, DIMACS Ser. Discrete Mathematics and Theoretical Computer Science, 6, American Mathematical Society, Providence, RI, 1991.
 [25] Lenstra, H.W. (Jr.) (1981). Integer programming with a fixed number of variables, Report 81–03, Mathematisch Instituut, Universiteit ban Amsterdam.
 [26] Marcotte, O., R.M. Soland (1986). An interactive branchandbound algorithm for multiple criteria optimization. Management Science 32 61–75.
 [27] Mavrotas, G., D. Diakoulaki (1998). A branch and bound algorithm for mixed zeroone multiple objective linear programming. European Journal of Operational Research 107 530–541.
 [28] Przybylski, A., Gandibleux, X. and Ehrgott, M. (2008). Two phase algorithms for the biobjective assignment problem. European Journal of Operational Research 185(2):509533, 2008.
 [29] ElSherbeny, N. (2001). Resolution of a Vehicle Routing Problem with Multiobjective Simulated Annealing Method, PhD thesis, Faculte Polytechnique de Mons, Belgium.
 [30] Steiner, S. and Radzik, T. 2008. Computing all efficient solutions of the biobjective minimum spanning tree problem. Comput. Oper. Res. 35, 1 (Jan. 2008), 198211.
 [31] Steuer, R.E. (1985). Multiple Criteria Optimization: Theory, Computation and Application. John Wiley & Sons, New York, NY.
 [32] Tsukiyama, S., Ide, M., Ariyoshi, H. and Shirakawa, I. (1979). A new algorithm for generating all maximal independent sets. SIAM J. Comput. 6, pp. 505–517.
 [33] Ulungu, E. and Teghem, J. (1995). The twophases method: An efficient procedure to solve biobjective combinatorial optimization problems, Foundations of Computing and Decision Sciences 20(2), 149–165.
 [34] Visée, M., Teghem, J., Pirlot, M., and Ulungu, E. L. (1998). Twophases Method and Branch and Bound Procedures to Solve the Biobjective Knapsack Problem. J. of Global Optimization 12, 2, 139–155.
 [35] Woods, K. and Yoshida, R. (2005). Short rational generating functions and their applications to integer programming , SIAG/OPT Views and News, 16 , 1519.
 [36] Zionts, S. (1979). A survey of multiple criteria integer programming methods. Annals of Discrete Mathematics 5, 389–398.
 [37] Zionts, S. and Wallenius, J. (1980). Identifying efficient vectors: some theory and computational results. Operations Research 23, 785–793.