Fast Computation of Pareto Set for Bicriteria Linear Programs with Application to a Diet Formulation Problem

In case of mathematical programming problems with conflicting criteria, the Pareto set is a useful tool for a decision maker. Based on the geometric properties of the Pareto set for a bicriteria linear programming problem, we present a simple and fast method to compute this set in the criterion space using only an elementary linear program solver. We illustrate the method by solving the pig diet formulation problem which takes into account not only the cost of the diet but also nitrogen or phosphorus excretions.


Introduction
Animal diet formulation is a very important problem from an economic and environmental point of view, so it is an interesting example in operations research.Many modern animal diet formulation methods tend to take into account nitrogen and phosphorus excretions that are detrimental from an environmental point of view.Following [1], it is appropriate to apply a tax on excretions so as to change the behavior of the producers in the swine industry.
These changes in behavior are studied using a formulation of the problem as a bicriteria problem and are obtained by the determination of the Pareto set of the problem.For linear models, this Pareto set is a simple polygonal line.This fact implies that changes in behavior of the producers are abrupt and correspond to particular values of the tax.In other words even in increasing the tax it can happen that there is no change in behavior.Behavior changes happend only at very particular values of the tax.We will see that these behaviors correspond to efficient extreme points of the Pareto set, and to every extreme point corresponds a tax interval so that any value of the tax in this interval leads to the behavior given by that extreme point.
The computation and visualization of the Pareto set, also known as the efficiency set, for bicriteria linear programming problems is a useful tool for decision makers.We could try to compute this set in the decision space [2]- [10], but due to the high dimension of this space, it can be a quite large and complicated set.Methods to obtain this set are also complicated, see for example [11].Fortunately, the geometric aspect of the Pareto set in the criterion (or outcome) space for bicriteria linear program is quite simple [12].
The outline of the paper is the following.The bicriteria problem is presented in Section 2. We will see in Section 3, that the Pareto set of a bicriteria linear problem is a simple polygonal line with L + 1 extreme points joined by L adjacent segments.Then in Section 4 we presents the link between the geometric structure of the Pareto set and the weighted-sums approach.Then an elementary algorithm to determine the Pareto set in the criterion space is suggested and its complexity is analyzed.Let us point out that this method uses only elementary result from a linear program solver, that is to say the optimal solution (values of the decision variables).This fact is an interesting property of the method.
Few methods exist for computing the Pareto set in the criteria space.One such method is presented in [13].The method requires information about the dual, assume the feasible set is compact, and determine the Pareto set with at most 2L + 4 calls to a linear program solver.Another simple method for bi-criteria problems is presented in [12] to obtain the Pareto set in the criterion space.The algorithm is based on information about the reduced costs of all nonbasic variables, which is equivalent to have information about the solution of the dual problem.For bi-criteria linear problems we could also use a parametric analysis to obtain the Pareto set [11] [14].The last two methods require that the software used to solve a linear program send information about the dual, reduced cost or postoptimal analysis, which is not always possible for a simple linear program solver.Unfortunately, even if it seems that those two methods require around 2L iterations, their complexities are nowhere analyzed.Moreover they can cycle as explained in [15] (pages 281-282), and [16] (pages 162-166).
Finally, in Section 5, we compute the Pareto sets for least cost diet formulation problems for pig, or any monogastric animal, taking into account the nitrogen and/or phosphorus excretions.Tax systems related to efficient extreme points of this problem are described.

Bicriteria Linear Programming Problem
Let us consider the standard form of the bicriteria linear programming problem [11] F. Dubeau where x is a column vector in n  , and the k c 's ( 1, 2 k = ) are two row vectors , , , where A is a ( ) , m n -matrix, and b is a column vector in m  .Let C be the ( ) It is well-known that  and c  are polyhedral sets in n  and 2  respectively.Throughout this paper we will suppose that the two criteria are lower bounded on  which means that for 1, 2 i = we have

Efficiency Set
A feasible solution x ∈  is an efficient solution if and only if it does not exist any other feasible solution x ∈  such that 1) ( ) ( ) ( ) ( ) . The set of all efficient solutions is called the efficiency set noted  , also called Pareto set.The corresponding set in the criterion space is the set c C =  .

Geometric Structure
Under the assumption that the two cost vectors and we consider the single criteria problem for From [11] we have Hence the efficiency set  in the decision space is a connected set and is the union of faces, edges and vertices of  .This set may be quite complex due to the high dimension of the decision space.On the other side c  , which is the image in 2  of  by a linear transform, is a much simpler set.
Since we have assumed that both criteria are lower bounded on  , it follows that c  is a simple compact polygonal line.Indeed in that case c  is the union of a finite number L of segments [ ] To each segment is associated a weight 1, l l λ − such that the vector ( ) , ,

Weak Efficiency Set
We will call weak efficiency set, or weak Pareto set, the set defined by .
In the criteria space we will have Geometrically in the criterion space 2  , this means we add to c  possibly a vertical segment or a ray from 0 Q in the positive direction of 2 z , ( ) and/or a horizontal segment or a ray from L Q in the positive direction of 1 z , ( ) where 0 η and L η are nonnegative finite or infinite values.They are the maximal values of η such that ( )

