Solving bilevel problems with products of upper- and lower-level variables
Sina Hajikazemi Florian Steinke

Photoelectron – residual-ion entanglement in streaked shake-up ionization of helium
    Hongyu Shi and Uwe Thumm
    September 9, 2024
===================================================================================



§ ABSTRACT 


Bilevel programming problems frequently arise in real-world applications across various fields, including transportation, economics, energy markets and healthcare.
These problems have been proven to be NP-hard even in the simplest form with linear upper and lower-level problems.

This paper addresses a specific type of bilevel programming problem where the upper-level is linear, and the lower level includes bilinear terms involving product of variables from both levels.
We propose a new iterative algorithm that addresses this specific class of bilevel problems by penalizing the duality gap and linearizing the bilinear terms.


The effectiveness of the algorithm is argued and demonstrated through a numerical example.






§ INTRODUCTION


Bilevel programming problems frequently arise in real-world applications across various fields, including transportation, economics, energy markets and healthcare.
This problem is NP-hard even in the simplest form with linear upper and lower-level problems<cit.>.

A common approach to solve the LP-LP bilevel programming problem to the global optimum is to replace the lower-level optimization problem with its KKT conditions, which turns it into a single-level optimization problem.
However, the complementary slackness constraints involve nonlinear terms.
These constraints can be linearized using either binary variables and a big-M value <cit.> or SOS1 constraints <cit.>.
Both approaches are typically very slow to solve and are not suitable for large problems.
In addition, finding an appropriate big-M value is very challenging;
<cit.> shows that finding an appropriate big-M value is as hard as solving the original problem.
Another method to solve this problem is the Penalty Alternating Direction Method (PADM) <cit.>. Unlike the previous methods
it only converges to a partial optimum, but it is computationally feasible for large problems and is shown to yield solutions close to the global optimum in a large fraction of sample problems. 


The problem studied in this research, is linear in the upper level and nonlinear in the lower level, with bilinear terms of the upper and lower-level variables in the constraints.
It is easy to prove that any LP-LP bilevel problem can be reduced to this problem by defining and fixing some lower level variables. Thus, the studied problem is harder than an LP-LP bilevel problem.
To the best of our knowledge, there is no prior research that addresses this specific class of problems. 

To solve the proposed, computationally hard problem,
we design a scalable, iterative algorithm called Penalty Adaptive Linearization Method (PALM).
This iterative algorithm penalizes the optimality gap of the lower-level problem 
similar to PADM <cit.>, balancing the optimality of the upper- and lower-level problems with increasing weight on the lower-level optimality gap during the course of the algorithm execution.
During this process, it also linearizes the bilinear terms in the constraints and iteratively readjusts the linearization.



The rest of the paper is organized as follows.
Section <ref> formulates the problem.
Section <ref> introduces the proposed algorithm.
Section <ref> demonstrates the algorithm with a small numerical example.
Finally, section <ref> discusses the results and outlines possible future work.



§ PROBLEM FORMULATION

The studied problem consists of two levels.
The upper-level problem is given in (<ref>), where A ∈ℝ^p× mn, B ∈ℝ^p× n, c ∈ℝ^mn, d ∈ℝ^n and a ∈ℝ^p are the coefficients.
The variables X ∈ℝ^m × n and y ∈ℝ^n correspond to the upper and lower-level problems, respectively.
Additionally, function vec is used to vectorize the elements of a matrix,
and S(X) denotes the point to set mapping from X to the optimal solutions of the lower-level problem.


    
    min_X,y      c^T vec(X) + d^Ty 
        A  vec(X) + B y ≥ a 
        y ∈ S(X)





The lower-level problem is given in (<ref>), where C ∈ℝ^m× n, b ∈ℝ^m and e ∈ℝ^n are the coefficients.
Constraint (<ref>) in the upper-level problem ensures that y attains the optimal value of the lower-level problem.

    
    min_y     e^T y 
    
        s.t.      (C+X)  y ≥ b







§ PROPOSED SOLUTION METHOD

To obtain a single-level formulation of the bilevel problem, the lower-level problem can be replaced by a set of optimality conditions.
For a linear programming problem, these can be primal feasibility, dual feasibility, and strong duality constraints. The strong duality constraint forces the objective function of the primal minimization problem to be less than or equal to the objective function of the dual problem.

