A Gröbner bases methodology for solving multiobjective polynomial integer programs

# A Gröbner bases methodology for solving multiobjective polynomial integer programs

V. Blanco and J. Puerto Departamento de Estadística e Investigación Operativa
December 10, 2008
###### Abstract.

Multiobjective discrete programming is a well-known family of optimization problems with a large spectrum of applications. The linear case has been tackled by many authors during the last years. However, the polynomial case has not been deeply studied due to its theoretical and computational difficulties. This paper presents an algebraic approach for solving these problems. We propose a methodology based on transforming the polynomial optimization problem in the problem of solving one or more systems of polynomial equations and we use certain Gröbner bases to solve these systems. Different transformations give different methodologies that are analyzed and compared from a theoretical point of view and by some computational experiments via the algorithms that they induce.

This research was partially supported by Spanish research grant numbers MTM2007- 67433-C02-01 and P06-FQM-01366.

## 1. Introduction

A multiobjective polynomial program consists of a finite set of polynomial objective functions and a finite set of polynomial constraints (in inequality or equation form), and solving that problem means obtaining the set of minimal elements in the feasible region defined by the constraints with respect to the partial order induced by the objective functions.

Polynomial programs have a wide spectrum of applications. Examples of them are capital budgeting [17], capacity planning [6], optimization problems in graph theory [3], portfolio selection models with discrete features [2, 16] or chemical engineering [20], among many others. The reader is referred to [18] for further applications.

Polynomial programming generalizes linear and quadratic programming and can serve as a tool to model engineering applications that are expressed by polynomial equations. Even those problems with transcendental terms such as sin, log, and radicals can be reformulated by means of Taylor series as a polynomial program. A survey of the publications on general nonlinear integer programming can be found in [9].

We study here multiobjective polynomial integer programs (MOPIP). Thus, we assume that the feasible vectors have integer components and that there are more than one objective function to be optimized. This change makes single-objective and multiobjective problems to be treated in a totally different manner, since the concept of solution is not the same.

In this paper, we introduce a new methodology for solving general MOPIP based on the construction of reduced Gröbner bases of certain ideals related to the problem and on solving triangular systems of polynomial equations given by those bases. Gröbner bases were introduced by Bruno Buchberger in 1965 in his PhD Thesis [7]. He named it Gröbner basis paying tribute to his advisor Wolfgang Gröbner. This theory emerged as a generalization, from the one variable case to the multivariate polynomial case, of the Euclidean algorithm, Gaussian elimination and the Sylvester resultant. One of the outcomes of Gröbner Bases Theory was its application to linear integer programming [8, 14, 22]. Later, Blanco and Puerto [5] introduces a new notion of partial Gröbner basis for toric ideals in order to solve multiobjective linear integer programs. A different approach for solving linear integer programs was developed by Bertsimas et al. [4] based on the application of Gröbner bases for solving system of polynomial equations. This alternative use of Gröbner bases is also used in the paper by Hägglof et al. [13] for solving continuous polynomial optimization problems. Further details about Gröbner bases can be found in [10, 11].

We describe different approaches for solving MOPIP using Gröbner bases which are based on reducing the problem to several optimality conditions: the necessary Karush-Kuhn-Tucker, the Fritz-John and the multiobjective Fritz-John optimality conditions.

The paper is structured as follows. In the next section we give some preliminaries in multiobjective polynomial integer optimization. We present in Section 3 our first algorithm for solving MOPIP using only the triangularization property of lexicographic Gröbner bases. Section 4 is devoted to two different algorithms for solving MOPIP using a Chebyshev like scalarization and the Karush-Kuhn-Tucker or the Fritz-John optimality conditions. The last algorithm, based on the multiobjective Fritz-John optimality condition, is described in Section 5. Finally, in Section 6, we compare the algorithms with the results of some computational experiments and its analysis.

## 2. The multiobjective integer polynomial problem

The goal of this paper is to the solve multiobjective polynomial integer programs (MOPIP):

 (\it MOPIPf,g,h) min(f1(x),…,fk(x))s.t.gj(x)≤0j=1,…,mhr(x)=0r=1,…,sx∈Zn+

with polynomials in and the constraints defining a bounded feasible region. Therefore, from now on we deal with and we denote , and . If the problem had no equality (resp. inequality) constraints, we would denote it by (resp. ), avoiding the nonexistent term.

However, can be transformed to an equivalent multiobjective polynomial binary problem. Since the feasible region is assumed to be bounded, it can be always embedded in a hypercube . Then, every component in , , has an additional, but redundant, constraint . We write in binary form, introducing new binary variables with values in , , substituting every in we obtain an equivalent problem.

