Using neural networks to solve linear bilevel problems with unknown lower level

Bilevel problems are used to model the interaction between two decision makers in which the lower-level problem, the so-called follower’s problem, appears as a constraint in the upper-level problem of the so-called leader. One issue in many practical situations is that the follower’s problem is not explicitly known by the leader. For such bilevel problems with unknown lower-level model we propose the use of neural networks to learn the follower’s optimal response for given decisions of the leader based on available historical data of pairs of leader and follower decisions. Integrating the resulting neural network in a single-level reformulation of the bilevel problem leads to a challenging model with a black-box constraint. We exploit Lipschitz optimization techniques from the literature to solve this reformulation and illustrate the applicability of the proposed method with some preliminary case studies using academic and linear bilevel instances.


Introduction
Bilevel optimization has been a highly active field of research in the last decades and has gained increasing attention over the last years. The main reason is that this class of optimization problems can serve as a powerful modeling tool in situations in which one has to model hierarchical decision making; see, e.g., the recent surveys by Beck et al. [18] and Kleinert et al. [2] as well as the annotated bibliography by Dempe [10] to get an overview of the many applications of bilevel optimization. However, this ability also renders bilevel optimization problems very hard to solve 1 3 both in theory and practice. For instance, even linear bilevel problems are strongly NP-hard [15].
As discussed in the survey by Beck et al. [2], bilevel optimization problems can also be subject to uncertainty. Moreover, the sources of potential uncertainties are even richer compared to usual, i.e., single-level, optimization. The reason is that not only the problem's data can be uncertain but also the (observation of the) decisions of the two players can be noisy or uncertain. For instance, it might be the case that the leader is not sure about whether the follower can solve the respective lower-level problem to global optimality and might want to hedge against possibly occurring near-optimal solutions; see, e.g., Besançon et al. [3,4]. In this short paper, we go even one step further and assume that the leader has no knowledge about an explicit formulation of the lower-level problem. Hence, the leader needs to solve a bilevel optimization problem and the lower level is unknown. What might sound impossible at a first glance can be done, at least approximately, if the bilevel game is repeatedly played so that past pairs of leader decisions and respective responses of the follower can be observed and collected. The obtained set of pairs of decisions can then be used as training data to train a neural network that learns the best-response function of the follower.
The downside of this approach is that the input-output mapping of the neural network can again not be stated using a closed-form expression. This leads to the field of optimization under black-box constraints [23]. Fortunately, there is some recent research on deep neural networks with specific activation functions that shows that the input-output mapping of the neural network is Lipschitz continuous and that Lipschitz constants can actually be computed using specifically chosen semidefinite programs [11,22]. This paves the way for applying the Lipschitz optimization method proposed in Schmidt et al. [25], see Schmidt et al. [26] as well, in order to solve the single-level reformulation of the bilevel problem with unknown lower level in which the optimal response of the follower is replaced by the mapping that describes the input-output behavior of the trained neural network. In the Lipschitzbased approach, the only nonlinearity, which is the one given by the input-output mapping of the neural network, is outer-approximated using the Lipschitz constant of the mapping. This outer approximation is tightened from iteration to iteration, leading to a union of polyhedra that form a relaxation of the graph of the nonlinearity, which enables the application of powerful state-of-the-art mixed-integer solvers.
The main contribution of this short note is the combination of the recent results about the Lipschitz continuity of special neural networks with recent publications on Lipschitz optimization. By carrying this out in a careful way and by slightly adapting the method proposed in Schmidt et al. [25], we obtain a convergent algorithm for this highly challenging situation. We further illustrate the applicability of our approach in a case study based on academic bilevel instances from the literature.
Although there have been some applications of bilevel optimization problems for machine learning (see, e.g., Table 1 in Khanduri et al. [17]), this is, to the best of our knowledge, the first paper that uses neural networks to solve general bilevel optimization problems. The only other paper we are aware of following this idea is the one by Vlah et al. [27], who apply deep convolutional neural networks to classic bilevel bidding problems in power markets. However, they focus on the specific application and the specific type of network and its training. In contrast, we focus on the general mathematical idea of replacing unknown lower levels with neutral networks. Furthermore, the work by Borrero et al. [7] presents a sequential learning method for linear bilevel problems under incomplete information. Provided that the strategic players interact with each other multiple times, the authors develop feedback mechanisms to update the missing information for the lower-level objective. Nevertheless, this way to address uncertainty and learning clearly differs from our neural-network based approach. The main idea of this short paper is to give a proof of concept for the proposed method and we, thus, make some simplifications such as that we only consider linear bilevel problems with a scalar upper-level variable and without coupling constraints. We discuss these assumptions in more detail in Sect. 2 and how the setting can be generalized in Sect. 6. Let us finally comment on that our approach is related to other methods for solving bilevel optimization problems that rely on using optimal-value or best-response functions (see, e.g., Lozano and Smith [20] and the references therein). However, our approach is different since we do not work with these functions themselves but with surrogates that we obtain from training a neural network.
The remainder of the paper is structured as follows. In Sect. 2, we introduce the class of bilevel problems that we consider and pose the main assumptions that are required in what follows. Afterward, in Sect. 3, we then review the recent literature about computing Lipschitz constants for deep neural networks, which is a prerequisite for the overall solution approach discussed in Sect. 4. Section 5 contains the case studies that illustrate the applicability of our approach for small instances from the literature. Finally, we conclude in Sect. 6 and discuss some topics for future research.

