Regional division and reduction algorithm for minimizing the sum of linear fractional functions

This paper presents a practicable regional division and cut algorithm for minimizing the sum of linear fractional functions over a polyhedron. In the algorithm, by using an equivalent problem (P) of the original problem, the proposed division operation generalizes the usual standard bisection, and the deleting and reduction operations can cut away a large part of the current investigated region in which the global optimal solution of (P) does not exist. The main computation involves solving a sequence of univariate equations with strict monotonicity. The proposed algorithm is convergent to the global minimum through the successive refinement of the solutions of a series of univariate equations. Numerical results are given to show the feasibility and effectiveness of the proposed algorithm.


Introduction
Consider the following class of fractional programming: where A = (a rj ) m×n is a real matrix, c i = (c ij ) 1×n and d i = (d ij ) 1×n are real vectors for each i = 1, . . . , N , b = (b r ) m×1 , c 0i , d 0i ∈ R. In addition, we assume that D = y ∈ R n | Ay ≤ b, y ≥ 0 is a nonempty compact set. Since by choosing a sufficiently large valueM i (i = 1, . . . , N ) with c i y + c 0i +M i (d i y + d 0i ) > 0 for any y ∈D. Therefore, in the following, without loss of generality, we suppose that c i y+c 0i > 0 and d i y + d 0i > 0, ∀y ∈D, for each i = 1, . . . , N . Problem (FP) is a well-known class among fractional programming problems. Theoretically, it is NP-hard [1,2]. The primary challenges in solving problem (FP) arise from a lack of useful properties (convexity or otherwise) and from the number of ratios. In general, problem (FP) possesses more local optimal solutions that are not globally optimal [3], and so problem (FP) owns major theoretical and computational difficulties. From an application point view, this problem has a large deal of applications; for instance, traffic and economic domain [4,5], multistage stochastic shipping problems [6], data envelopment analysis [7], and queueing-location problems [8]. The reader is referred to a survey to find many other applications [4,5,[9][10][11].
Many algorithms have been proposed to solve problem (FP) with a limited number of ratios [4,[12][13][14][15][16][17][18][19]. For instance, Wang and Shen [4] give an efficient branch and bound algorithm by using a transformation technique and the linear relaxation programming of the objective function. By applying Lagrangian duality theory, Benson [9] presents a simplicial branch and bound duality-bounds algorithm. Carlsson and Shi [10] propose a linear relaxation algorithm with lower dimension which is performed on a 2N -dimensional domain instead of the original n-dimensional one, the computational time is long with larger N . Ji and Zhang [16] consider a deterministic global optimization algorithm by utilizing a transformation technique and a linearizing method. Jiao et al. [17,18] present the branch and bound algorithms for globally solving sum of linear and generalized polynomial ratios problems, by solving a sequence of linear relaxation programming problems. In short, most of them (see [4,9,[16][17][18] for example) are branch and bound algorithms.
The key idea behind such algorithms mentioned above is that the branch and bound operator is performed on an N -dimensional region (or the correspondence relating to N and n) rather than the native n-dimensional feasible set, that is, they all work on a space whose dimension increases with the number N of ratios.
In this article, a new division and reduction algorithm is proposed for globally solving problem (FP). To solve problem (FP), an equivalent optimization problem (P), whose objective function is just a simple univariate, is first presented by exploiting the feature of this problem. Then, in order to design a more efficient algorithm, several basic operations: division, deleting, and reduction, are incorporated into a similar branch and bound framework by utilizing the particular structure of problem (P). Compared with the usual branch and bound (BB) methods (e.g., [4,5,9,10,16]) mentioned above, the goal of this research is three fold. First, the proposed bounding operation is simple since the lower bound of the subproblem of each node can be achieved easily only by arithmetic computations, distinguishing it from the ones obtained by solving convex/linear programs in the usual BB methods. Second, the reduction operation that does not appear in other BB methods is used to tighten the range of each variable, such that the growth of the branch tree can be suppressed. Moreover, the main computational cost of the algorithm is to implement the reduction operation which solves univariate equations with strict monotonicity. Third, the problem in this paper is more general than the others considered in [10,11], since we only request d i y + d 0i = 0 for each i. Further, the proposed adaptive division operation both generalizes and is superior to the usual standard bisection in BB methods according to the numerical computational result in Sect. 5. Also, the computational results of the problem with the large number of ratio terms can be obtained to illustrate the feasibility and validity of the proposed algorithm.
This paper is summarized as follows. In Sect. 2, by using a conversion strategies, the original problem (FP) is transformed into an equivalent optimization problem (P). The adaptive division, deleting, and reduction operations are shown in Sect. 3. In Sect. 4 we give the proposed algorithm and its convergence. Some numerical results demonstrate the feasibility and availability of the algorithm in Sect. 5.

Equivalent problem
For solving problem (FP), we first convert the primary problem (FP) into an equivalent optimization problem (P), in which the objective function is a single variable and the constraint functions are the difference of two increasing functions. To see that such a reformulation is possible, let us denote, for each i = 1, . . . , N , Clearly, L i , U i > 0 for each i. Additionally, by introducing some additional variables w = (w 1 , w 2 , . . . , w N ) ∈ R N , from (FP) we then obtain the following equivalent problem: For simplicity, we denote Define a box as follows: [y,ȳ] = y ∈ R n | y l j min y∈D y j ≤ y j ≤ y u j max y∈D y j , j = 1, . . . , n .
Based on the above notations, by introducing an extra variable y 0 ∈ R + , we can convert the above problem (FP1) into the form: Then problems (FP2) and (FP) are equivalent in the sense of the following result.

Proposition 2.1 y * ∈ R n is a globally optimal solution for problem (FP) if and only if
(y * 0 , y * , w * ) ∈ R n+N+1 with y * 0 ∈ R and w * ∈ R N is a globally optimal solution for problem (FP2), where y * 0 = N i=1 w * i (c i y * + c 0i ) and w * i = 1 d i (y * )+d 0i for each i = 1, . . . , N .
Proof The proof of this result is obvious.
From Proposition 2.1, notice that, in order to globally solve problem (FP), we may globally solve problem (FP2) instead. Moreover, it is easy to see that the global optimal values of problems (FP) and (FP2) are equal.
In addition, for convenience of the following discussion, let us denote x = (y 0 , y, w) ∈ R n+N+1 with y ∈ R n , w ∈ R N . Then problem (FP2) can be rewritten as problem (P): (P) :

being increasing functions given by
Note that both problem (FP) and problem (P) are equivalent according to Proposition 2.1, hence for globally solving problem (FP), the algorithm to be presented concentrates on how to solve problem (P).

Essential operations
For solving problem (P), we first give the following concept about an approximate optimal solution. Definition 3.1 For given ε, η ≥ 0, let a point x ∈ D ε is said to be an ε-feasible solution to problem (P). If an ε-feasible solution x = (x 0 ,x 1 , . . . ,x n+N ) to problem (P) satisfies x is called an (ε, η)-optimal solution to problem (P).
For seeking an (ε, η)-optimal solution of problem (P), a division and cut algorithm to be developed includes three essential operations: division operation, deleting operation, and reduction operation.
First, the division operation consists in a sequential box division of the original box D 0 following in an exhaustive subsection principle, such that any infinite nested sequence of division sets generated through the algorithm reduces to a singleton. This paper takes an adaptive division operation, which extends the standard bisection in the normally used exhaustive subsection principle. Second, by using overestimation of the constraints, the deleting operation consists in eliminating each subbox D generated by the division operation, in which there is no feasible solution. In addition, the reduction operation is used to reduce the size of the current partition set (referred to a node), aiming at tightening each subbox which contains the feasible portion currently still of interest. For for convenience, we will use the following functions throughout this paper: where α, β ∈ (0, 1), p i ≤ p i ≤ q i , and e i denotes the ith unit vector, i.e., a vector with 1 at the ith position and 0 everywhere else, i = 0, 1, . . . , n + N . At a given stage of the proposed algorithm for problem (P), let V represent the best current objective function value to problem (P). Next, we will show these detailed operations.

Deleting operation
In this subsection, we will give a suitable deleting operation, which offers a possibility to remove a subbox D of D 0 without feasibility. Toward this end, we take into account a subproblem of problem (P) over a given box Given η > 0, for solving P(D), we need to seek out a feasible solutionx = (x 0 ,x 1 , . . . , x n+N ) ∈ D of P(D) such thatx 0 ≤ Vη, or to draw a conclusion that there exists no suchx, where V is the best objective value to problem P(D) known so far. Obviously, if p 0 > Vη, there exists nox ∈ D withx 0 ≤ Vη. If p 0 ≤ Vη and G k (p) ≤ 0 for each k = 0, . . . , m + N , then updatex = p. Therefore, without loss of generality, in the following discussion we shall suppose that Clearly, ifF(x) > ε for any x ∈ D, there exists no feasible solution on D, and so D is deleted for further discussion. However, since the judgement satisfyingF(x) > ε is not easy, we introduce an auxiliary problem of P(D) as follows: Observe that the objective and constraint functions are interchanged in P(D) and Q(D). Let V (P(D)) and V (Q(D)) be the optimal values of problems P(D) and Q(D), respectively. Then we give the following results about problems P(D) and Q(D). Theorem 3.1 Let ε > 0, η > 0 be given, and let D be a box with D ⊆ D 0 . We have the following result: Proof (i) This result is obvious, and here it is omitted.
(ii) By utilizing the assumption V (Q(D 0 )) > 0, i.e., we have the following conclusions: Consequently,x is an (ε, η)-optimal solution of P(D), and this accomplishes the proof. Theorem 3.1 illustrates that by utilizing problem Q(D) one can know whether or not there exists a feasible solutionx = (x 0 ,x 1 , . . . ,x n+N ) of P(D) with improving the current objective function value V , i.e.,x 0 ≤ Vη. Further, if V (Q(D 0 )) ≥ ε > 0, then an (ε, η)optimal solution of P(D 0 ) can be obtained, or it can be confirmed that problem P(D 0 ) has no feasible solution.
Additionally, if V (Q(D)) > ε, it is easy to see that a fathomed box D cannot contain the feasible solutions for problem (P). Thus, D is excluded from further consideration by the search. Unfortunately, solving problem Q(D) may be as difficult as solving the original problem P(D). Hence, a lower bound LB(D) of V (Q(D)) is required for eliminating the part of D 0 which does not contain the solution to problem (P). Clearly, if LB(D) ≥ ε, D is excluded from further consideration by the search. Since G + k (x) and Gk (x) are all increasing, an apparent lower bound is Although quite straightforward, the bound is sufficient to confirm convergence of the algorithm, as we will see immediately. For strengthening the computational validity in solving problem P(D), some tighter lower bounds can be obtained by utilizing the following Theorem 3.2. Especially, what is more important is that such tighter lower bounds can be gained only by simple arithmetic computations, which is different from the ones in the usual branch and bound methods for solving convex or linear programs [5,6,11,12,14].
Then a lower bound LB(D) satisfying LB(D) ≤ V (Q(D)) for problem Q(D) is given by Proof (i) If α = 0, this result is obvious.
Notice that LB(D) in Theorem 3.2 satisfies

Adaptive division
The division operation repeatedly subdivides an (n + N + 1)-dimensional box D 0 into (n + N + 1)-dimensional subboxes. This operation helps the algorithm confirm the position of a global optimal solution in D 0 for problem (P). Throughout this algorithm, we take a new adaptive subdivision principle as follows.
Adaptive subdivision For given η > 0, consider (iii) By usingx t , let us denote Based on the above division operation, D is divided into two new boxes D 1 and D 2 . Especially, when α = 0, the adaptive subdivision simply reduces to the standard bisection. As we will see from numerical experiments in Sect. 5, the adaptive subdivision is superior to the standard bisection. Moreover, the subdivision can confirm the convergence of the algorithm, and we have the following results.
According to Ref. [18], this adaptive division ensures the existence of an infinite sub- Thus we can obtain that In addition, by assumption (3.1), since u s t 0 = p s t 0 ≤ Vη, it holds that which implies thatû is feasible to Q(D 0 ), then we get Consequently, the limitation LB(D s t ) → V (Q(D 0 )) (t → +∞) holds and then we accomplish the proof.

Reduction operation
At a given stage of the proposed algorithm for problem (P), let D = [p, q] ⊆ D 0 be a box generated by the division operation and still of interest. Clearly, the smaller this box D is, the closer the feasible solution will be to the (ε, η)-optimal solution to problem (P). Hence, to effectively tighten the variable bounds in a particular node, a valid range reduction strategy is introduced by overestimation of the constraints, and by applying the monotonic decompositions to problem (P). Based on the above discussion, for any box D = [p, q] ⊆ D 0 generated by the division operation and still of interest, we intend to identify whether or not the box D contains a feasible solutionx of P(D) such thatx 0 ≤ Vη. Consequently, seeking such a pointx can be confined to the setD ∩ D, wherê The reduction operation is based on special cuts that exploit the monotonic structure of problem (P), and it aims at substituting the box D = [p, q] with a smaller box D without losing any valid point x ∈D ∩ D, i.e.,  Proof (i) The proof of this result is easy, here it is omitted.
(ii) The former of the conclusion is apparent, we only need to give the proof of the latter. If there exists some i ∈ {0, 1, . . . , n + N} so that max{f i k (0)|k = 0, 1, . . . , m + N} > ε, we can obtain Therefore, we can acquire R[p, q] = ∅, and this accomplishes the proof.  Proof For any given x = (x 0 , . . . , x n+N ) T ∈ [p, q], we first confirm that x ≥p.
Assume that x p, that is to say, there exists some i ∈ {0, 1, . . . , n + N} such that Then we discuss as follows: Ifα i ∈ (0, 1), from Rule (i) and the definition ofα i , there must exist some k such that f i k (α i ) = ε, i.e., In addition, through Rule (i) and the definition of Gk (x), it holds from (3.6) and (3.7) that Consequently, According to the above discussions, we have demonstrated that x ≥p, i.e., x ∈ [p, q]. Next we will show that x ≤q.

From Theorem 3.5, by Rules (i) and (ii) the main computational effort for deriving R[p, q]
is to solve some univariate equations about the variablesα i k andβ i k , which is easy to solve, for example, by using the bisection approach. What is more, as seen below, the cost of main computation in the proposed algorithm is also to form R[p, q].

Algorithm and its convergence
According to the above discussions, the proposed algorithm is shown as follows.

Algorithm statement
Step ( Step (2) DenoteM t to be the collection of boxes that results from M t after accomplishment of Step (1), and setN t = N t ∪M t .
(a) IfN t = ∅, stop: x * is an (ε, η)-optimal solution of problem (P). (b) IfN t = ∅, Theorem 3.4 is applied to each D ∈N t , then eliminate D and update x * = p, V = p 0 if necessary.
Step (3) DenoteN t to be the collection of boxes after accomplishment of Step (2). Choose D t = [p t , q t ] ∈ argmin{LB(D) | D ∈N t }, and denote LB t = LB(N t ). If LB t > 0, then terminate: the conclusion is the same as situation (a) of Step (2); otherwise go to Step (4). Step Step (5) Divide D t into two subboxes by the adaptive division, and set M t+1 to be the collection of these two subboxes of D t . Let N t+1 =N t \{D t }. Set t = t + 1, and return to Step (1).
Remark By utilizing a local solver in Step (4) of the proposed algorithm, instead of evaluating one point, we may acquire a point with the better objective function value V , and the iteration count of the algorithm may be reduced. However, because the computational cost increases rapidly with the quality of the objective value V , the running time is not always decreasing. So a trade-off must be made, practically just one evaluating point as in the above algorithm is used. Moreover, to implement the algorithm, all that is required is the ability to solve univariate equations with monotonicity and to execute simple algebraic steps.
The convergence of the algorithm is shown by the following theorem. Theorem 4.1 For given tolerances ε, η > 0, the above algorithm always terminates after finitely many iterations, obtaining an (ε, η)-optimal solution of problem (P), or a demonstration that the problem is infeasible.
Proof Since any feasible solution x = (x 0 , x 1 , . . . , x n+N ) to problem (P) with x 0 ≤ Vη must exist in some box D ∈N t , case (a) of Step (2) implies that there exists no such solution x. In addition, if LB t > 0 occurs in Step (3), then we can obtain V (MQ(D)) > 0. Consequently, by Theorem 3.5 the conclusion is true if one of the following cases occurs: that is to say, for extremely large t Steps (4) and (5) cannot appear. It remains to illustrate that eitherN t = ∅ or LB t > 0 must take place for extremely large t.
By contradiction, assume that the algorithm is infinite. Because each appearance of Step (4) reduces the current best objective function value V = x * 0 at least by η > 0 ifF(s t ) < ε, this conflicts with the fact that x 0 is bounded below if Step (4) takes place infinitely. In other words,F(s t ) < ε in Step (4) cannot occur infinitely often. Therefore, for all t extremely large, x * is unaltered, andF(s t ) ≥ ε while LB t ≤ 0, and then Step (5) must be performed infinitely. By the division operation of the algorithm, as t → +∞, q tp t → 0, we have According to (3.5), it holds that thus, it follows that which conflicts with lim t→+∞F (s t ) =F(x) ≥ ε > 0. Consequently, the algorithm must be finite, and the proof is finished.
Next, we are particularly interested in instances of problem (FP) with a small number of variables and a lager number of ratios. The reason is that this case has many applications, especially in layered manufacturing and material layout [20][21][22]. In the following, we will show computational results of experimenting on the proposed algorithm for randomly generated problems, which are of the following form: where the elements of the matrix A ∈ R m×n , c i , d i ∈ R n (i = 1, . . . , N ) are randomly generated in the unit interval [0, 1] and c ∈ R. As for the subdivision operation, except for the proposed operation, we also test the usual bisection operation in branch and bound methods (e.g., [4,5,[16][17][18]23]) for comparison, where their respective computer programs were referred to as adaptive division and bisection. During running procedures, η is fixed on 10 -5 , then we obtain a few data up and down; ultimately we choose the average result for running 10 times.
The comparison results between the adaptive division and bisection are shown in Figs. 1-4. Figure 1 illustrates the variedness of average computational time (in seconds) acquired by each item with c changing from 10 to 100 in 10 increments when n = 3, N = 50, and ε = 0.05, by which one can see that the influence of the value of c on computational time is slight, hence we let c be constant 3 in other experiments. Figure 2 gives the alteration of the average computational time (in seconds) gained by each term with N ranging from 50 to 400 when (n, ε) = (3, 0.05), by which one can observe that when N > 150 the running time increases rapidly, so the effect of variedness of the number of ratios on computational time is biggish. Figure 3 demonstrates the change of the average computational time acquired by each item with ε ranging from 0.01 to 0.07 in 0.01 increments. It is obvious that this program is sensitive to changes in ε. From numerical results, we can see that the adaptive division acquires less time required by bisection for each program.
For the adaptive division and bisection the computational time of solving examples with larger N ranging from 600 to 1500 are listed in Table 2. We can see that even if it takes about seventeen minutes to solve problem (P) of size (n, N) = (3, 1500), the adaptive division takes less time required by bisection for each N , which confirms the feasibility and availability of the proposed algorithm. We can draw a conclusion that the algorithm has more than enough performance, at least for three dimensions.
How does the algorithm behave for instances with n > 3? Unfortunately, the performance of the adaptive division rapidly deteriorates with increasing n, as shown in Fig. 3. For the same set of instances as in Fig. 1, Fig. 4 shows the variation of average computational time obtained by each program with n varying from 10 to 100 when (N, ε) is fixed at (50, 0.05), from which one can know that even though the running time is large, but the impact of variedness of the number of variables is mild with increasing running time in a nearly linear mode with increasing N .
Based on the above results, it is easy to see that the adaptive division has many advantages over the usual bisection in the computation. Additionally, when n is not larger than 3, the algorithm can rapidly solve problems, even if N takes on values as high as 1500. On the other hand, when the number of variables is as high as 100, with increasing n, we need to solve more equations in reduction operation, therefore, it is reasonable that the computational time is larger. Moreover, the lower bound is acquired only by simple arithmetic computations, that is, by executing simple algebraic steps, which is different from the ones used in solving convex or linear programs in common branch and bound methods. Finally, the computation for obtainingα i k andβ i k is easy by solving the equations with univariate and monotonicity in reduction operation, which is also the main computational cost of the algorithm.

Results and discussion
In this article, a new division and reduction algorithm is proposed for globally solving problem (FP). First, the original problem (FP) is converted into an equivalent optimization problem (P), in which the objective function is a single variable and the constraint functions are the difference of two increasing functions. Second, several basic operations (i.e., division, deleting, and reduction) are presented for designing a more efficient algorithm to problem (P). Finally, the numerical computational results show the feasibility and efficiency of the proposed basic operations, compared with the usual branch and bound (BB) methods (e.g., [4,5,9,10,16,17]). Additionally, as further work, we think the ideas in this article can be extended to the sum of nonlinear ratios optimization problems; for example, the numerator and denominator of each ratio in the objective function to problem (FP) are replaced with a generalized polynomial function, respectively.