Then, from now on, without loss of generality, we restrict ourselves to multiobjective polynomial binary programs () in the form:

 (\it MOPBPf,g,h) min(f1(x),…,fk(x))s.t.gj(x)≤0j=1,…,mhr(x)=0r=1,…,sx∈{0,1}n

If the problem had no equality (resp. inequality) constraints, we would denote the problem by (resp. ), avoiding the nonexistent term.

The number of solutions of the above problem is finite, since the decision space is finite. Thus, the number of feasible solutions is, at most .

It is clear that 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 Pareto-optimal) points.

A feasible vector is said to be a nondominated (or Pareto optimal) solution of if there is no other feasible vector such that

 fj(y)≤fj(ˆx)∀j=1,…,k

with at least one strict inequality for some . If is a nondominated solution, the vector is called efficient.

We say that a feasible solution, , is dominated by a feasible solution, , if for all and . We denote by the set of all nondominated solutions for and by the image under the objective functions of , that is, . Note that is a subset of (decision space) and is a subset of (space of objectives).

From the objective functions , we obtain a partial order on as follows:

 x≺fy:⟺f(x)f(y) or x=y.

Note that since , the above relation is not complete. Hence, there may exist incomparable vectors.

In the following sections we describe some algorithms for solving MOPIP using tools from algebraic geometry. In particular, in each of these methods, we transform our problem in a certain system of polynomial equations, and we use Gröbner bases to solve it.

## 3. Obtaining nondominated solutions solving systems of polynomial equations

In this section we present the first approach for solving multiobjective polynomial integer programs using Gröbner bases. For this method, we transform the program in a system of polynomial equations that encodes the set of feasible solutions and its objective values. Solving that system in the objective values, and then, selecting the minimal ones in the partial componentwise order, allows us to obtain the associate feasible vectors, thus, the nondominated solutions.

Through this section we solve . Without loss of generality, we reduce the general problem to the problem without inequality constraints since we can transform inequality constraints to equality constrains as follows:

 (1) g(x)≤0⟺g(x)+z2=0,z∈R.

where the quadratic term, , assures the nonnegativity of the slack variable and then, less than or equal to type inequality. Initially, we suppose that all the variables are binary. In Remark 3.1 we describe how to modify the algorithm to incorporate the above slack variables.

This approach consists of transforming to an equivalent problem such that the objective functions are part of the constraints. For this transformation, we add new variables, to the problem, encoding the objective values for all feasible solutions. The modified problem is:

 (2) min(y1,…,yk)s.t.hr(x)=0r=1,…,syj−fj(x)=0j=1,…,kxi(xi−1)=0i=1,…,ny∈Rkx∈Rn

where integrality constraints are codified as quadratic constraints so, is a polynomial continuous problem.

The algorithm consists of, first, obtaining the set of feasible solutions of Problem (2) in the variables; then, selecting from that set those solutions that are minimal with respect to the componentwise order, obtaining the set of efficient solutions of . The feasible solutions in the -variables associated to those efficient solutions correspond with the nondominated solutions of .

Then, first, we concentrate in describing a procedure for solving the system of polynomial equations that encodes the feasible region of Problem (2), i.e. the solutions of

 (3) hr(x)=0 for all r=1,…,syj−fj(x)=0 for all j=1,…,kxi(xi−1)=0 for all i=1,…,n.

For analyzing the system (3) we use Gröbner bases as a tool for solving systems of polynomial equations. Further details can be found in the book by Sturmfels [21].

The set of solutions of (3) coincides with the affine variety of the following polynomial ideal in :

 I=⟨h1(x),…,hm(x),y1−f1(x),…,yk−fk(x),x1(x1−1),…,xn(xn−1)⟩.

Note that is a zero dimensional ideal since the number of solutions of the equations that define is finite. Let denote the affine variety of . If we restrict to the family of variables (resp. ) the variety (resp. ) encodes the set of feasible solutions (resp. the set of possible objective values) for that problem.

Applying the elimination property, the reduced Gröbner basis for , , with respect to the lexicographical ordering with gives us a method for solving system (3) sequantially, i.e., solving in one indeterminate at a time. Explicitly, the shape of is:

1. contains one polynomial in :

2. contains one or several polynomials in : .

3. contains one or several polynomials in : .

4. contains one or several polynomials in : .

Then, with this structure of , we can solve, in a first step, the system in the variables using those polynomial in that only involve this family of variables as follows: we first solve for in , obtaining the solutions: . Then, for fixed , we find the common roots of getting solutions and so on, until we have obtained the roots for . Note that at each step we only solve one-variable polynomial equations.