Since the lower-level problem (<ref>) is linear with respect to the primal variables y, its dual is given by (<ref>), where λ∈ℝ^m denotes the dual variables.


    max_λ     b^T λ
    
    s.t.      (C+X)^T λ = e 
       λ≥ 0


Replacing (<ref>) with the above mentioned optimality conditions yields the following single-level reformulation of (<ref>): 


    min_X,y,λ      c^T vec(X) + d^Ty       
    
    s.t.       (<ref>)       upper-level constraints
        (<ref>)         primal constraints
        (<ref>)-(<ref>)       dual constraints
         e^Ty ≤ b^T λ      strong duality constraint


Constraint (<ref>) is the strong duality constraint.
Including this constraint along with the primal and dual feasibility constraints strictly enforces the optimality of the lower-level problem.
Problem (<ref>) is a non-convex quadratically constrained problem which are NP-hard to solve globally in general <cit.>. 

Even finding a feasible point is difficult due to the strictness of the strong duality constraint (<ref>). 



Hence, we proceed iteratively and  move the strong duality constraint into the objective function as a penalty term, gradually increasing its weight to progressively enforce the optimality of the lower-level problem, similar to the PADM method.
Shifting the strong duality constraint to the objective function as a penalty term with weight μ results in the following formulation:

	
    min_X,y,λ      c^T vec(X) + d^Ty + μ[ e^Ty - b^T λ] 
    
    	s.t.       (<ref>)       upper-level constraints
        (<ref>)         primal constraints
        (<ref>)-(<ref>)       dual constraints


Constraints (<ref>)  and (<ref>) both contain bilinear terms which makes problem (<ref>) non-convex and difficult to solve.
Our approach is to linearize the bilinear terms in (<ref>) and solve the problem iteratively.


To linearize the Xy term in (<ref>), we approximate it around a point (X̅,y̅).
Let dX ∈ℝ^m× n and dy ∈ℝ^n be defined as dX = X - X̅ and dy = y -y̅.
Then the approximation is as follows:

    X y     = X̅ y̅ + X̅  dy + dX  y̅ + dX   dy 
       ≈X̅ y +  dX  y̅

Assuming that dX and dy are small, the term dX   dy is ignored in the second step. 

This approximation helps to linearize the bilinear terms in (<ref>).
The same approximation can be applied to linearize the X^Tλ terms in (<ref>).

Replacing the bilinear terms with the linear approximation changes the primal form of the lower-level problem (<ref>) to (<ref>) and the dual form of the lower-level problem (<ref>) to (<ref>).


	
    min_y      e^T y 
        (C+X̅) y + dX  y̅≥ b




	
    max_λ      b^T λ
        (C+X̅)^T λ   + dX^T  λ̅= e 
       λ≥ 0







Substituting the primal and dual feasibility constraints (<ref>) and (<ref>) in problem (<ref>) with (<ref>) and (<ref>),
yields the following formulation:


    min_dX,y,λ      c^T vec(X̅+dX) + d^Ty + μ[ e^Ty - b^T λ]  
    
    s.t.       (<ref>)  
        (<ref>)   
        (<ref>)-(<ref>)


The problem given in (<ref>) is linear. Let S^L(X̅,y̅,λ̅) denote the set of optimal values of (<ref>).



The proposed algorithm is presented in <ref> and consists of two loops.
The first (outer) loop (starting at line <ref>) evaluates the optimality of the lower-level problem. 
At the end of each iteration, the penalty coefficient μ is doubled, thereby increasing the priority of the optimality of the lower-level problem over the upper-level objective function.
This loop terminates when the optimality of the lower-level problem is satisfied, that means a feasible solution to the upper-level problem is found.
If it does not terminate after a certain number of iterations, it means that either the problem is infeasible or the algorithm is not able to find a feasible solution.
The second loop (inner) (starting at line <ref>), solves the linear problem <ref>, updates the values of X̅, y̅ and λ̅ and reapeats until these values converge.

Fixing the variable X to X̅ in problem (<ref>) results in a linear problem. Moreover, the constraints involving y and λ in are separable in (<ref>), so the problem can be decomposed into two subproblems (<ref>) and (<ref>), which are solved independently.

    min_y       d^Ty + μ  e^Ty   
    
    	s.t.       A  vec(X̅) + B y ≥ a  
        (C+X̅)  y ≥ b


    max_λ       b^T λ
    
    	s.t.        (C+X̅)^T λ = e 
       λ≥ 0


