A multi-parametric programming algorithm for special classes of non-convex multilevel optimization problems

A global solution strategy for multilevel optimization problems with special non-convexity formulation in the objectives of the inner level problems is presented based on branch-and-bound and multi-parametric programming approach. An algorithm to such problems is proposed by convexifying the inner level problem while the variables from upper level problems are considered as parameters. The resulting convex parametric under-estimator problem is solved using multi-parametric programming approach. A branch-and-bound procedure is employed until a pre-specified positive tolerance is satisfied. Moreover, a ε−convergence proof is given for the algorithm.


Introduction
In many real-world problems decisions have been made in a hierarchical order where individual decision makers have no direct control upon the decisions of the others, but their actions affect all other decision makers.Further, higher levels (or leaders) of the hierarchy have the power to strongly influence the performance and strategies of the decision makers at lower levels (or followers) [1].Multilevel Optimization Problems (MLOP) are mathematical programs which have a subset of their variables constrained to be an optimal solution of other programs parameterized by their remaining variables [2].It's implicitly determined by a series of optimization problems which must be solved in a predetermined sequence.
133 ⃝IJOCTA where k represents the number of levels in the hierarchy, and the decision vector partitioned into the k decision levels.Two-level problems (which are known as bilevel programming problems) and 3-level problems (usually referred to as trilevel programming problems) are common in applications [2].A trilevel optimization problem, where k = 3 in equation (1), for instance, comprises of three subproblems, one at each optimization level with the following basic definitions of sets: is called a feasible set for the third level.
• The set of solutions (5) is called the rational reaction set for the second level, for X i ⊆ R n i , i = 1, 2, 3.
One can easily see the parametric nature of the rational reaction sets, in equations ( 3) and (5), which describe the dependence of the decisions taken at the upper levels on the decisions taken at lower levels.
Because of this sequential nature of the rational reaction sets, MLOP are generally very complex and are difficult to solve.It has been shown that even the linear bilevel programming problem is strongly NP-Hard [2,3].Moreover, the complexity of the problem increases as significantly as the number of levels increase [4].
However, several algorithmic approaches have been developed that can solve convex bilevel and trilevel programming problems.For linear MLOPs vertex enumeration methods have been used in [5,6].For nonlinear MLOPs many of the methods try to transform the lower level problems by using the Karush-Kuhn-Tucker (KKT) conditions or penalty functions [7,8].However, the relaxed KKT conditions may result in nonoptimal reactions from the lower levels or in ambiguity to choose one solution among the optimal reaction set.Moreover, it seems too difficult to extend such an approach beyond two levels, because of the non-convex constraints introduced due to complementarity conditions.Recently, Tilahun et al. [9] have proposed a meta-heuristic algorithm to solve a multilevel problem of general form.The algorithm takes the variables from other levels as fixed values and solve the problem at each level (from top to bottom) using (1+1)evolutionary strategy.However, the comparison made between the previous values and the solution from each iteration, in the adaptation procedure, depends only on the objective value of each level without considering the reaction from other levels.This may create non optimal Stackelberg solutions, as the decision makers at each level optimizes its own objective function without taking the reaction from other levels into consideration.
By considering the upper level variables as parameters in the lower level problems Faísca et.al., [10,11] have proposed a Multi-parametric Programming (MPP) approach to solve MLOPs, when the lower level problems are convex.Using this approach, one can convert multilevel problems as sequential multi-parametric optimization problems.If the lower levels are convex, their parametric solutions are unique in each region, where the solution is stable.But many practical problems that are modeled using MLOP may contain non-convex terms in their lower level problems [8].In this case, even though MPP approach produces a mathematical program without equilibrium constraints depending on the upper level problems, the rational reaction set of the inner level problem (with non-convexity formulation) may be a disconnected set.So, it turns out that the upper level optimization problem may not have a solution in the rational reaction set, even if the upper level is a linear programming problem.Moreover, the global optimal solution cannot be efficiently computed and the behavior of a local solution is hard to analyze for the inner level optimization problems.
In this paper we apply the process of convexification of the lower level problems to underestimate them by convex functions (if they are nonconvex) at each iteration and use MPP approach to propose a branch-and-bound algorithm to find a global approximate solution for multilevel problems with non-convexity at their inner levels.The paper is organized as follows; in Section 2, multiparametric nonlinear programming problems are described and the respective critical regions are defined.Section 3 presents the proposed algorithm for bilevel and trilevel programming problems and convergence proof is also presented.Illustrative examples are given in Section 4. We conclude the paper by giving a conclusive remark in Section 5.