We denote by the above set of solutions in vector form

 Ω={(^y1,…,^yk):pk(^yk)=0,p1k−1(^yk−1,^yk)=0,…,pmk−1k−1(^yk−1,^yk)=0,…p11(^y1,^y2,…,^yk)=0,…,pm11(^y1,^y2,…,^yk)=0}.

As we stated above, is the set of all possible values of the objective functions at the feasible solutions of . We are looking for the nondominated solutions that are associated with the efficient solutions. From , we can select the efficient solutions as those that are minimal with respect to the componentwise order in . So, we can extract from the set of efficient solutions, :

 YE={(y∗1,…,y∗k)∈Ω:∄(y′1,…,y′k)∈Ω with y′j≤y∗j% for j=1,…,k and (y′1,…,y′k)≠(y∗1,…,y∗k)}

Once we have obtained the solutions in the variables that are efficient solutions for , we compute with an analogous procedure the nondominated solutions associated to the -values in . It consists of solving the triangular system given by for the polynomial where the -variables appear once the values for the -variables are fixed to be each of the vectors in .

A pseudocode for this procedure is described in Algorithm 1.

###### Theorem 3.1.

Algorithm 1 either provides all nondominated and efficient solutions or provides a certificate of infeasibility whenever .

###### Proof.

Suppose that . Then, has exactly one element, namely . This follows from the observation that is a polynomial ideal in one variable, and therefore, needs only one generator.

Solving we obtain every . Sequentially we obtain extending to the partial solutions in and so on.

By the Extension Theorem, this is always possible in our case.

Continuing in this way and applying the Extension Theorem, we can obtain all solutions in . These vectors are all the possible objective values for all feasible solutions of the problem. Selecting from those solutions that are not dominated in the componentwise order in , we obtain .

Following a similar scheme in the - variables, we have the set encoding all efficient (in the first coordinates) and nondominated (in the last coordinates) solutions.

Finally, if , then, the ideal coincides with , indicating that is empty (it is the set of the common roots of all polynomials in ). Then, we have an infeasible integer problem. ∎

###### Remark 3.1.

In the case when we have added slack variables, as explained in (1), we slightly modify the above algorithm solving first in the slack variables and selecting those solutions that are real numbers. Then continue with the procedure described in Algorithm 1.

###### Remark 3.2.

The Gröbner basis, , computed for solving the system of polynomial equations can be computed with respect to any other elimination ordering. The only conditions that are required for that ordering is that it allows to separate the family of variables from the family of -variables and such that the system of polynomials given by that basis allows solving first for the -variables and then for the -variables sequentially.

## 4. Obtaining nondominated solutions by the Chebyshev norm approach

In this section we describe two more methods for solving MOPIP based on a different rationale, namely scalarizing the multiobjective problem and solving it as a parametric single-objective problem. We propose a methodology based on the application of optimality conditions to a family of single-objective problems related to our original multiobjective problem. The methods consist of two main steps: a first step where the multiobjective problem is scalarized to a family of single-objective problems such that each nondominated solution is an optimal solution for at least one of the single-objective problems in that family; and a second step that consists of applying necessary optimality conditions to each one of the problems in the family, to obtain their optimal solutions. Those solutions are only candidates to be nondominated solutions of the multiobjective problem since we just use necessary conditions.

For the first step, the scalarization, we use a weighted Chebyshev norm approach. Other weighted sum approaches could be used to transform the multiobjective problem in a family of single-objective problems whose set of solutions contains the set of nondominated solutions of our problem. However, the Chebyshev approach seems to be rather adequate since it does not require to impose extra hypothesis to the problem. This approach can be improved for problems satisfying convexity conditions, where alternative well-known results can be applied (see [15] for further details).

For the second step, we use the Fritz-John and Karush-Kuhn-Tucker necessary optimality conditions, giving us two different approaches. In this section we describe both methodologies since each of them has its own advantages over the other.

For applying the Chebyshev norm scalarization, we use the following result that states how to transform our problem to a family of single-objective problems, and how to obtain nondominated solutions from the optimal solution of those single-objective problems. Further details and proofs of this result can be found in [15].

###### Theorem 4.1 (Corollary 11.21 in [15]).

Let be feasible. is a nondominated solution of if and only if there are positive real numbers such that is an image unique solution of the following weighted Chebyshev approximation problem:

 (Pω) minγs.t.ωi(fi(x)−^yi)−γ≤0i=1,…,kgj(x)≤0j=1,…,mhr(x)=0r=1,…,sxi(xi−1)=0i=1,…,nγ∈Rx∈Rn

where is a lower bound of , i.e., for all feasible solution and .

