This paper addresses the bilevel programming problems (BPPs) with the quadratic objective functions at the upper and the lower levels. The new solution method for such BPPs is developed. The main feature of the approach proposed is employment of the original A.S. Strekalovsky's Global Search Theory. Numerical testing of the method on specially constructed instances demonstrated the efficiency of the approach.
Introduction variables of both levels which are in cooperation: where > 0 is a penalty parameter and D := {(x; y; v) | Ax ≤ b; D1y + d1 + xQ + vB1 = 0; v ≥ 0; A1x + B1y ≤ b1}. We can prove the standard results about the connection between problems (DCC) and (DC( )) [Dempe, 2002, Strekalovsky et al., 2010]. In particular, if the triplet (x( ); y( ); v( )) is a solution of problem (DC( )) and ⟨v( ); b1 − A1x( ) − B1y( )⟩ = 0, then (x( ); y( ); v( )) turns out to be a solution of problem (DCC). Moreover, it can be shown that there exists a finite value of with ⟨v( ); b1 − A1x( ) − B1y( )⟩ = 0.
Note that for a fixed problem (DC( )) has a bilinear structure and belongs to the class of d.c. minimization problems [Strekalovsky, 2003, Strekalovsky, 2014] with a convex feasible set. It is easy to see that the objective function of (DC( )) can be represented as a difference of two convex functions [Strekalovsky, 2003, Strekalovsky, 2014].
We will employ, for example, the following d.c. representation based on the known property of a scalar product:
Φ(x; y; v) = g(x; y; v) − h(x; y; v);
1 1
where g(x; y; v) = 2 ⟨x; Cx⟩ + ⟨c; x⟩ + 2 ⟨y; Dy⟩ + ⟨d; y⟩ + ⟨b1; v⟩ + 4 ∥A1x − v∥2 + 4 ∥B1y − v∥2,
h(x; y; v) = 4 ∥A1x + v∥2 + 4 ∥B1y + v∥2. Note that the so-called basic nonconvexity in Problem (DC( )) is
provided by the function h
So we can apply the Global Search Theory in d.c. minimization problems [Strekalovsky, 2003, Strekalovsky, 2014] to seek a global solution to Problem (DC( )) with a fixed . 3
Global Optimality Conditions and Solution Algorithm The Global Search Procedure is based on the Global Optimality Conditions (GOCs) developed by A.S. Strekalovsky [Strekalovsky, 2003, Strekalovsky, 2014]. In particular, the necessary Global Optimality Conditions have the following form in terms of Problem (DC( )).
Theorem 1. [Strekalovsky, 2003, Strekalovsky, 2014]. If the feasible point (x∗; y∗; v∗) is a (global) solution to Problem (DC( )), then ∀(z; u; w; ) ∈ Rm+n+q+1 : h(z; u; w) = − ; := Φ(x∗; y∗; v∗); g(x; y; v) ≤ ≤ sup (g; D);
x;y;v (QBP) (1) (2) the following inequality holds g(x; y; v) − ≥ ⟨∇h(z; u; w); (x; y; v) − (z; u; w)⟩ ∀(x; y; v) ∈ D: (3)
The conditions (2)-(3) possess the so-called algorithmic (constructive) property: if the GOCs are violated, we can construct a feasible point that will be better than the (critical) point in question (the value of the objective function in the new point will be less) [Strekalovsky, 2003, Strekalovsky, 2014, Orlov, 2017].
The constructive property is a foundation of global search algorithms for nonconvex problems [Strekalovsky, 2003, Strekalovsky, 2014]. Taking into account the d.c. representation (1), on the basis of Theorem 1, the Global Search Algorithm (GSA) for quadratic bilevel problems can be formulated in the following way.
Let there be given a starting point (x0; y0; v0) ∈ D, numerical sequences { k}; { k}. ( k; k > 0; k = 0; 1; 2; :::; k ↓ 0; k ↓ 0; (k → ∞)), a set Dir = {(z¯1; u¯1; w¯1); :::; (z¯N ; u¯N ; w¯N ) ∈ IRm+n+q |(z¯i; u¯i; w¯i) ̸= 0; i = 1; :::; N }, the numbers − =△ inf(g; D) and + =△ sup(g; D), and the algorithm’s parameters M and .
Step 0. Set k := 0; (x¯k; y¯k; v¯k) := (x0; y0; v0); i := 1; := −; △ = ( + − −)=M .
Step 1. Proceeding from the point (x¯k; y¯k; v¯k) by a local search method, build a k-critical point (xk; yk; vk) ∈ D to Problem (DC( )). Set k := Φ(xk; yk; vk):
Step 2. Using (z¯i; u¯i; w¯i) ∈ Dir, construct a point (zi; ui; wi) of the approximation Ak = {(z1; u1; w1); :::; (zN ; uN ; wN ) | h(zi; ui; wi) = − k; i = 1; :::; N } of the level surface U ( k) = {(x; y; v) | h(x; y; v) = − k} of the convex function h(x; y; z), such that h(zi; ui; wi) = − k.
Step 3. If g(zi; ui; wi) > + , then i := i + 1 and return to Step 2.
Step 4. Find a k-solution (x¯i; y¯i; v¯i) of the following linearized problem: g(x; y; v) − ⟨∇h(zi; ui; wi); (x; y; v)⟩ ↓ xm;yi;nv; (x; y; v) ∈ D: (PLi)
Step 5. Starting at the point (x¯i; y¯i; v¯i), build a k-critical point (xˆi; yˆi; vˆi) ∈ D to Problem (DC( )) by means of the local search method.
Step 6. If Φ(xˆi; yˆi; vˆi) ≥ Φ(xk; yk; vk); i < N; then set i := i + 1 and return to Step 2.
Step 7. If Φ(xˆi; yˆi; vˆi) ≥ Φ(xk; yk; vk); i = N and < +, then set := + △ , i := 1 and go to Step 2.
Step 8. If Φ(xˆi; yˆi; vˆi) < Φ(xk; yk; vk); then set := −, (x¯k+1; y¯k+1; v¯k+1) := (xˆi; yˆi; vˆi), k := k + 1, i := 1 and return to Step 1.
Step 9. If Φ(xˆi; yˆi; vˆi) ≥ Φ(xk; yk; vk); i = N and = +, then stop. (xk; yk; vk) is the obtained solution of the problem.
It can be readily seen that this algorithm is not an algorithm in the conventional sense, because some of its steps are not specified. For example, we do not know how to construct a starting point and the set Dir, how to implement a local search, how to solve the problem (PLi) etc. These issues will be considered below. 4 Implementation of the Global Search Algorithm First, to construct a feasible starting point, we used the projection of the chosen infeasible point (x0; y0; v0) onto the feasible set D by solving the following quadratic programming problem: ∥(x; y; v) − (x0; y0; v0)∥2 ↓ xm;yi;nv; (x; y; v) ∈ D: (PR(x0; y0; v0)) The solution to Problem (PR(x0; y0; v0)) was taken as a starting point (x0; y0; v0) ∈ D. In this work (x0; y0; v0) = (0; 0; 0). The value of the penalty parameter was chosen experimentally: = 10.
The local search (see Steps 1 and 5) can be based on the consecutive solution of the following convex quadratic (QP) and linear programming (LP) problems derived from Problem (DC( )): 1 1 2 ⟨x; Cx⟩ + ⟨c; x⟩ + 2 ⟨y; Dy⟩ + ⟨d; y⟩ − ⟨vsA1; x⟩ − ⟨vsB1; y⟩ ↓ mx;iyn;
Ax ≤ b; A1x + B1y ≤ b1; D1y + d1 + xQ + vsB1 = 0:
⟨b1 − A1xs − B1ys; v⟩ ↓ mvin; D1ys+1 + d1 + xsQ + vB1 = 0; v ≥ 0;
}
(QP(vs)) (LP(xs; ys)) where (xs; ys; vs) ∈ D is a feasible point in Problem (DC( )). Such local search methods show their efficiency in optimization problems with bilinear structure [Orlov et al., 2016, Orlov & Strekalovsky, 2016, Strekalovsky & Orlov, 2007, Strekalovsky et al., 2010, Strekalovsky, 2014]. Note that the accuracy for auxiliary problems is s = 10−7. The accuracy of the local search is k = 10−5.
The key element of the above GSA consists in constructing an approximation of the level surface of the convex function h(·), which generates the basic nonconvexity in the problem under consideration. For Problem (DC( )) the approximation Ak = A( k) has been constructed with the help of a special set of directions Dir = {((x; y) + el; v + ej ); l = 1; :::; m + n; j = 1; :::; q}; where el ∈ IRm+n; ej ∈ IRq are the Euclidean basis vectors of the corresponding dimension, (x; y; v) is a current critical point.
Dir is the standard direction set for problems with a bilinear structure according to our previous experience [Orlov & Strekalovsky, 2005, Orlov, 2008, Orlov et al., 2016, Strekalovsky et al., 2010]. Unfortunately, we cannot theoretically guarantee the global optimality of the point generated by the GSA with the set Dir. But numerically we obtain global solutions in most cases. In addition, we apply a special technique for reducing the approximation, because the number of points in the original set is rather large (especially when the dimension of the problem grows).
Consider an arbitrary set of directions with q(m + n) points:
Dir0 = {(z¯i; w¯j ) | (z¯i; w¯j ) ̸= 0; i = 1; :::; m + n; j = 1; :::; q}: In the matrices A1 and B1, find rows and columns with the maximum sum of their elements. Denote them as iA; iB and jA; jB, respectively. Then the new direction set will have the following form:
Cut(Dir0) = {(z¯iA ; w¯j ); (z¯iB ; w¯j ); j = 1; :::; q; (z¯i; w¯jA ); (z¯i; w¯jB ); i = 1; :::; m + n}: The number of points in the approximation based on the set Cut(Dir0) is equal to 2(q + m + n). Therefore, as the dimension of the problem grows, it increases slower than the number of points in the approximation based on the set Dir0. Note that here we use the matrices A1 and B1, because they are included in the definition of the function h(·) which generates the basic nonconvexity of the problem in question.
We employ the following technique to construct approximation points Ak = A( k) on the basis of the direction
set: the triples (zi; ui; wi) are found in the form (zi; ui; wi) = i(z¯i; u¯i; w¯i); i = 1; :::; N , where i ∈ IR are
computed using the condition h( i(z¯i; u¯i; w¯i)) = − k. In that case the search of i can be performed analytically
Further note that the selection of the algorithm parameters M and can be carried out on the basis of our previous experience in solving problems with a bilinear structure [Orlov & Strekalovsky, 2005, Orlov, 2008, Orlov et al., 2016, Strekalovsky et al., 2010]. The parameter is responsible for the accuracy of the inequality (2) from the GOCs (in order to diminish the computer rounding errors) [Orlov & Strekalovsky, 2005, Orlov, 2008, Orlov et al., 2016, Strekalovsky et al., 2010] (Step 3). Different values of the parameter M are responsible for splitting the interval [ −; +] into a suitable number of parts to realize a passive one-dimensional search along . Here we use the following sets: 1) M = 2, = 0:0; 2) M = 5, = 0:02; 3) M = 33, = 0:1. If it is required that an approximation to the global solution to Problem (DC( )) be found rapidly, we may use option 1). With increase of M and (options 2) and 3)), the algorithm gains in precision but loses in performance rate.
To compute the segment [ −; +] for one-dimensional search according to the GOCs, we need to solve two problems: on the minimum and maximum of a convex quadratic function g(·). The minimum problem can be solved by any quadratic programming method and an appropriate software subroutine. By the way, the same is true for the linearized problem (PLi) at Step 4 ( k = 10−5). To tackle the maximum problem, we can employ a known global search strategy for convex maximization problems [Strekalovsky, 2003]. But in this case the computational process does not require an exact knowledge of these bounds. It is sufficient to have comparatively rough estimates [Strekalovsky, 2003]. Therefore, here we use: − := 0:0; + := (m + n + l) ∗ :
Finally, Steps 6-9 represent verification of the main inequality (3) from the GOCs, the stopping criteria, and looping. 5
Test Problem Generation Method and Numerical Computations
One of the most important issues in the testing of new numerical methods is the selection or
construction of test cases. In the present work, we use the method for generation of bilevel test cases proposed in
[Calamai & Vicente, 1994]. The idea of such generation is based on constructing bilevel problems of an arbitrary
dimension with the help of the so-called kernel problems, which are one-dimensional bilevel problems with known
local and global solutions
In generation of quadratic bilevel test problems of the type (QBP), we used the kernel problems of the form F (x; y) = ∑i=1( 12 xi2 − xi + 12 yi2) ↓ mx;iyn; r
x ∈ X = {x ∈ IRm | xi ≥ 0; i = 1; :::; r}; y ∈ Y∗(x) = Arg myin{∑i=1( 21 yi2 − xiyi) | y ∈ Y (x)}; Y (x) = {y ∈ IRn | xi − yi ≤ 1; xi + yi ≤ i; − xi − yi ≤ −1}; where i ∈ {1; 1:5; 2; 3}, i = 1; :::; r.
Let cl1; cl2; cl3; cl4 be the number of kernel problems of each class included into the “big” problem. Then the dimension of the “big” problem will be m = n = p = cl1 + cl2 + cl3 + cl4; q = 3 ∗ (cl1 + cl2 + cl3 + cl4): Note that we can compute the number local and global solutions to the “big” problem.
Proposition 1 [Calamai & Vicente, 1994]. Problem (5) has 2cl3 global solutions, and it has additional local solutions only when cl2 + cl4 > 0. In this case the number of additional local minima to problem (5) equals 2cl2+cl3+cl4 − 2cl3.
Let us now describe the numerical testing of the global search method on series of problems generated by the method described above. The software that implements the method developed was coded in MATLAB 7.11.0.584 R2010b. To run the software, we used the computer with Intel Core i5-2400 processor (3.1 GHz) and 4Gb RAM. In total, 333 problems of dimension from 2 up to 200 were generated and solved.
The most interesting and typical results are presented in Table 1 with the following denotations: N N = m = n = p is the dimension of the problem; M= is the number of the set of parameters by which we can find a solution; GIt is the number of iterations of the global search method; Loc stands for the number of start-ups of the local search procedure performed to find the approximate global solution to the problem; LP , QP are the numbers of the LP and QP problems solved, respectively; T is the operating time of the program (in seconds); Glob=Loc are the numbers of global solutions in the “big” problem and local solutions which are not global, respectively.
First of all, note that all generated problems were solved with the prescribed accuracy " = 10−3. Further, pay attention to a huge number of global and local solutions in some cases. It is interesting that the hardness of a problem does not depend on these values directly. Here we can see easy problems which are solved at the local search stage. And we can also see hard problems which take more than one hour to find a global solution. We propose a hypothesis that the hardness of problems depends on the set of classes of kernel problems in the “big” problem. At the end of the paper we present an overall statistics about complexity of classes combinations in the “big” problem.
Let us introduce a special complexity coefficient (CC): CC = (Loc=N N )=Cnt, where Loc is the number of the local search procedures for a given classes combination; N N = m = n = p is the dimension of the problem with the given classes combination; Cnt is the total number of problems solved with the given classes combination.
The diagram about values of CC in different cases is presented in Figure 1, where the denotation 0x0x, for example, means that in the “big” problem we use kernel problems of classes 2 and 4; at the same time, the denotation 00x0 means that in the “big” problem we use the kernel problems of class 3 only etc. (4) (5)
So, we can conclude that the hardest problems either consist of all four classes of kernel problems or
comprise 2, 3, and 4 classes. This conclusion will be used in our future works when we will address more
complicated test problems with a special random transformation
Acknowledgements This work has been supported by the Russian Science Foundation (Project no. 15-11-20015).