Preliminary Concepts
In MLOPs, the lower level problems can be considered as parametric optimization problems, where the decision vectors for upper levels are the parameters.To solve such kind of problems the following notions are helpful.

Multi-parametric nonlinear programming
Consider the following multi-parametric nonlinear programming problem as the inner level problem of a bilevel programming problem, where the parameter vector θ represents the vector of upper level optimization variables and x is the optimization variable of the current level: where f , g 1 , . . ., g p and h 1 , . . ., h q are twice continuously differentiable functions in x and θ.Assume also that f is a convex function and the functions g 1 , . . ., g p , and h 1 , . . ., h q define convex sets.
If x is a feasible solution to problem (6) for a given θ, we classify the constraint functions as: active constraints, where the set of constraints that satisfy the property g i (x, θ) = 0, for some i, and as inactive constraints, those that satisfy g i (x, θ) < 0 for some i.Definition 1.The active set A(x, θ) of the inequality constraints of problem (6) is the set of constraints' indices of the active constraints, that is, Definition 2. For an active set A(x, θ), we say that the linear independence constraint qualification condition holds if the set of active constraint gradients are linearly independent.
The first-order KKT optimality conditions for (6) are given for the Lagrangian: by the following conditions.
In [11] Dua et al., have proposed an algorithm to solve equation ( 8) in the entire range of the varying parameters for general convex problems.This algorithm is based on approximations of the nonlinear optimal expression, x = γ ⋆ (θ) by a set of first order approximations, as given by Corollary 1 below.
Corollary 1 (First order estimation of x(θ), λ(θ), µ(θ), near θ = θ 0 [15]).Under the considerations of Theorem 1, a first order approximation of where The space of θ where this solution (9) remains optimal is defined as the critical region (CR), and can be obtained by using feasibility and optimality conditions.Feasibility is ensured by substituting x(θ) into the inactive inequalities given in problem (6), whereas the optimality condition is given by λ(θ) ≥ 0, where the multiplier λ(θ) corresponds to the vector of active inequalities, resulting in a set of parametric constraints.Each piecewise linear approximation is confined to regions defined by feasibility and optimality conditions [11].If ǧ corresponds to the non-active constraints, and λ to the Lagrangian multipliers of the active constraints, we have: Consequently, the explicit expressions are given by a conditional piecewise linear function [11]: where C i 's are column vectors and K i 's are real matrices, whereas CR R i ⊆ R m are critical regions and note that CR R i denotes the i th critical region.