Problem statement
We consider linear bilevel problems of the general form where Ψ(x) is the set of optimal solutions of the x-parameterized problem with c ∈ ℝ n x , d, f ∈ ℝ n y , A ∈ ℝ m×n x , C ∈ ℝ ×n x , D ∈ ℝ ×n y , a ∈ ℝ m , and b ∈ ℝ . Problem (1) is called the upper-level or leader's problem. The decision variables of (1a) min x∈X,y the leader are x ∈ ℝ n x . Problem (2) is called the lower-level or follower's problem and has the decision variables y ∈ ℝ n y . Note that we consider the optimistic bilevel problem. Hence, the leader also optimizes over y if the lower-level problem is not uniquely solvable. The set Ω ∶= {(x, y) ∈ X × Y ∶ (1b), (1c), (2b)} is called the shared constraint set and the set F ∶= {(x, y) ∈ Ω ∶ y ∈ Ψ(x)} is called the bilevel feasible set. With F x , we denote its projection onto the x variables. For what follows, we make the ensuing assumptions.
(ii) For all upper-level decisions x for which Ψ(x) is non-empty, we have |Ψ(x)| = 1 . Hence, if the lower level is feasible and bounded, then it has a unique optimal solution y. Consequently, we can write y = Ψ(x) instead of y ∈ Ψ(x). (iii) The solution-set mapping Ψ(⋅) is Lipschitz continuous.
According to Assumption 1, Ψ(x) is a singleton, leading to a one-to-one correspondence between follower's optimal responses and the upper-level decisions. Moreover, the mapping Ψ(⋅) is polyhedral, i.e., its graph is a finite union of polyhedra; see Theorem 3.1 in Dempe [9]. Furthermore, the bilevel feasible set is connected because we have no coupling constraints and, thus, Ψ i (x) is a Lipschitz continuous and piecewise linear function in x for all i ∈ [n y ] ∶= {1, … , n y }.
Lower-level uniqueness is particularly important when it comes to training a neural network to learn the optimal response of the follower for a given upper-level decision. It ensures that, during supervised learning, there is only a single output y to be learned for a given input x. The assumption of a scalar leader's decision can be generalized and is mainly taken for the ease of presentation. We will discuss this in more detail in our conclusion in Sect. 6.