Link to Parametric Analysis
The parametric analysis is based on the weighted-sum given by ( ) ( ) ( ) µ ∈ +∞ , and the value function in this case is defined by P λ , we could consider the single criteria problem for ) ( ) , , In many applications, the parameter µ is in fact a tax over the the second criteria (for a minimization problem).Interesting enough is to observe that the behavior change (extreme point) only for the critical values 1, l l µ − of the parameter µ .Indeed when µ increases and its value passes through 1, l l µ − , the optimal point, extreme point, move from

Preliminaries
Let us associate to any ( ) Hence we have the following results.
Theorem 4.1.[12] Let ( ) Theorem 4.3.[17] The function ( ) ϕ λ is continuous, piecewise linear and concave.The abscissae of slope changes are the increasing values Let us observe that the slope associated to We deduce the next results.
Theorem 4.4.[12] Let i Q′ and j Q′ be two distinct points on c Theorem 4.5.Let i Q′ and j Q′ be two distinct points on c  and , and consider the following two lines , the point of intersection of ( ) ( ) ; , ,

Algorithm
In this section we consider both criteria upper bounded on  .In the forthcoming algorithm we initialize the process with the two points 0 Q and  .Then we gradually obtain a sequence of points { } 0 a sequence of intervals associated to these points At the end of the process I L = and we have . We get the initial point ( ) which as the same first coordinate as 0 Q , and which as the same second coordinate as L Q .Those two points might not be on c  , but are STEP 1.As long that there exists an index i such that such index i and do: (C) Update the list of points { } 0 We modify as follow: a) for In the sequel no more point will be generated on [ ] , insert the point and modify intervals as follows (see Theorem , hence we modify , hence we modify For any i such that 0 Let us observe that this process use only optimal solutions of ( ) optimal values of the decision variables, which is easily obtained from any elementary linear program solver.
Remark 4.8.This algorithm produces at each iteration an inner and an outer approximation.The inner approximation is the polygonal line joining the i Q′ for 0, , i I =  .The outer approximation is the polygonal line joining the points 0 Q′ , ( ) , , Q′ , as long as the ( ) 's are well determined (see Theorem 4.7).At the end of the algorithm the two approximations agree.

Complexity
In this section we are going to determine the maximum number of calls to a linear program solver to completely determine the Pareto set, or equivalently its 1 The result is given in the last theorem of this section and says that it takes at most 2 3 L + calls to a linear program solver to generates the 1 We will use the following ordering on Theorem 4.9.The algorithm generates at most 3 points on [ ] and two of these points are Proof.Let us remark that the algorithm will eventually find a point in This first point can be generated at STEP 0, an initial point, if Otherwise, it is generated through STEP 1-C-II, with Then this point is included in the list, and there are three cases to study:  and we have ( ) , ). 2) ,  and we have ( ) and we have ( ) Q is generated through STEP 1-C-II, and it is included in the list.There are two cases to study: 1)  with ** λ modified as in the preceding case.Moreover if ( ) we will modify the upper bound to get ( ) 2) [ )  with ** λ modified as in the preceding case.Moreover if ( ) we will modify the lower bound to get ( ) Otherwise the two points are the extreme points There are two cases to study: 1) We have only one extreme point the segment in the list and ( ) As in the preceding paragraph, we will introduce it in the list, and depending of the case, by passing through STEP 1-C-II, 3 . Moreover, American Journal of Operations Research we will have respectively ( ) ( ) 2) We already have two extreme points 1 l Q − and l Q in the list.In that case We pass through STEP 1-C-I and we modify the intervals to get ( ) ( ) is not added to the list.
In the sequel, the algorithm generate no more point on [ ] we have two points or else, if we have three points, ( ) Theorem 4.10.
is eventually removed of the list without any supplementary call to the linear program solver.
Proof.When 0 Q is introduced in the list, there is no supplementary call for . Similarly for the interval ,

A Real World Application: Pig Diet Formulation
To illustrate our method of computation of the Pareto set we consider the pig diet formulation problem taking into account not only the cost of the diet but also environmental considerations, such as the reduction of nitrogen or phosphorus excretions.One way to analyze this problem is to rewrite the problem as bicriteria problem.Hence the Pareto set indicates the effect of the reduction of excretions, nitrogen or phosphorus, on the cost of the diet.This information is certainly useful for a decision maker which have to choose a diet which decrease the excretions without being too expensive [1].Even if in thispaper we describe the problem for the swine industry, the method could be applied to any monogastric animal: pig, rabbit, chicken, etc.

Classical Model
The least cost diet problem, introduced in [18], is a classical linear programming problem [19] [20] [21].A decision variable j x is assigned to each ingredient and represents the amount (in kg) of the j th ingredient per unit weight (1 kg) of the feed.Together, they form the decision vector The constraints impose some bounds on the quantity of the different ingredients in the For example a unit of feed is produced (a 1 kg mix), expressed by the constraint is the amount of digestible phosphorus contained in a unit of ingredient j.