Convex relaxation of bilinear and concave terms
If the multi-parametric nonlinear optimization problem ( 6) contains bilinear and concave terms, the KKT conditions ( 7) may not produce parametric optimal solutions.To assure that KKT optimality conditions are both necessary and sufficient for obtaining the global parametric optimum of the inner level problem, the occurrence of any non-convex term must be underestimated by a convex envelope to approximate it by a convex function.
The convex envelope of bilinear terms b ij y i y j taken over the rectangle R = {(y i , y j ) : and can be obtained as follows: Theorem 2 ([16]).Let b ij , for i = 1, 2, . . ., n 2 − 1 and j = i + 1, . . ., n 2 , be a real number, then the convex envelope of a bilinear term b ij y i y j can be defined as: Theorem 2 gives us a tighter lower bound for each bilinear terms.On the other hand the over-estimator of such terms can be obtained by adding the maximum separation to the lower bound and can be stated in Corollary 2.

Corollary 2. The maximum separation between b ij y i y j and VexR
Proof.See [17] Similarly, any occurrence of univariate concave functions can be trivially underestimated by their linearizations at the lower bounds of the variable ranges [18].Thus the convex envelope of the concave function c(y) over the interval [y L , y U ] is the linear function of y obtained by: The maximum separation between the concave function and its convex under-estimator can be found by minimization problem [19] as: where, Vex c is the convex under-estimator of c(y).

Algorithm for Multi-level Optimization with Special Non-convexity Formulation at Inner Levels
Applying the formulations in subsection 2.1, it is possible to solve convex MLOP as described in [10] and [11].However, if the inner level problems are non-convex, we need to convexify them as described in subsection 2.2 and underestimate the given problem by an approximate convex problem.Using this approach we will propose an algorithm to solve nonlinear MLOP in terms of branch-and-bound procedure.In this paper, special form of non-convexity, where only bilinear and concave terms appear in the objective of the most inner level problem, will be considered and the algorithm will be tested for bilevel and trilevel problems of such a form.

Algorithm for bilevel programming problem with special non-convexity formulation at the inner problem
Consider the following non-convex bilevel problem: where g 2 (x, y) forms a convex polyhedron, h 1 (x) and h 2 (y) are linear with respect to x and y re- and X and Y are compact convex sets.Note that, the inner level is a nonconvex minimization problem.
The approach begins with rewriting the inner level of problem (10) as a MPP: As discussed in subsection (2.2) each bilinear and concave terms can be under-estimated by their respective tighter convex terms VexR[b ij y i y j ] and Vex c respectively.Thus, the under-estimator problem of ( 11) is formulated as: 12) is a linear MPP problem, which can be solved by using Linear MPP algorithm described in [14], resulting in: where, i = 1, 2, . . ., I 2 , with I 2 being the number of critical regions.Substituting the expression given in (13) into the objective function of problem (12) we obtain the parametric lower bound Ž1 i (x) for the solution of problem (11) within the corresponding critical regions, CR i .Adding the maximum separations, δ ij and M axSe, for bilinear and concave terms respectively on Ž1 i (x) gives the parametric over-estimator Ẑ1 i .Each of the lower bounds Ž1 i (x) is then compared to Ẑ1 i and the region of x where Ẑ1 i (x) − Ž1 i (x) ≤ ϵ are fathomed for a given small positive tolerance ϵ.If each of Ž1 i (x) are within ϵ of Ẑ1 i (x) in the space of x, then the expressions in (13) can be incorporated into the leader's problem of (10).Otherwise, the initial region of y is partitioned into two by bisecting the variable that has the longest side from among those that resulted in non-convexity in the problem, as in [17].After branching in y, the corresponding under-estimator subproblems in each subregion can be formulated as follows: and, where VexR1[b ij y i y j ] and VexR2[b ij y i y j ] are under-estimators of the bilinear terms over the region y L ≤ y ≤ y U new and y L new ≤ y ≤ y U respectively and Vex c 1 and Vex c 2 are convex underestimators of the concave terms over the region y L ≤ y ≤ y U new and y L new ≤ y ≤ y U respectively.Each of the problems ( 14) and ( 15) are again solved using linear MPP algorithm and one can substitute each of these parametric solutions into the objective functions of ( 14) and ( 15) respectively to obtain the parametric lower bounds Ži 21 (x) and Ži 22 (x) for the solution of problem (11).Similarly, one can get the parametric upper bound by adding the maximum separation to each of the lower bounds Ži 21 (x) and Ži 22 (x), to obtain the parametric over-estimators Ẑ21 i and Ẑ22 i .Now, by comparing each parametric upper bounds with the former upper bound for each subproblem we choose the least upper bound within the corresponding critical regions; and we shall denote it by Z u (x).Similarly, comparing the lower bounds and taking the maximum of them, we update the lower bound as well and denote it by Z l (x).Consequently, convergence test is performed by comparing the updated parametric upper bound with the updated lower bounds.If the convergence test is satisfied the branching procedure stops there.
Otherwise, the process continues in the same manner as discussed above until the lower bound and the upper bound are separated by sufficiently small positive tolerance ϵ.Following this procedure, the bilevel optimization problem will be converted to a single-level optimization problem: Due to the partitioning of the parametric region, the solution of the first-level optimization problem depends on the number of critical regions obtained in the inner problem.Based upon the above discussion an algorithm for the solution of a bilevel programming problem is presented in Table 1 Table 1