Using neural networks to approximate the follower's response
Stating the bilevel problem (1) relies on the knowledge about an explicit formulation for the lower-level problem (2). In the case in which the follower's problem (2) is not known by the leader but past pairs (x, y) of leader and follower decisions are available, we propose using this historical data to learn the optimal response y = Ψ(x) of the follower using neural networks.
Such (x, y)-pairs naturally arise in many applications. For instance, for pricing models, the upper-level variable x is the price set by the leader for a certain good and y is the amount of this good bought by the follower. Both the price and the bought goods can be observed and collected to obtain (x, y)-pairs to train a neural network. Similar situations appear in many other fields such as in bilevel optimization for market design, transportation, or security. Obviously, it is required that the game modeled by the bilevel problem is played repeatedly so that large enough training sets can be collected over time.
In what follows, we use a neural network to learn Ψ(⋅) , which will then replace the lower level and turn the bilevel model into the single-level problem where g i is the function corresponding to the neural network for the ith follower's response In what follows, we exploit some recent results from the literature to show that g i is Lipschitz continuous and that Lipschitz constants can indeed be computed. This property of the used neural networks is vital for the decomposition method we use to solve Problem (3); see Sect. 4.

Lipschitz constants of neural networks
In this paper, we use the LipSDP-Neuron method published in Fazlyab et al. [11] to compute Lipschitz constants of the neural network functions g i , i ∈ [n y ] . That is, we compute a constant L i ≥ 0 that satisfies and for all i ∈ [n y ] . The main idea in Fazlyab et al. [11] is to replace the nonlinear activation functions at the nodes of a neural network by so-called incremental quadratic constraints, which then allows to state the problem of estimating Lipschitz constants as a semidefinite program (SDP). It is worth noting that the most complex and, hence, the most accurate version of LipSDP in Fazlyab et al. [11] is shown to be wrong in Pauli et al. [22]. Thus, we use the second of the three versions of LipSDP, namely LipSDP-Neuron.
We now describe the quadratic constraints that we use to replace all activation functions (x) = [ (x 1 ), … , (x n )] ⊤ ∶ ℝ n → ℝ n in a network layer, where the same sloperestricted function ∶ ℝ → ℝ is applied to each component of . Here and in what follows, we call a function slope-restricted w.r.t. 0 ≤ < < ∞ if holds for all x, y ∈ ℝ . This definition states that the slope of the line connecting any two points x and y on the graph of is bounded by and . It is easy to see that the Rectified Linear Unit (ReLU) activation function defined as (x) ∶= max{0, x} is slope-restricted with = 0 and = 1 ; see Goodfellow et al. [13]. Furthermore, if the activation function is slope-restricted w.r.t.
The following lemma shows that the slope property (4) can be written as a quadratic constraint. We use the notation n for the set of all symmetric n × n matrices.
Lemma 1 (Based on Lemma 1 in Fazlyab et al. [11]) Suppose ∶ ℝ → ℝ is sloperestricted w.r.t. and . Moreover, we define the set of diagonal matrices T with nonnegative entries. Then, for any T ∈ T n , the function The statement of the lemma above is based on Lemma 1 in Fazlyab et al. [11] and has been modified according to the corrections published in Pauli et al. [22]. There, the authors give a counterexample that shows that the original lemma in Fazlyab et al. [11] is wrong. To be more specific, they show that there exists a matrix T built according to the definition of the set T n given in Fazlyab et al. [11] that violates (6).
Assuming that all activation functions ∶ ℝ → ℝ are the same, a feed-forward neural network f (x) ∶ ℝ n 0 → ℝ +1 can be written compactly as is the concatenation of the input values at every layer of the network, is the number of layers, x i ∈ ℝ n i for all i ∈ {1, … , } , is the vector-valued function, which applies the activation function to every entry of the input vector, and holds, where W i is the weight matrix connecting layer i with layer i + 1 , and I n i is the n i × n i identity matrix. The next theorem is the central result in Fazlyab et al. [11] (5) and shows that the Lipschitz constant of a neural network is the solution of an SDP in which the matrix T defined in Lemma 1 serves as a decision variable.
Theorem 2 (Theorem 2 in Fazlyab et al. [11]) Consider an -layer and fully connected neural network given by (7). Let n = ∑ k=1 n k be the total number of hidden neurons and suppose that the activation functions are slope-restricted w.r.t. and . Define T n as in (5), A and B as in (8), and consider the matrix inequality with As a result of Theorem 2, a Lipschitz constant for multi-layer networks can be computed by solving where ( , T) ∈ ℝ + × T n are the decision variables. Furthermore, M( , T) is linear in and T and the set T n is convex. Hence, Problem (10) is a semidefinite program.