According to the above result, every nondominated solution of is the unique solution of () for some . We apply, in the second step, necessary optimality conditions for obtaining the optimal solutions for those problems (taking as parameters). These solutions are candidates to be nondominated solutions of our original problem. Actually, every nondominated solution is among those candidates.

In the following subsections we describe the above-mentioned two methodologies for obtaining the optimal solutions for the scalarized problems () for each .

### 4.1. The Chebyshev-Karush-Kuhn-Tucker approach:

The first optimality conditions that we apply are the Karush-Kuhn-Tucker (KKT) necessary optimality conditions, that were stated, for the general case, as follows (see e.g. [1] for further details):

###### Theorem 4.2 (KKT necessary conditions).

Consider the problem:

 (4) minf(x)s.t.gj(x)≤0=1,…,mhr(x)=0r=1,…,sx∈Rn

Let be a feasible solution, and let . Suppose that and , for , are differentiable at , that , for , is continuous at ,and that , for , is continuously differentiable at . Further suppose that , for , and , for , are linearly independent (regularity conditions). If solves Problem 4 locally, then there exist scalars , for , and , for , such that

 (KKT) ∇f(x∗)+m∑j=1λj∇gj(x∗)+s∑r=1μr∇hr(x∗)=0λjgj(x∗)=0for j=1,…,mλj≥0for j=1,…,m

From the above theorem the candidates to be optimal solutions for Problem (4) are those that either satisfy the KKT conditions (in the case where all the functions involved in Problem (4) are polynomials, this is a system of polynomial equations) or do not satisfy the regularity conditions. Note that these two sets are, in general, not disjoint.

Regularity conditions can also be formulated as a system of polynomial equations when the involved functions are all polynomials. Let be a feasible solution for Problem (4), does not verify the regularity conditions if there exist scalars , for , and , for , not all equal to zero, such that:

 (Non-Regularity) ∑j∈Iλj∇gj+s∑r=1μr∇hr=0

The above discussion justifies the following result.

###### Corollary 4.1.

Let be a nondominated solution for . Then, is a solution of the systems of polynomial equations (5) or (6), for some .

 (5) 1−k∑i=1νi=0k∑i=1νiωi∇fi(x)+m∑j=1λj∇gj(x)+s∑r=1μj∇hr(x)+n∑l=1βlδil(2xi−1)=0νiωi(fi(x)−^yi)−γ=0,i=1,…,k⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭

for some , for , , for , and , for .

 (6) k∑i=1νi=0k∑i=1νiωi∇fi(x)+m∑j=1λj∇gj(x)+s∑r=1μj∇hr(x)+n∑l=1βlδil(2xi−1)=0ωi(fi(x)−^yi)−γ≤0,i=1,…,k⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭

with such that , for , , for
, for , and , for .

Let denote the set of solutions, in the -variables, of system (5) and let denote the set of solutions, in the -variables, of system (6) (the problem is solved avoiding inequality constraints, then every solution is evaluated to check if it satisfies the inequality constraints).

For solving these systems (Chebyshev-KKT and Non-Regulariry), we use a Gröbner bases approach. Let be the ideal generated by the involved equations.

Let us consider a lexicographical order over the monomials in such that . Then, the Gröbner basis, , for with this order has the following triangular shape:

• contains one polynomial in :

• contains one or several polynomials in : .

• contains one or several polynomials in : .

• The remaining polynomials involve variables and at least one , , , or .

We are interested in finding only the values for the -variables, so, we avoid the polynomials in that involve any of the other auxiliary variables. In general, we are not able to discuss about the values of the parameters , , , and . Needless to say that in those cases when we can do it, some values of may be discarded simplifying the process. We denote by the subset of that contains only polynomials in the -variables. By the Extension Theorem, is a Gröbner basis for .

Solving the system given by and checking feasibility of those solutions, we obtain as solutions those of our KKT or Non-Regularity original systems.

It is clear that the set of nondominated solutions of our problem is a subset of , since either a solution is regular, and then, KKT conditions are applicable or it satisfies the non regularity conditions. However, the set may contain dominated solutions, so, at the end we must remove the dominated ones to get only .

The steps to solve Problem using the Chebyshev-KKT approach are summarized in Algorithm 2.

###### Theorem 4.3.

Algorithm 2 solves Problem in a finite number of steps.

### 4.2. The Chebyshev-Fritz-John approach

Analogously to the previous approach, once we have scalarized the original multiobjective problem to a family of single-objective problems, in this section we apply the Fritz-John (FJ) conditions to all the problems in this family. The following well-known result justifies the use of FJ conditions to obtain candidates to optimal solutions for single-objective problems. Proofs and further details can be found in [1].

###### Theorem 4.4 (FJ necessary conditions).

Consider the problem:

 (7) minf(x)s.t.gj(x)≤0=1,…,mhr(x)=0r=1,…,sx∈Rn

Let be a feasible solution, and let . Suppose that and , for , are differentiable at , and that , for , is continuously differentiable at . If locally solves Problem (7), then there exist scalars , for , and , for , such that

 (FJ) λ0∇f(x∗)+m∑j=1λj∇gj(x∗)+s∑r=1μr∇hr(x∗)=0λjgj(x∗)=0for j=1,…,mλj≥0for j=1,…,m(λ0,λ,μ)≠(0,0,0)

Note that, in the FJ conditions, regularity conditions are not required to set the result.

###### Corollary 4.2.

Let be a nondominated solution for . Then, is a solution of the system of polynomial equations (8) for some , for , , and .

 (8) λ0−k∑i=1νi=0k∑i=1νiωi∇fi(x)+m∑j=1λj∇gj(x)+s∑r=1μj∇hr(x)+n∑l=1βlδil(2xi−1)=0νi(ωi(fi(x)−^yi)−γ)=0,i=1,…,kλjgj(x)=0,j=0,…,m⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭

where , for , , for and , for , denotes the Kronecker delta function, .

Let denote the set of solutions, in the -variables, that are feasible solutions of and solutions of system (8).

The set of nondominated solutions of our problem is a subset of , since every nondominated solution is an optimal solution for some problem in the form of (8), and every solution of this single-objective problem is a solution of the FJ system.

However, dominated solutions may appear in the set of solutions of (8), so, a final elimination process is to be performed to select only the nondominated solutions.

The steps to solve using the Chebyshev-FJ approach are summarized in Algorithm 3.

###### Theorem 4.5.

Algorithm 3 solves in a finite number of steps.

The last part of the section is devoted to show how to solve the Chebyshev-FJ system using Gröbner bases.

Consider the following polynomial ideal

in the polynomial ring .

Let us consider a lexicographical order over the monomials in such that . Then, the Gröbner basis, , for with this order has the following triangular shape:

• contains one polynomial in :

• contains one or several polynomials in :

• contains one or several polynomials in :

• The remainder polynomials involve variables and at least one of , , , or .

We are interested in finding only the values for the -variables, so, we avoid the polynomials in that involve any of the other auxiliary variables. We denote by the subset of that contains only all the polynomials in the -variables. By the Extension Theorem, is a Gröbner basis for .

Solving the system given by , we obtain as solutions, those of our FJ original system.

###### Remark 4.1 (Convex Case).

In the special case where both objective functions and constraints are convex, sufficient KKT conditions can be applied. If the feasible solution satisfies KKT conditions, and all objective and constraints functions are convex, then is a nondominated solution. As a particular case, this situation is applicable to linear problems.

In this case, we may choose a linear scalarization instead of the Chebyshev scalarization. With this alternative approach, the scalarized problem is

 mink∑s=1tsfs(x)s.t.gj(x)≤0j=1,…,mhr(x)=0r=1,…,sxi(xi−1)=0i=1,…,n

for .

Then, by Corollary 11.19 in [15], and denoting by the feasible region, if is convex, then each is a nondominated solution if and only if is a solution of Problem 4.1 for some .

Using both results, necessary and sufficient conditions are given for that problem and the removing step is avoided.

###### Remark 4.2 (Single-Objective Case).

The same approach can be applied to solve single-objective problems. In this case, KKT (or FJ) conditions can be applied directly to the original problem, without scalarizations.

## 5. Obtaining nondominated solutions by multiobjective optimality conditions

In this section, we address the solution of by directly applying necessary conditions for multiobjective problems. With these conditions we do not need to scalarize the problem, as in the above section, avoiding some steps in the process followed in the previous sections.

The following result states the Fritz-John necessary optimality conditions for multiobjective problems.

###### Theorem 5.1 (Multiobjective FJ necessary conditions, Theorem 3.1.1. in [19]).

Consider the problem:

 (9) min(f1(x),…,fk(x))s.t.gj(x)≤0=1,…,mhr(x)=0r=1,…,sx∈Rn

Let a feasible solution. Suppose that , for , , for and , for , are continuously differentiable at . If is a nondominated solution for Problem 9, then there exist scalars , for , , for , and , for , such that

 (MO-FJ) k∑i=1νi∇fi(x∗)+m∑j=1λj∇gj(x∗)+s∑r=1μr∇hr(x∗)=0λjgj(x∗)=0for j=1,…,mλj≥0for j=1,…,mνi≥0for i=1,…,k(ν,λ,