. A Parametric programming algorithm for bilevel programming problems with non-convexities in the inner problem
Step Description 1 Consider the inner problem of (10) as a MPP problem, with the upper level variables are being considered as parameters; 2 Initialize the current parametric upper bound as Zu (x) = ∞, current parametric lower bound as Zl (x) = −∞, a space of upper level optimization variable x, as a parameter space (CR) determined by the lower and upper bounds x L and x U respectively, a space of inner level optimization variable y determined by the lower and upper bounds y L and y U respectively, and tolerance, ϵ; 3 For a given region of x, CR, and the corresponding space of y, convexify the inner problem and solve using linear MPP algorithm and obtain the parametric lower bound, Ž(x) and the parametric upper bound, Ẑ(x) for the solution of the inner problem as discussed in Section 3.1; 4 Compare Zu (x) and Ẑ(x) as described in Appendix B and update the current parametric upper bound as Zu (x) = min( Zu (x), Ẑ(x)); 5 Compare Zl (x) and Ž(x) as described in Appendix B and update the current parametric lower bound as Zl (x) = max( Zl (x), Ž(x)); 6 If, Z u (x) − Z l (x) ≤ ϵ, and Z u (x) ≤ Z l (x) fathom the corresponding space of parameters and the relaxed convex optimization problems are infeasible for some rectangle, then in such cases fathom those regions and the corresponding space of parameters.7 If there is any more space of parameter x and a space of optimization variable y to be explored go to Step 8. Otherwise, go to Step 12.
Step Description 8 Branch on y and formulate a convex underestimator in each subrectangles of y; 9 Solve the convex under-estimator problems in each subrectangles using MPP approach and obtain the parametric lower and upper bounds for the solution in each subrectangles.10 Compare the lower bounds of each subrectangles within the corresponding critical regions and take, Ž(x) to be the minimum within the corresponding critical region; 11 Compare the upper bounds of each subrectangle within the corresponding critical regions and take Ẑ(x) to be the minimum within the corresponding critical region and go to Step (4).12 Substitute each of the solutions of the inner problem into the leader's problem and formulate one-level optimization problem; 13 Solve each single-level problem using suitable global optimization method; 14 Compare the leader's optimal solutions and select the best as one needs.
In order to prove the convergence of the algorithm stated in Table 1 we need to justify the following statements.Theorem 3. Consider problem (10) and assume that the optimization variables are bounded.If the reformulated inner level problem has a solution y(x) ∈ 2 R n at each iteration, then the functional sequence χ k (x) = Ẑk (x) − Žk (x), is a decreasing sequence, where Ẑk (x) is the k th least upper bound and Žk (x) is the k th greatest lower bound.Proof.To prove the theorem we need the following two arguments.
(2) Secondly, we need to show that Žk (x) is an increasing sequence.Let Žp be the lower bound for p th iteration.As a result of the algorithm, we have, Žp (x) = max{max{ Ži (x)} p−1 i=1 , Žp } and similarly at p + 1 we have, , Žp } and as the size of the rectangular domain decreases, the maximum separation between the original non-convex function and its respective convex under-estimator function decreases.This implies that Žp (x) ≤ Žp+1 (x).
From the discussion above we can realize that the parametric upper bound is found by adding the maximum separation between the underestimator and the original non-convex function over the parametric lower bound.This shows that χ k is independent of the parameter vector x.
Therefore, since Ẑk (x) is a decreasing sequence and Žk (x) is an increasing sequence, we can conclude that the difference χ k is a decreasing sequence and hence a convergent sequence.

Proof.
(1) Since CR i ⊆ X for all i and To show the backward inclusion let x ∈ X and assume that x / ∈ CR Q .Then, there exists an index i such that gi As a direct consequence of the above two theorems we have the following corollary.
k=1 is a decreasing functional sequence, then the algorithm, given by Table 1, converges.
Note that one can extend the idea discussed here above also for trilevel programming problems as discussed here below.