Solution approach
The neural network constraint (3c) turns model (3) into a challenging problem. In particular, for complex and large neural networks, it cannot be expected to get a closed-form expression for g i , i ∈ [n y ] , that has reasonable properties required for optimization. In any case, the resulting constraints will be nonlinear, nonconvex, and nonsmooth. However, we can evaluate these constraints and we can compute their Lipschitz constants. Thus, it is reasonable to make the following assumption.
Assumption 2 An oracle is available that evaluates g i (⋅) for all i ∈ [n y ] and all g i are globally Lipschitz continuous on x ≤ x ≤x with known global Lipschitz constant L i .
To solve Problem (3), we use the decomposition method published in Schmidt et al. [25] but modify it slightly so that we can apply it to our setting. We assume Lipschitz continuity of the problematic nonlinearities and use the corresponding Lipschitz constants to build a mixed-integer linear problem (MILP) that is a relaxation of the original model (3). We refine this relaxation in every iteration until a satisfactory solution is found or until the problem is shown to be infeasible. A satisfactory solution is formally defined to be an -feasible point, i.e., a point that solves where ≥ 0 a user-specified tolerance. Note that this relaxation of Problem (3) only affects the nonlinear functions g i , i ∈ [n y ].

Core ideas and notation
The main idea of the decomposition method is to relax the graphs of g i with help of a set Ω i . The set Ω i is given by linear constraints built around g i using the corresponding Lipschitz constant L i . The resulting relaxation Ω i is then refined in each iteration k such that

3
Using neural networks to solve linear bilevel problems with… for j ∈ J k i . Visually speaking, (12) states that the polytopes Ω k i (j) are all quadrilaterals with two vertices x k,j−1 i and x k,j i on the graph of g i ; see Fig. 1. To understand how a point is added to X i , we discuss the two problems solved in each iteration of the algorithm.

The master problem
The master problem (M(k)) is solved in each iteration over the sets Ω k i . The problem is formally given by with = (x, y) and i ∶= (x, y i ) . Note that i is not the ith entry in , since ∈ ℝ 1+n y is the vector containing x and y while i ∈ ℝ 2 contains x and y i . If the solution k with x k ∈ ℝ and y k ∈ ℝ n y of (M(k)) is -feasible, then Problem (3) is approximately solved. According to Schmidt et al. [25], the master problem (M(k)) can be modeled as the mixed-integer linear problem

The subproblem
On the other hand, if the solution k is not -feasible, the relaxation Ω k i of the graph of g i needs to be refined. This is achieved by solving a subproblem to find a point on the graph of g i , which is as close as possible to the solution k i ∈ Ω k i of the master problem (M(k)). We then use this point to refine the relaxation.
Our subproblem differs from the original method described in Schmidt et al. [25]. There, an optimization problem over nonlinearities and other linear constraints is solved. This is not possible in the setting we consider here, or is at least extremely challenging, due to the nonconvex and nonsmooth nature of neural networks.
The subproblem is solved only for the particular polytope . That is, we look for a point on the graph of g i , which is contained in Ω k i (j k i ) . Additionally, the feasible set of the subproblem is further reduced to only allow for a solution in an inner-approximation of the respective polytope. This way we ensure that newly found points do not accumulate at an already existing value in X k i . For a given j ∈ J k i , this subset is defined as