So for a given X̅, the solutions to the subproblems (<ref>) and (<ref>) determine the values of y̅ and λ̅, respectively.
Let S^P(X̅) and S^D(X̅) be the optimal solution sets of the subproblems (<ref>) and (<ref>), respectively.
If the solutions of (<ref>) and (<ref>) are not unique for all X̅, then y̅ and λ̅ may jump between different alternative solutions within the sets S^P(X̅) and S^D(X̅).
This can cause the algorithm to fail to converge in the inner loop. To mitigate this, we select the optimal solution that is closest to the previous one. 
To find it we first compute the optimal objective value by solving the LP. We then resolve it with the objective value fixed to the optimal value and minimizing the distance to the prior solution. Depending on whether the 1-norm or 2-norm is used, this results in a linear or convex quadratic programming problem, respectively. The values of y̅ and λ̅ are then updated according to lines <ref> and <ref>.
Similarly, when choosing among the alternative optimal solutions in S^L(X̅,y̅,λ̅), we choose the solution that results in the smallest change in X̅.



The algorithm needs an initial value X_0 for the upper-level variable, for which a feasible primal and dual solution exists to start with.
Let Ω be the set of feasible solutions of problem (<ref>).
Then the initial values for y and λ are computed as in line <ref>.
Note that this problem is linear since X is fixed to X̅.



Since the bilinear terms are omitted in the linearization of both the primal and dual constraints, the solution to the linear problem (<ref>) may not be feasible for the original problem (<ref>).
Nonetheless, the step size dX^* tends to increase as the penalty parameter μ is raised, allowing indirect control over dX^* through adjustments to μ.
If the solution to the linear problem is not feasible for the original problem, then slowing the rate of increase of the penalty parameter  μ can help reduce the step size dX^* and improve the likelihood of achieving a feasible solution.



§ A MINIMAL EXAMPLE

In this section a minimal numerical example is presented to show how the algorithm works.
The notation used here aligns with that in Section <ref> except the upper-level variables which is denoted by x instead of a matrix X.
The upper-level problem is as follows: 


    min_x,y_1,y_2      |x|  
    
    s.t.       y_2 ≤ 1.5
        y_1, y_2 ∈ S(x)



The absolute value function can easily be replaced with an auxiliary variable and two constraints.
The lower-level problem is as follows:


    min_y_1,y_2      y_1 + y_2  
    
    s.t.       0.5y_1 + y_2 + x y_1 ≥ 3
        y_1 + 0.5 y_2 - x y_1 ≥ 3 
        y_1, y_2 ≥ 0




The problem is solved to the global optimum using a non-convex QCQP solver.
The optimal solution is y^*=(2.5, 1.5) and x^*=0.1.

Then the Penalty Adaptive Linearization Method is applied to the problem.
The initial value of the upper-level variable is set to zero and the optimal solution of the lower-level problem is y^*=(2.0, 2.0) for this inital value.
The algorithm successfuly converges to the optimal solution of the problem, which is y^*=(2.5, 1.5) and x^*=0.1.
Figure <ref> shows the convergence of the duality gap of the lower-level problem to zero which means that the algorithm converges to a feasible solution of the main problem.
As expected, the optimal value y̅_̅2̅ remain constant at 1.5 throughout the iterations.










§ CONLCUSION


In this paper we designed a new algorithm for solving bilevel programming problems with product of upper and lower-level variables in the lower-level problem.

We showed that the algorithm was able to iteratively approximate the bilinear terms in the sample problem and solve it to the global optimum.


The most crucial aspect missing in this algorithm is to prove that the inner loop is guaranteed to converge.
It’s possible that convergence only occurs under specific conditions.
In such cases, identifying non-convergent scenarios can help clarify the algorithm’s limitations, highlight areas for improvement, and specify the conditions under which convergence to a partial minimum is guaranteed.

From a practical perspective, exploring the algorithm’s performance in solving a wide range of real-world problems is of great interest.
However, unlike LP-LP bilevel programming problems, using KKT conditions or strong duality conditions along with the big-M method does not linearize the problem, making it challenging to determine the global optimal solution.




§ ACKNOWLEDGEMENTS

This research was funded by the German Federal Ministry of Education and Research Infrastructure in project RODES (grant number 05M22RDA). 