Algorithm for trilevel programming with non-convexity formulation in the inner problems
Consider the following non-convex trilevel optimization problem, where the second and third level objectives can have bilinear and concave terms.
where the constraint functions g 1 , g 2 and g 3 each forms a convex polyhedron, h is linear with respect to all its variables, and c(z) is a concave function.
Here again the solution procedure starts by convexifying the third-level problem and formulating a convex multi-parametric programming problem, where the variables determined by the two upper level decision makers are considered as parameters.Then the global parametric solution for the third level problem is found by the same approach as discussed in Section 3.1 for the solution of the inner level problem.
Substituting the solutions into the second-level problem, we obtain a multi-parametric programming problem, where the variables determined by the leader are the parameters.Note that since the parametric solutions for the third level problem are obtained in linear form, the complexity of the second level problem does not increase.
If the second level problem also has special non-convexity formulation, it is necessary to apply the same procedure as discussed in Section 3.1 for the inner level problem and substitute the solution into the leader's problem.Finally, the resulting single-level problem is solved using suitable global optimization methods.
On the other hand, if the second level problem is a convex multi-parametric problem, the KKT conditions are sufficient.Thus, the convexification process and the branch-and-bound procedure for the second level problem is not necessary in the above discussion.
Generally, the algorithm in Table 1 will be applied twice if the second level problem also has the same type of non-convexity as the third level.Otherwise, the relaxation step in Table 1 could be omitted for the second stage.

Algorithm for k-level programming with non-convexity formulation in the inner problems
Consider problem (1) with special non-convexity formulation at each (or some of the) optimization level(s).The k th -level optimization problem can be rewritten as a multi-parametric programming problem where, the upper levels optimization variables are considered as parameters.
The resulted problem can be solved parametrically using the algorithm described in Table 1 (from Step 2 -Step 11).The parametric solution can be incorporated into the (k − 1) thlevel.In the same manner the (k − 1) th -level optimization problem can be reformulated as multiparametric programming problem which also be solved globally using the algorithm described in

Illustrative Examples Example 1: A bilevel problem
Consider the following bilevel programming problem, where the inner problem contains a bilinear term: First of all, consider the inner level problem of (18) as a MPP problem with x a parameter vector of the optimization problem: Under-estimate problem (19) by using the formula discussed in subsection 2.2 to get: Solving problem (20) using linear MPP algorithm one can get parametric solutions within each of the corresponding critical regions Substitute the parametric solutions into the objective function of problem (20) to get the parametric lower bound of the solution with in the corresponding critical regions.These are Ž1 = 0.56x 1 − 0.06x 2 − 1.6 in CR R 1 , Ž2 = 0.56x 1 − 0.06x 2 + 1.75 in CR R 2 and Ž3 = 0.56x 1 − 0.06x 2 + 0.5 in CR R 3 .The parametric upper bound can be found by adding the maximum separation (as described in Corollary 2) δ ij = 0.125 to the lower bounds and can be described as follows: Ẑ1 = 0.56x 1 − 0.06x 2 − 1.475 in CR R 1 , Ẑ2 = 0.56x 1 − 0.06x 2 + 1.875 in CR R 2 and Ẑ3 = 0.56x 1 − 0.06x 2 + 0.625 in CR R 3 .Compare upper bounds with lower bounds within the corresponding critical regions, CR R s .But since the difference, 0.125, is greater than the pre-specified tolerance ϵ = 2.8275 × 10 −6 , branching on the optimization variable is performed.Hereby, the algorithm terminates after fifteen branching steps with parametric solutions and there is no rest region to be explored.So, one can incorporate the rational reaction set into the upper level problem and problem ( 18) becomes a single level optimization problem: Solving this problem globally one obtains an optimal solution (y 1 , y 2 , x 1 , x 2 ) = (0.1, 1, −1, −1) and an optimal value as f 2 = 0.1 and f 1 = −0.9.