3
Using neural networks to solve linear bilevel problems with… with and d k,j is the length of the corresponding segment on the x-axis. The constant 0.25 can be replaced by any value in (0, 0.5). The left plot in Fig. 1 indicates the set Ω k i with vertical dashed lines. To solve the subproblem, we sample points x i on the x-axis segment corresponding to Ω k i (j k i ) and evaluate g i at these points. Then, the solution of the subproblem is given by the computed point ̃i ∶= (x i , g i (x i )) that is closest to the solution k i = (x k , y k i ) of the master problem w.r.t. the Euclidean distance; see Fig. 1. In other words, we choose from a finite set of points in Ω k i (j k i ) the one that is on the graph of g i and as close as possible to the solution k i of (M(k)). Note that we solve the subproblem only for those i ∈ [n y ] that are not -feasible, i.e., if |g i (x k ) − y k i | > holds. The solution x k i of the subproblem is added to the set X k i and thus refines the relaxation of g i ; see Fig. 1 (right). While x k in k i is the solution of (M(k)) and thus shared by all i ∈ [n y ] , the solution x k i of the subproblem unique to i ∈ [n y ].

Algorithm and convergence properties
The Lipschitz decomposition method is formally given in Algorithm 1. There, the master problem (M(k)) is solved in each iteration k and the algorithm checks if the solution k is -feasible. If this is not the case for all i ∈ [n y ] , the polytope Ω k i (j k i ) containing the solution k i is identified using the index of the variables z k,j i in the MILP (13) for all -infeasible i ∈ [n y ] . Then, a new point x k i is added to X k i , refining the approximation of g i . Notice that while x is scalar, the index i in x k i signals that the new point found on the x-axis in iteration k is specific to the relaxation of function g i . Moreover, the number of points in X k i is at most k in every iteration, which means that |J k i | ≤ k for all i ∈ [n y ].
We remark that the new point x k i found by the subproblem splits the original quadrilateral Ω k i (j k i ) into two smaller ones; see Fig. 1. For the two new polytopes, we use the notations Ω k i (j k 1 ) and Ω k i (j k 2 ) . Hence, the union of polytopes Ω k+1 i that approximates g i in iteration k + 1 is given by and we have the following lemma. Similarly to Theorem 1 in Schmidt et al. [25], Theorem 4 follows from the fact that, according to Lemma 3, the volume of Ω k i decreases by a positive value k in each iteration that is uniformly bounded away from zero.
However, due to the fact that Ω i is refined in each iteration to better approximate g i , the master problem (M(k)) grows linearly over the course of the iterations, leading to one additional binary variable z k,j i and some additional linear constraints. Consequently, the computational effort increases in every iteration.

3
Using neural networks to solve linear bilevel problems with…

Case study
In this section, we illustrate the applicability of Algorithm 1 on the basis of a set of exemplary case studies. We first discuss some implementation details and then present the results of the algorithm.