Modelling of Nitrogen and Phosphorus Excretions
Nitrogen and phosphorus excretions are directly related to the excess of amounts of protein (amino acids) and phosphorus in the diet.Hence, we have to establish the protein and the phosphorus contents of the diet and take into account the parts that are actually assimilated.
The protein content of a diet , where j pr is the amount of protein per unit of ingredient j.The total excretion of protein ( ) pr r x is then given by the amount in protein of the diet from which we remove the amount of protein effectively digested given by As for the nitrogen, the amount of phosphorus of a unit weight diet , where j ph is the amount of phosphorus per unit of ingredient j.The amount * ph b is the the amount of phosphorus which is actually digested.In this way the phosphorus excretion ( ) ph r x is given by the phosphorus content of the diet from which we remove the amount of which is actually digested Hence, decreasing the phosphorus excretion ( ) ph r x is equivalent to decreasing the phosphorus content pr q x of the diet while maintained fixed the needs * ph b in phosphorus.

Data
The ingredients and their corresponding variables are described in Table 1.
Table 2 contains the entire model together with the values of the technical coefficients of the model.

Software
The algorithm was programmed in MATLAB, which includes in its standard library the linear program solver called Linprog.This software can use the simplex method or an interior point method.

Two Criteria Models and Results
At first we analyse the relation between the cost of the diet and the two different excretions (nitrogen and phosphorus).As a curiosity, we also consider the interactions between the two kind of excretions: nitrogen and phosphorus.

Cost and Excretions
We have considered two separate bicriteria models.We look for least cost diets while taking into account the nitrogen excretion for the first model and the phosphorus excretion for the second model.For each of these two bicriteria problems, the Pareto curve indicates the diet cost increase caused by an excretion decrease.
While considering the nitrogen excretion, the problem is : ( ) From its associated weighted-sum model given by subject to ,  , and the change in the behavior will happend only when the taxation level µ passes through the extremities l µ or l µ of this interval A similar analysis holds for the second bicriteria problem with phosphorus excretion.Indeed, for the phosphorus excretion problem, the model is: From its associated weighted-sum model given by ( ) ( )  shows the opposite effect of trying to reduce simultaneously both excretions.
Minimizing one excretion leads to an increse in the other excretion.For this problem, the algorithm detects 5 L = segments and 6 extreme points.A total of 12 calls to the linear program solver was required (the predicted maximum is 2 3 13 L + = ).
For each 0, ,5 l =  , the value function is   Let us observe that the last line of Table 3 ( 10 l = ) corresponds to the first line of Table 5 ( 0 l = ) and the last line of Table 4 ( 22 l = ) corresponds to the last line of Table 5 ( 5 l = ).

Conclusion
In this paper we have considered bicriteria linear programming problems and have presented an elementary and efficient algorithm to compute the Pareto set in the criterion space.We have illustrated the method on a real important application.This application also suggests that it could be interresting to extend the method to three-criteria problems.Moreover it could be interesting to compare our method to other methods to find the Pareto set in the criterion space, but it is out of the scope of this paper and could be a nice subject for a future research.

(
λ and µ are related by the formulae

STEP 3 .
End of the process (and I L = ).The output is the list

Theorem 4 . 11 .
The algorithm generates the extreme points { } 0L l l Q = of the Pareto set in at most 2 3L + calls to a linear program solver.Proof.The initialization STEP 0 requires 2 calls.For STEP 1, as we generate there is at most 2 1 L + calls.Hence the algorithm requires at most 2 3 L + calls.

==
in our model.The model's objective function is the diet cost.A vector of unit costs is American Journal of Operations Research used, where each j c represents the unit cost of the j th ingredient (euro/kg or $/kg).Thus the total cost of a unit of weight (1 kg) of diet which must be minimized over the set of feasible diets denoted by  .The classic least cost animal diet formulation model is:

.
Some ingredients, or combinations of ingredients, can be imposed on the diet.These restrictions give rise to equality constraints (=) or inequality constraints (≥ or ≤).More specifically, to satisfy protein requirements, the following constraints are introduced for the L groups of amino acids contained in the ingredients.We set() of digestible amino acid l contained in a unit of ingredient j and * l b is the minimum amount of digestible amino acid l required.Finally, the diet must satisfy the digestible phosphorus requirements equivalent to decrease the protein content pr q x of the diet while maintained fixed the needs * pr b in protein.
the set of efficient extreme points of the Pareto set in the criterion space, and the Pareto curve is sketched in Figure 1.For this problem, the algorithm detects 10 L = segments and 11 efficient extreme points .A total of 22 calls to the linear program solver was required (the predicted maximum is 2 3 23 L + = ).
all value of the parameter λ in the interval , in the behavior will happend for values of the parameter λ corresponding to the extremities l λ ou l λ of this interval.
, M. E. N. Habingabwa DOI: 10.4236/ajor.2018.85019325 American Journal of Operations Research one of the two points used to generate2

Table 1 .
List of available ingredients.

Table 3 .
Efficient extreme points in the criterion space

Table 4
presents the efficient extreme points in the criterion space while the Pareto curve is sketched in Figure2.For this problem, the algorithm detects

Table 4 .
Efficient extreme points in the criterion space