Example 2: A trilevel problem
Consider the following trilevel programming problem with the occurrence of a concave term in the third level: The solution procedure starts with reformulation of the third level problem as a MPP problem by considering the optimization variables of the two upper levels as parameters and can be described as follows: min Convexify the concave term as discussed in Subsection 2.2 and reformulate problem (23) as follows: After solving problem (24) using linear MPP algorithm, one gets a parametric solution 1 Substitute the solution into the objective function of problem (24) to get the parametric lower bound Ž(x) = 3  14 x 1 + 29 28 .Consequently, add the maximum separation (MaxSe = 0.0625) over Ž(x) to obtain the parametric upper bound as Ẑ(x) = 3  14 x 1 + 29 28 + 0.0625.Compare upper bounds with lower bounds within the corresponding CR R .But the difference is greater than the pre-specified tolerance ϵ = 0.002, then branching on the optimization variable is performed.Hereby, the algorithm terminates after three branching steps with parametric solution: Compare the upper bound ẑ(x) = 0.357x 1 + 1.0595x 2 + 0.0019 with lower bound Ž(x) = 0.357x 1 + 1.0595x 2 and the difference (0.0019) is less than the pre-specified tolerance 0.002.Now we can incorporate these results into the second level problem and obtain the second level problem as: min 25) is a linear parametric programming problem, there is no need of convex relaxation and branching procedures.Hence, one can solve it using linear MPP algorithm for global optimality and get a parametric optimal solution for the second level, x 2 (x 1 ) = 2x 1 −1, and critical region CR R = {0 ≤ x 1 ≤ 1, lastly one can incorporate the result above into the leaders problem and solve it for x 1 to get the following results (x 1 , x 2 , z 1 , z 2 ) = (0.6, 0.2, 0.271, 0.047) with (f 1 , f 2 , f 3 ) = (−3.2,0.447, 0.16853)

Conclusion
The difficulty and complexity of the solution approach for MLOP is easily confirmed by looking its simplest version, what we call it linear MLOP.Especially when nonconvexity appear in the inner level problem, most of the existing algorithms fail to work.However, in this paper we have described a global strategy for the solution of MLOP with non-convexity formulation in the inner problems based on the combination of MPP approach and a branch-and-bound algorithm.The proposed algorithm is suitable for problems involving only special non-convex and linear terms in the objective functions as well as in the constraint sets.The same procedure can be applied successively to solve any multilevel problem with the same form as discussed in subsection 3.3 for k-level case.
CR rest = CR IG − CR R .For the sake of simplicity we consider the case when only two parameters, θ 1 and θ 2 , with CR R defined by three inequalities, {g 1 ≤ 0, g 2 ≤ 0, g 3 ≤ 0} where g 1 , g 2 , g 3 are linear in θ.The procedure consists of considering one by one the inequalities which define CR R .For example, consider inequality g 1 ≤ 0, the rest of the region can be addressed by reversing the sign of inequality g 1 ≤ 0 and removing redundant constraints in CR IG , which is CR rest

Appendix B: Comparison of parametric solutions
A method has been proposed in [14] for the comparison of two parametric solutions, Z(θ) 1  and Z(θ) 2 which are valid in the critical regions CR R 1 and CR R 2 respectively.The comparison process needs two steps.The first step is to define a region In the second step, check whether CR int is empty or not.If CR int is empty, then there is no comparison to be performed, otherwise a new constraint Z(θ) 1 ≤ Z(θ) 2 is formulated and a constraint redundancy check is made for the new constraint in CR int .This constraint redundancy test results in three cases which are analyzed as follows: Case 1: If the new constraint is redundant (see [14]), then Z(θ) 1 ≤ Z(θ) 2 , ∀θ ∈ CR int Case 2: If the new constraint is infeasible (see [13]), then Z(θ) 1 ≥ Z(θ) 2 , ∀θ ∈ CR int Case 3: If the new constraint is non-redundant, then

Table 2 .
θ 2 ≤ θ U 2 }.Thus by considering the rest of the inequalities, the total of the rest region is given by, CR rest = {CR rest Definition of the rest regions