Implementation details
We now explain how we generate the (x, y)-pairs for the considered instances. We also discuss the used neural networks and their training, the sampling of points for the subproblem, and how we verify the solutions obtained with our algorithm. We generate the (x, y)-pairs as follows. First, we solve the high-point relaxation of the bilevel problem but with a different objective function in which we either minimize or maximize the upper-level variable x to get the set F x . Then, we equidistantly sample in this interval and solve the parametric lower-level problem for the given values of x, obtain the lower-level problem's solution y, and store the point (x, y) with y = Ψ(x) in our data set. The resulting data set entirely consists of bilevel feasible points. The obtained (x, y)-pairs are then used to train neural networks to learn the optimal follower's response. Note that using the modified high-point relaxation to obtain the x-interval used for sampling is not possible in reality if the lower-level problem is not known. However, a generation procedure would not be required in reality at all since the set of (x, y)-pairs for training the neural networks consist of observed data from the past.
The lower-level problem of every instance is then assumed to be unknown. Instead, for every considered problem, the optimal follower's response for a given leader's decision x, i.e., Ψ i (x) , is learned with feed-forward neural networks with ReLU activation functions at every node, for all i ∈ [n y ] ; see, e.g., Goodfellow et al. [13]. The functions learned with ReLU networks are inherently piecewise linear and, thus, Lipschitz continuous. Moreover, for all instances the weights of the network are updated via the Adam optimizer [24]. We vary the learning rate, the number of epochs, and the network architecture depending on the instance at hand.
The bilevel feasible set F and its projection F x onto the x-variables also depend on the unknown lower level. Hence, we attempt to deduce F and F x from the available (x, y)-pairs. Due to the fact that the interval [x,x] with x ≤ x ≤x can be larger than F x , we do not use x and x to initialize Ω i ; see (12). Instead, we use the closed interval generated by the smallest and the largest x-values in the available training set as a proxy for F x .
Given the trained networks g i , we solve the master problem and, subsequently, the subproblems as described in Sect. 4. To solve the subproblem, we evaluate the functions g i at p = 100 points on the corresponding x-axis segment. If s points have already been evaluated on this segment in earlier iterations, then only p − s new equidistantly distributed points (x, g i (x)) are computed.
Finally, Algorithm 1 stops with the indication that Problem (3) is infeasible or with an -feasible solution, where = 10 −5 is used in our experiments. The solution computed by our algorithm is compared with the solution obtained by solving the mixed-integer KKT reformulation with sufficiently large big-M constants; see, e.g., Kleinert et al. [18].
We implemented the LipSDP-Neuron method using the cvxpy 1.1.13 package and the SDP solver MOSEK 9.3.20. All occurring linear or mixed-integer linear problems (in particular, (M(k))) are solved using Gurobi 9.5.1. All neural networks have been trained using the Python library torch 1.11.0. All computations have been executed on a Intel © Core TM i7-10510U CPU with 8 cores of 1.8 GHz each and 32 GB RAM.

Discussion of the results
We apply the Lipschitz decomposition method to 6 instances of linear bilevel problems from the literature. All instances have a scalar upper-level decision variable x so that Assumption 1 is satisfied. For all 6 instances we use 60 % of the (x, y)-pairs for training and the rest for validation. Alternative proportions of training and validation set sizes could be used as well if appropriate. Furthermore, given that Ψ i (⋅), i ∈ [n y ] , is piecewise linear in our context, rather small data sets are often already enough to find good solutions for the considered instances. Table 1 shows the training set sizes, the configurations of the neural networks, and the used learning rates. Table 2 contains the corresponding Lipschitz constants computed for the networks described in the previous table as well as the time required to compute them. Finally, Table 3 displays the solutions obtained with Algorithm 1 (as well as the required number of iterations) next to the ones computed by solving the KKT reformulation. All computation times are given in seconds.

The performance of Algorithm 1
As discussed in Sect. 4.4, the computational effort of solving the master problem naturally increases over the course of the iterations. This can also be seen clearly in Fig. 2 for 4 exemplarily chosen instances. Furthermore, the computed Lipschitz constants determine the volume of the polytopes Ω k i , i ∈ [n y ] . Thus, according to Lemma 3 and Theorem 4, the convergence of Algorithm 1 directly depends on the constants L i .
One way to keep Algorithm 1 most efficient is to compute Lipschitz constants L i that are as small as possible. To illustrate this, we consider the instance [21] where Ψ(x) is the set of optimal solutions of the x-parameterized linear lower-level problem Using neural networks to solve linear bilevel problems with… Moore and Bard [21] 3.02 0.022 Table 3 Algorithm 1 solves the instances above using the networks in Table 1 and the corresponding Lipschitz constants in Table 2 Instance  The optimal solution of Problem (14) is (0, 1.5). As explained in Sect. 5.1, we use the smallest and the largest x-values in the available data set to initialize [x,x] , which is required in Line 1 of Algorithm 1. For this instance, the unknown lower level is substituted by a neural network that has 2 hidden layers with 5 nodes each and is trained on a set of 30 points. The learning rate is approximately 0.074. Using the trained network's weights, we compute a Lipschitz constant of about 3.02. Finally, Algorithm 1 finds the approximate solution (0, 1.4999) in 1.21 seconds and 31 iterations for = 10 −5 . Table 4 captures how different Lipschitz constants influence the computation time and the number of iterations of Algorithm 1 when applied to Problem (14). For all computations, we keep the tolerance = 10 −5 . This table clearly shows that the algorithm performs better, the closer the used Lipschitz constant is to the true value of 2.5.
Moreover, constants smaller than 2.5 seemingly perform even better. For instance, using 0.5, Algorithm 1 finds a solution in only 10 iterations. Nevertheless, using values smaller than the Lipschitz constant of function g i (⋅) can potentially lead to false infeasibilities reported by the algorithm.

The accuracy of Algorithm 1
According to (11), an -feasible solution of an instance is -close to the graph of all functions g i (⋅) . This does not mean that it is also necessarily close to the true optimal response Ψ i (⋅) , i ∈ [n y ] . Consequently, the accuracy of the method depends on the ability of the neural network to approximate the true optimal responses as good as possible. Figure 3 illustrates the importance of a good approximation. There you can see in the bottom plot that the -feasible solution, in this case (0.03, 1.67), obtained for instance (14), is -close to (0.03, g(0.03)) but not to (0.03, Ψ(0.03)).
In this context, it is also interesting that we observed that larger training sets can lead to better solutions. Table 5 shows how the computed Lipschitz constants and the solutions obtained for the instance in Bard [1] change with increasing data sets.

Conclusion
In many practical situations, the leader of a bilevel optimization problem is not aware of an explicit formulation of the follower's problem. To cope with this issue, we proposed a method that uses neural networks to learn the follower's optimal reaction from past bilevel solutions. After training of the network, we compute Lipschitz constants for the learned functions and use a Lipschitz decomposition method to solve the reformulated, single-level problem with neural-network constraints. Fig. 3 Two neural networks trained on the same set with different learning rates. Consequently, the learned function g (blue curves) is different and Algorithm 1 computes two different solutions (blue stars) for instance (14) using these functions Using neural networks to solve linear bilevel problems with… This short paper should serve as a proof of concept for the ideas sketched above. However, many aspects can be improved. First of all, the assumption of a scalar leader's decision seems rather strong. Very recently, a follow-up paper [14] of Schmidt et al. [25] appeared, in which a similar Lipschitz decomposition method is developed that can tackle the multi-dimensional case. This method can now be used to also consider bilevel optimization problems with an unknown follower problem and multiple decision variables of the leader. Second, we restricted ourselves in this paper to the case in which we have no coupling constraints. Fortunately, it is rather straightforward to use Algorithm 1 also for linear bilevel problems with coupling constraints as well. Even if the bilevel feasible set is disconnected due to the presence of coupling constraints, the follower's optimal response function Ψ i (⋅) is learned to be a piecewise linear function by the corresponding neural network. Due to the fact that solutions are found at vertices of the feasible set, at least one feasible point is always available to serve as solution, even if the line connecting two points is not bilevel feasible. On the other hand, this is getting more complicated if nonlinear bilevel problems are considered. Hence, the setting of nonlinear problems with coupling constraints is a topic of future research. Third, the overall idea should be reasonable also in the mixed-integer case, at least if no coupling constraints are present. Fourth, there has been some recent work [6] on the relation between multilevel mixed-integer linear optimization problems and multistage stochastic mixed-integer linear optimization problems with recourse. Hence, it might be possible to exploit these relations to carry over learning-based techniques for two-stage stochastic optimization to bilevel optimization. Fifth and finally, it would be very interesting to see an application of the concept discussed in this paper to a real-world situation, which is, however, out the of the scope of this short paper.
Acknowledgements The authors thank the DFG for their support within RTG 2126 "Algorithmic Optimization".
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission