Learning to Solve Bilevel Programs with Binary Tender
Authors: Bo Zhou
Ruiwei Jiang
Siqian Shen
Venue: ICLR 2024 poster (2024)
Keywords: Deep Learning
Bilevel Program
Binary Tender
Enhanced Sampling
Input Supermodular Neural Network
Primary area: optimization
TL;DR: We develop a enhanced sampling method and a novel input supermodular neural network to solve bilevel programs with binary tender

Abstract: Bilevel programs (BPs) find a wide range of applications in fields such as energy, transportation, and machine learning. As compared to BPs with continuous (linear/convex) optimization problems in both levels, the BPs with discrete decision variables have received much less attention, largely due to the ensuing computational intractability and the incapability of gradient-based algorithms for handling discrete optimization formulations. In this paper, we develop deep learning techniques to address this challenge. Specifically, we consider a BP with binary tender, wherein the upper and lower levels are linked via binary variables. We train a neural network to approximate the optimal value of the lower-level problem, as a function of the binary tender. Then, we obtain a single-level reformulation of the BP through a mixed-integer representation of the value function. Furthermore, we conduct a comparative analysis between two types of neural networks: general neural networks and the novel input supermodular neural networks, studying their representational capacities. To solve high-dimensional BPs, we introduce an enhanced sampling method to generate higher-quality samples and implement an iterative process to refine solutions. We demonstrate the performance of these approaches through extensive numerical experiments, whose lower-level problems are linear and mixed-integer programs, respectively.

Published as a conference paper at ICLR 2024
LEARNING TO SOLVE BILEVEL PROGRAMS WITH BINARY TENDER
Bo Zhou, Ruiwei Jiang, Siqian Shen
Department of Industrial and Operations Engineering University of Michigan, Ann Arbor, MI 48109, USA (bozum, ruiwei, s sigian) ammannmm@emsee www.ch.ede.
ABSTRACT
Bilevel programs (BPs) find wide range of applications in fields such as energy. transportation, and machine learning. As compared to BPs with continuous (lin- ear/convex) optimization problems both levels, the BPs with discrete decision variables have received much less attention, largely due to the ensuing computa- tional intractability and the incapability of gradient-bas based algorithms for handling discrete optimization formulations. In this paper, develop deep learning tech- niques to address this challenge. Specifically, we consider BP with binary ten- der, wherein the upper and lower levels are linked via binary variables. We train a neural network to approximate the optimal value of the lower-level problem, as function of the binary tender. Then, we obtain single-level reformulation of the BP through a mixed-integer representation of the value function. Furthermore, we conduct a comparative analysis between two types of neural networks: gen- eral neural networks and the novel input supermodular neural networks, studying their representational capacities. To solve high-dimensional BPs, we introduce an enhanced sampling method to generate higher-quality samples and implement an iterative process to refine solutions. We demonstrate the performance of these ap- proaches through extensive numerical experiments, whose lower-level problems are linear and mixed-integer programs, respectively.
INTRODUCTION
Bilevel programs (BPs) are appealing for modeling problems that involve sequential decision inter- actions from two or multiple players. Their hierarchical decision processes arise in a wide range of real-world applications, including energy (Zhou et al., 2023), security (Bhuiyan al., 2021), transportation (Santos et al., 2021), and market design (Nasiri et al., 2020). Recently, BPs have also shown strong modeling power machine learning problems, including meta learning (Wang al., 2021), actor-critic reinforcement learning (Hong al., 2020), hyperparameter optimization (Bao et al., 2021), and deep learning (Chen et al., 2022). Next, we provide a brief introduction to BP, its generic formulation, and solution methods.
Bilevel Program. In a BP, leader and a follower solve their own decision making problems in an interactive way: The leader's s decisions made in the upper level will affect the follower's problem solved at the lower level (e.g., the leader's decisions are involved in the objective function and/or constraints of the follower's problem) and vice versa. This can be described formally as
min f(r, (*)
(1a)
E X(v")
(1b)
where y" E arg max g x)
(1c)
y Y(x),
(Id)
where æ represents the leader's decision and y represents the follower's. Here, the leader's decision
x will affect the follower's 's objective function g(y. x) and feasible region Y(x) C IR.", On the other
hand, the leader's 's objective function fxx y²) depends on both her own decision and the follower's
optimal decision y", which a function of x defined through (1c)-(ld). Additionally, the leader's
feasible region X ((v*) C RB can also depend on y"

---

Published as a conference paper at ICLR 2024
Single-Level Reformulation. To solve the BP (1), general approach is incorporate the lower- level decisions and their feasible region into the upper-level problem, giving rise to
min fx,
(2a)
s.t. (())
(2b)
y (E)()..
(2c)
Formulation (2) relaxes the optimality of y, and therefore, its solution only provides a lower bound for the BP (1). To retrieve optimality, one can incorporate an optimality condition
gly d(x),
where (xx := (x) (a)) g(y, 2)) represents the optimal value of the lower-level problem as a function of z. Hence, combining (2) and (3) produces single-level reformulation of the BP (1). In case the lower-level problem is continuous (i.c., y consists of continuous decision variables only), constraint (3) be made explieit using strong duality or the KKT conditions. However, such luxury is immediately lost if the lower-level problem involves discrete decision variables (e.g.. is binary or mixed-binary). In that case, the closed-form expression of ((x) is either non-existent or highly intractable, prohibiting solving the BP effectively.
In light of this challenge, we study BPs with binary tender. Here, "tender" is defined as the linking variables between the upper- and lower-level problems and "binary tender" means that all tender variables are binary. Note that the problems in both levels can involve general decisions (e.g,, continuous and/or integer variables), and we only assume that the entries of x appearing in the lower- level formulation are binary-valued. Such BPs arise many applications, including energy system expansion planning (Kabirifar et al., 2022), charging station planning (Li et al., 2022), competitive facility location (Qi et al., 2022), and network interdiction (Smith & Song, 2020). In addition, we assume the leader's feasible region is independent of y² i.e., X (v²)) X, which is special case of (I). this case, we can obtain an upper bound for the BP (1) from any feasible solution x and the corresponding follower's 's optimal solution y" (e.g.. by solving the follower's 's problem using z). Our main contributions are three-fold.
We employ neural networks to learn and approximate the value function 6(x). Then, derive a closed-form representation of the learned value function. This yields a single-lvvel, mixed-integer formulation for the BP (I), which be readily solved by off-the-shelf optimization solvers.
Motivated by the fact that (xx) supermodular for a large class of BPs, we design an input supermodular neural network (ISNN) that ensures supermodular mapping from input to output. We analyze the representability of both general neural networks (GNN) and ISNN provide guidance for network architecture selection.
To solve high-dim -itbe-iimnsional BPs, we propose an enhanced sampling method for generating higher-quality samples and training neural networks to better approximate ((x) around the optimal solution. Building upon this enhanced sampling method, we execute an iteration process to improve the accuracy of the derived solution.
The remainder of the paper is organized as follows. In Section 2, review related works in the lit- crature. Section 3, we elaborate our methodology of learning to solve BPs. We conduct numerical experiments in Section 4 and draw conclusions in Section 5.
2 RELATED WORKS
The theories and algorithms for finding optimal solutions to BPs depend on the structure and proper- ties of the lower-level problem, as well as the coupling between the two levels (Kleinert al., 2021; Beck et al., 2022).
Bilevel linear/convex programs. The earlier studies on BPs focused on linear or convex lower- level problems. Consequently, one can use the KKT conditions of the lower-level problem or strong duality to reformulate the BP as a single-level problem with complementarity constraints (Fortuny- Amat & McCarl, 1981) or bilinear terms (Zare et al., 2019; McCormick, 1976). respectively. In both cases, solving BPs boils down to solving the resulting single-level, nonconvex, and nonlinear reformulation (Colson et al., 2005; Bard, 1998).
2

---

Published as a conference paper at ICLR 2024
Bilevel (mixed-integer) nonconvex programs. Recently, more studies focused BPs having discrete decision variables or nonconvex objectives/consrrinns. In such problems, we notice that, because of the existence of integer decision variables or nonconvexity, neither KKT conditions nor strong duality approaches may be able to capture the (parametric) global optimal so- lutions at the lower level. To achieve global optimum, instead need to exploit special structures of BPs (Qi et al., 2022) or construct the aforementioned relaxation (2) for single-level reformu- lation. For BPs with nonconvex objectives/constrains. approximation methods were employed to approximate the lower-level problem and reformulate the BP as single-level formulation, such as polyhedral approximation (Sato al., 2021) and other gradient-biae baddentt-baed methods in machine learning (Hong et al., 2020; Chen et al., 2022). However, when discrete decision variables are incorporated, gradient-based methods are no longer applicable. To address this issue, integer programming tech- niques are employed and various cutting plane-based algorithms, such as the MiBS solver, have been developed (Ralphs et al., 2015; Zeng & An, 2014). However, learning techniques have not been explored in solving such BPs, which is the focus of this paper.
general
Embedding neural networks into optimization. Due to its strong power in data analysis and fit- ting, machine learning, such as neural networks, has drawn wide interest in function approximation (Fujimoto et al., 2018; Liang & Srikant, 2017; Ferrari & Stengel, 2005). Furthermore, when the adopted activation function is a piecewise linear function, such as ReI.U, the output ((x) of a neural network can be reformulated as mixed-integer linear form (Fischetti & Jo, 2018; Serra et al., 2018).
Consequently, in case ((x) unknown hard access, one can use o(xx) to replace ((x) and di- rectly incorporate the mixed-integer representation of ((z) into optimization formulations. Molan & Schmidt (2022) considered this idea in BPs with unknown lower-level problems. They utilized a neural network approximate the mapping from x to the optimal y" the lower level, and ex- ploited Lipschitz optimization techniques reformulate the activation function (ReLU). Different from Molan & Schmidt (2022), we propose to learn and approximate ((x) in BPs.
Input convex neural network (ICNN). Another related stream of literature proposed ICNNs, which ensure that the output approximation @(x) is convex function. In this case, the target optimization formulation remains convex program if it was before incorporating (zx). ICNN was first pro- posed by Amos et al. (2017), who derived sufficient conditions on the neural network architecture and parameter settings for the convexity of ((x). ICNNs were then applied to, e.g., optimal control (Chen al., 2019) and energy optimization (Bünning et al., 2021). Different from these works, we study BPs with binary tender, and thus the input æ of the neural network is binary-valued. Ac- cordingly, consider supermodularity. the counterpart of convexity discrete optimization, and establish ISNN that guarantees output a supermodular ((x). In Section we derive neural network architecture for ISNN and show its representability.
METHODOLOGGY
We neural networks (see Fajemisin al. (2023)) to obtain an approximate value function 5(x)) of <(z) and solve the BP (1) by incorporating a closed-form expression of 6(x) into formulation (2)-
(3). Figure 1 illustrates the proposed method, including sampling, training, and solving.
Sampltng (Section 3.1)
Samples
10
0o
#00 2 ((rw)
CHOCKEI
(
Solving (Section 3.3)
Closad-form reprisontaioon
Training (Section 3.2)
min (x,y) s.t.x X
Xg
y (x))
300
#(v.x) 2 ((x)
O(*)
min (8)
#Co1
Figure 1: The overview of the proposed method.
3

---

Published as a conference paper at ICLR 2024
3.1 SAMPLING
We consider set of the upper-level decision variables á and compute d()) = maxy {gly, #) y ( to obtain sample-label pairs (i, d())). When the dimension n of x is small, we can enumerate all possible sample-label pairs. As n increases, the number of all sample-label pairs increases exponentially and enumeration becomes numerically prohibitive. In this case, we can sample sufficiently large set of & instead. However, a naive sampling may frequently produce infeasible solutions, let alone obtaining high-quality ones. To this end, we propose an enhanced sampling method to quickly find feasible samples and improve the quality of samples when the number samples is limited.
Enhanced Sampling. A basic idea in enhanced sampling is "sampling via optimization." This involves confining the sample selection within the feasible region of the BP with respect a random but tractable objective function. This is mathematically specified as
min x Qu + h
(4a)
s.t. x E X, E
(4b)
where Q and are randomly generated n X n positive semi-definite matrix and n x 1 vector, respec- tively. Accordingly, we efficiently generate samples by solving (4). We note that the reason for solving quadratic program here is to avoid repeated samples and enhance sampling efficiency.
Another basic idea is to sample more often in the vicinity of the optimal solution the BP. A more
accurate approximation Φ(x) to the true ((x) in this vicinity, even at the cost of inaccuracy for those x far from optimum, enhances the tightness of approximating the BP (1), because most of the
effort in solving (1) is spent on comparing the near-optimal æ's. To this end, we let fab represent a known upper bound of the BP (1), which can be obtained from any feasible solution. Then, we strengthen (4) be
min x xQQr= + hTxx s.t. f(z, y) S fub
(5a)
(5b)
E X, E (()),
(5c)
where constraint (5b) restricts the search space to elevate quality of samples. Here we note that though (5) is mixed-integer quadratic program, we only need a feasible solution to (5) for sampling, which does not incur too much computational burden and saves the overall computational time.
Furthermore, by sampling feasible solutions & to the BP, we can compute their objective value f (²"), 8" representing the corresponding optimal solution from the lower level, and iteratively strengthen fub if f(i 5* < fub- This allows to continuously refine the search space and obtain higher-quality samples. Algorithm summarizes the above sampling process.
Algorithm Enhanced Sampling
1: INPUT: sample size N., initial upper bound fui, maximum number of updates Nub-
2: Initialize the set of samples Ω <- 2, count of samples k2 <- 0, count of updates kigb 4- 0.
3: repeat
4: Randomly generate an x n positive semi-definite matrix Q and an n X vector h.
5: Solve (5a)-(5c) and store an optimal solution a.
6: if â & Ω then
7:
Solve maX) ( (g.b) .). and store the optimal value ((+) and an optimal solution g". Augment + Q {á] and k, ke
8:
9:
end if
if
kibb < Nus and f(t 0*) < fub then Update fab & f(i
11:
kig + Rig + 1. end if
13:
14:
until hi No
15:
OUTPUT: sample- label pairs (z, ((x)) for all Ω.
3.2
TRAINING
Using the sample-label pairs from Algorithm 1, we adopt supervised learning to train a neural net- work to fit the mapping (xx:). In this step, we consider two types of neural networks, i.e., GNN and
4

---

Published as a conference paper at ICLR 2024
ISNN, which use the same architecture shown in Figure 2. There are K hidden layers and one out- put layer in the architecture and we employ passthrough to enhance the representability of the neural network. We note that the input of the architecture, 2, be different from (see Sections 3.2.1- 3.2.2 for details). In addition, & denotes the output of the neural network, 2 denotes the output of the kth hidden layer, / D& and by are the weights and biases of the kth layer, respectively, and at)) denotes the activation function.
Weights
We Dv
Waut
Activation
D2
Dzss
Output
Impott
Outputt of hickem bayers
Figure 2: Neural network architecture.
Specifically, the neural network defines ((z) through
o (W1 2 + b1)
(6a)
24 <W + by D&t), Ak 2...., K,
(6b)
3 WK+IZK + bKCH + DK+1I.
(6c)
3.2.1 GENERAL NEURAL NETWORK
For GNN, we define r. Proposition 1 shows the representability of the GNN.
Proposition 1. Consider an arbitrary é {0, 13" + RR. Then, the following statements hold. i. (Universal Approximation) there exists GNN with the architecture in Figure such that o fits á exactly, Le., ((x) (xx) for all E 1)
ii. (Maximum Representability) Suppose that the hidden layers of the GNN consist of Nar many neurons, which may be distributed arbitrarily among these layers. If we measure the GNN repre- sentability by the number of parameters trains, then the representability maximized when GNN consists of two hidden layers, i.e., when K = 2. In particular, the first layer consists of Nur / (if Nar is even) or (Nnr - 11/22 neurons (if Nar is odd), and the second layer consists of Nur/2 (if Nur is even) or (Nar + 1)/2 neurons (if Nur is odd).
ii. (Sufficient Fitting) Given N, sample-label pairs, the GNN able to output an approximation
that fits these pairs exactly if
Nur
V(2n 2)2 4N, + -2n 3)
We give a proof of Proposition in Appendix A.I. According to Proposition 1, when we adopt GNN for training, incorporate hidden layers the GNN architecture, each with Nnr/2 neurons, where Nar is the smallest even number that satisfies (7).
3.2.2 INPUT SUPERMODULAR NEURAL NETWORK
In light of the special class of BPs studied in this paper, we are motivated design an ISNN ar- chitecture that guarantees to output a supermodular approximation Φ(ă) of the true value function ((x), where we set ह : [xr (1 x)) TJT in the case of ISNN. Long et al. (2023) and Chen et al.
(2021) provide some sufficient conditions for parametric optimization ((x) to be supermodular.
First, we recall the definition of supermodularity.
Definition Consider subxet D of R", a function f D + R. and for x. E D define V y [max(f1, V1 t.)..... max En, yo T1T and x Ay= [min v1). vi)..... min En, VollT. Then, is called supermodular if f(x) f(y) < f(x + A for all x, y.y E D. In addition, a set S is called lattice if V E S and x Ay E for all x, S.
mixed- valued argument and the domain D of a supermodular function f needs not to be (0, 13
We remark that this definition of supermodularity applies
Then, ready present the main result of this section.

---

Published as a conference paper at ICLR 2024
Proposition 3. The function ((3) output from ISNN is supermodular in E if W W1((+11) and D2.K are non-negalive, and a(.) is convex and non-decreasing (e.g., ReLUU.
The proof is given Appendix A.2. By Proposition 3, we simply need to add appropriate sign constraints to some weight parameters when training the ISNN.
Next, Proposition shows the representability of ISNN.
Proposition 4. Consider arbitrary (not necessarily supermodular) Φ {0,1)" R. Then, the following statements hold.
i. (Universal Approximation) there exists ISNN with the architecture in Figure 2 such that 6(x.1 1 2) )(x) for all (0, 1)".
(Maximum Representabiltty) With fixed number of neurons Nar in its hidden layers, the ISNN's architecture does not affect its representability (measured by the number of parameters the ISNN trains).
ii. (Sufficient Fitting) Given Ng sample-laber pairs, the ISNN able to output an approximation
that fits these pairs exactly if
N 2n +
Nur 2
(8)
We give proof of Proposition Appendix A.3. Since the ISNN's architecture does affect its representability. to better compare with GNN, we incorporate two hidden layers in ISNN, each with Nur/2 neurons, where Nzr is the smallest integer that satisfies (8).
Thanks to supermodularity of δ(ā) promised by Proposition we can recast the approximate optimality condition g x) 2 (()) as linear inequalities, in lieu of linearizing 6(ī) using auxiliary binary variables and big-M coefficients. We place the theoretical analysis for supermodular ((7) in Appendix B.
3.3
SOLVING
3.3.1
MIXED-INTEGER REPRESENTATION OF ((2)
We represent the epigraph of o(i)) as linear inequalities using auxiliary binary variables, to be em- bedded into the BP reformulation (2)- To this end, we recall that the ReLU activation function takes the form
x, x 00
(x) max 03 =
(9)
0, <
Then, (zx) can be rewritten as
0 < (xx < M8
(10a)
x S (x) S x + M(1-5)
(10b)
where 5 is an auxiliary binary variable and M is a sufficiently large positive number. We note that, in computation, the value of M does not have be arbitrarily big and it can be iteratively strengthened by solving relaxations of (2)-(3); see, e.g., Qiu et al. 2014. Likewise, the defition of ((i) in (6) can be represented layer by layer as follows:
of W bj
0 2 $ M8
(11a)
o s 2 S + M(1 81)
1 Wh + by. + DAI
0 2 $ M5x
Vk = 2, 2...,
(11b)
2 S 2 $ 2 + M(1 5k))
WK+ + bx+1 + Dx+1ū,
(1lc)
where 2t and 62 auxiliary continuous and binary variables for all 1 < k K, respectively. Plug- ging this representation into (3) yields a mixed-integer programming approximation of the BP (1):

---

end for Published as a conference paper at ICLR 2024
min f(z,8) v) 2.3
(12a)
s.t. x E X, y (x))
(12b)
gly, r) 2 @(0). (11)
(12c)
2 if using GNN.
(12d)
[x
[ if using ISNN.
Here we note that existing algorithms and highly-efficiee off-the-shelf solvers can be adopted to handle the single-level program (12), which is yet out of the scope of this paper.
3.3.2 NEURAL BILEVEL ALGORITHM
Combining the above sampling, training, and solving processes, we conclude our neural bilevel algorithm in Algorithm 2. By iteratively conducting the enhanced sampling, we can find new and higher-quality samples and improve the accuracy of the found solution.
Algorithm Neural Bilevel Algorithm
1: INPUT: maximum number of iterations Nueration, initial upper bound fub := +o0.
2: Solve (2) and store an optimal solution x=.
3: Store y* E argmax, (xY()*) g(,xz))
4: Update fub +- min .y)".
5: for i = -iierrrvttttoooo do
6: Conduct Algorithm to obtain the pairs (④, (()) and their corresponding f(2, 9).
7: Train a neural network 2(z) (GNN ISNN) using (#, (())).
8: Solve (12) using the trained e and store an optimal solution x².
9: Store a y" E arggaaxyy (Vv
10: Update Sus <- min(Jus, [(x", ,y" ((,), ((,,0))
12: OUTPUT: Current (i.e., best) upper bound fub and its corresponding solution (I, y).
EXPERIMENTS
4.1
ILLUSTRATIVE EXAMPLES
We first illustrative example to show the effectiveness of our proposed methods. We consider a BP with a nonconvex and nonlinear program in the lower level:
min 2 1 + x2 3y
s.t. 11,42 E {0,
(13a)
where y² E arg max 2)2
s.t. 0 - - 1 2 x22
(13b)
This BP admits the following optimality condition for the lower-level problem:
(y 2)² 6(11,
(14)
where 22) max(-(v 2)2 (13b)) In this example, (())))2) admits closed- form expression. Indeed, by the structure of the lower-level objective function, we notice that y" = min(1 + 2|x - x22, 2} and ő(z1, x2) (min 2|x1 x2), 2} 2)2, which is de- picted in Figure 3(a). Incorporation of this expression yields the optimal solution (xi, = (0, 1). In Figure 8(a) of Appendix D.I, we illustrate the feasible region of (I1, x2, 22,8) y) described by con- straints (13a)-(13 and why (14) ensures optimality.
On the other hand, following the steps and methods described in Section 3, we use neural networks to approximate the value function ((x1, Thanks to the simplicity of this example, we find
closed-form expressions of these approximations. (0( 22) = I1 + I2 1 1- 20(1; + 22 1) for GNN and @1s(21, x2) x1 (1 x2) + 2a((1 mn) + D2 1)) for ISNN, which are

---

Published as a conference paper at ICLR 2024
04
(a) Exact ((21, 22).
(b) GNN (11,
(c) ISNN 0rs(21, T2))
Figure 3: Value function.
depicted in Figures 3(b) and 3(c), respectively. Figures 8(b) and 8(c) in Appendix D.1 visualize their corresponding optimality cuts. Both approximations de (Z1, x2) and (21, lead the (true) optimal solution x² = (0, 1)" From these, we observe that incorporating the approximate
((2), 22) from either GNN or ISNN correctly produces the optimal solution.
4.2
RANDOMLY GENERATED INSTANCES
4.2.1 SETUP
We randomly generate classes of instances for numerical experiments, with n 10, 20,30, 40, 50, 60. The details of instance generation is reported in Appendix C. We allow the lower-level problem be linear program (I.P) or a mixed-integer linear program (MILP). By de- fault, we generate 1,000 samples using Algorithm I and use them in training neural networks. We adopt ReLU as the activation function and design the architecture neural networks by Propo- sition and by Propositions 3-4 for GNN and ISNN, respectively. We the Adam algorithm (Kingma & Ba, 2014) for training for 1000 epochs and set the learning rate as 0.001 with the decay 0.001.
4.2.2 COMPARISON WITH STATE-OF-THE-ARR SOLVERS
We compare method with the state-of-the-an solver for bilevel problems, MiBS (Tahernejad et al., 2020) and use its solution within I-hour time limit as the benchmark to calculate the objective differences (i.e., the gap between the objectives from method and MiBS). The negative objective difference means that our method provides a better solution than the benchmark. For each instance, replicate method for 10 times, considering the randomness in sampling and training.
For the instances with a LP lower level, we report the objective difference of the best upper bound found by Algorithm in Figure 4 and the average computational time in Figure 5. In all figures, the legend, for example, "GNN-2" means that we use GNN to fit samples and Niceatioon = 2. Figure 4(a) and Figure 4(b) show the average and minimum objective difference in the 10 replications. The distribution of the objective difference of the 10 replications is provided in Figure 9 of Appendix D.2. The computational time reported in Figure is the average time of the 10 replications. The itemized time in sampling, training, and solving is provided in Figure 10 of Appendix D.2.
--NN-1 GAN-2 -
G-W-1 GNN-2 GN-
SAN-11 - B-N-2 -4-8AN-3
IAN-1
+ 1SN-22 - --N---
3
20
20
10
(a) Average objective difference of 10 replications.
(b) Minimum objective difference in 10 replications.
Figure 4: objective difference in different instances with linear lower level.
8

---

Published as a conference paper at ICLR 2024
2000
2500
MEBS
GAN-1
1.500
2000
GAN 2
188N-1
-G-ANN-
JSNN
18NN -3
ISNN-3
2000
500
20
30
Figure 5: Computational time in different in-
Figure 6: Computational time in different in-
stances with a LP lower level.
stances with a MILP lower level.
20
30
dbjettve
30
-
ISNN-2 1-NN-3
(a) Average objective difference of 10 replications.
(b) Minimum objective difference of 10 replications.
Figure 7: objective difference in different instances with a MILP lower level.
From the results, solution in all replications. This because we can enumerate all feasible solutions x in sampling for 10. Yet during the sampling process, there exists much repeated sampling, for which the computational time is slightly longer than that of n 20. When n 2 20, the computational time increases gets larger, which mainly results from the longer sampling time, and is almost linear with Niteration- However, the computational time is significantly shorter than that of MiBS when 40. In addition, the objective difference of instances with n 30. 50, 60 larger than 15%. The increase of Niteration results in the decrease of average objective difference for both GNN and ISNN. means that the influence of sampling gets significant and validates that enhanced sampling helps reduce average objective difference. Yet as increases, the marginal improvement of enhanced sampling becomes smaller, which is caused by the increasing difficulty finding high- quality samples. Comparing GNN and ISNN with the same Niteration. we observe that ISNN always outperforms GNN in both average and minimum objective difference. Using ISNN can reduce objective difference by more than 10% when 40, which validates the effectiveness of ISNN in improving accuracy.
see that when 10, both GNN and ISNN produce true optimal
For the instances with MILP lower level, we report the objective difference of the best upper bound found by Algorithm in Figure and the average computational time in Figure 6. More details provided in Figures 11 11-12 of Appendix D.2. We observe that these results produce similar insights as those of the LP lower-level problems. Notably, our method has excellent accuracy in most instances (achieving objective difference less than 5%, only except when = 30) and outperforms MiBS when n = 20, 40, 50. This is because when the lower level is a MILP, it becomes significantly more challenging
find high-quality feasible solutions, and these instances even MiBS reports incorrect optimal solutions. In contrast, the proposed method is able to produce better feasible solutions in (dramatically) shorter computational time, validating effectiveness.
CONCLUSIONS
We considered machine learning methods for solving -ittger, nonconvex BPs with binary tender. We developed an enhanced sampling algorithm to find high-quality samples, designed GNN and ISNN approximate the lower-level value function (in terms of the upper-level decisions), and incorporated the results optimality cuts into single-level reformulation of the BP. We validate the effectiveness the proposed approaches using an illustrative example and larger-scale, randomly generated instances. Through these experiments, reduce average objective difference and ISNN always outperforms GNN. The computational time of using either GNN or ISNN is significantly shorter than that of (a oftte-of-the--o art solver) MiBS when the dimension of the bianry tender exceeds 30. Notably, in most instances with a MILP lower level, the proposed method produce better solutions than MiBS in shorter computational time.
demonstrated that the enhanced sampling helps

---

Adejuyigbe O. Fajemisin, Donato Maragno, and Dick den Hertog. Optimization with constraint learning: A framework and survey. European Journal of Operational Research, 2023. doi: ott.://o.orgr0..0001eeoooo.00.00 Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
The authors would like to thank anonymous reviewers for their detailed and helpful comments. Ruiwei Jiang is supported in part by the U.S. Air Force Office of Scientific Research under the grant FA95550-231100233.
REFERENCES
Brandon Amos, Lei Xu, and Zico Kolter. Input convex neural networks. In International Confer- ence on Machine Learning, pp. 146-155. PMLR, 2017.
Fan Bao, Guoqiang Wu, Chongxuan Li, Jun Zhu, and Bo Zhang. Stability and generalization of bilevel programming in hyperparameter optimization. Advances in Neural Information Process- ing Systems, 34:4529 4541, 2021.
J. F. Bard. Practical Bilevel Optimization: Algorithms and Applications. Kluwer Academic Pub- lishers, Norwell, MA, 1998.
Yasmine Beck, Ivana Ljubić, and Martin Schmidt. A survey on bilevel optimization under uncer- tainty. In preparation, 2022.
Tanveer Hossain Bhuiyan, Hugh R Medal, Apurba K Nandi, and Mahantesh Halappanavar. Risk- averse bi-level stochastic network interdiction model for cyber-security risk management. Inter- national Journal of Critical Infrastructure Protection, 32:100408, 2021.
Felix Bünning, Adrian Schalbetter, Ahmed Aboudonia, Mathias Hudoba de Badyn, Philipp Heer, and John Lygeros. Input convex neural networks for building MPC. In Proceedings of the 3rd Conference on Learning for Dynamics and Control, pp. 251-262, 2021.
Can Chen, Xi Chen, Chen Ma, Zixuan Liu, and Xue Liu. Gradient Bradient-based bi-level optimization for deep learning: A survey. arXiv preprint 2rXiv:2207.117199 2022.
Xin Chen, Daniel Zhuoyu Long, and Jin Qi. Preservation supermodularity in parametric opti- mization: Necessary and sufficient conditions constraint structures. Operations research, 69 (1):1-12, 2021. doi: 10. 10. wtt.1288pr..2220.11
Yize Chen, Yuanyuan Shi, and Baosen Zhang. Optimal control via neural networks: A convex approach. In International Conference on Learning Representations, 2019.
B. Colson, P. Marcotte, and G. Savard. Bilevel programming: A survey. 4OR: A Quarterly Journal of Operations Research, 3(2):87-107, 2005.
William Cook, Albertus MH Gerards, Alexander Schrijver, and Êvva Tardos. Sensitivity theorems in
integer linear programming. Mathematical Programming. 34:251-264 1986.
Ferrari and R.F. Stengel. Smooth function approximation using neural networks. JEEE Transac- fions on Neural Networks, 16(1):24-38, 2005.
Matteo Fischetti and Jason Jo. Deep neural networks and mixed integer linear optimization. Con- straints, 23(():296-300) 2018.
J. Fortuny-Amat and B. McCarl. A representation and economic interpretation of two-level pro- gramming problem. 32(9)783-772, 1981.
Scott Fujimoto, Herke Hoof, and David Meger. Addressing function approximation enror in actor-critic methods. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 1587-1596, 2018.
10

---

Published as a conference paper at ICLR 2024
Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A two-timescale framework
for bilevel optimization: Complexity analysis and application to actor-critic. arXiv preprint arXiv:2007.05170, 2020.
Milad Kabirifar, Mahmud Fotuhi-Firuzabad, Moein Moeini-Aghtaie, Niloofar Pourghaderi, and Payman Dehghanian. A bi-level framework for expansion planning in active power distribution networks. IEEE Transactions on Power Systems, 37(4):2639-26535 2022. doi: 10.1109/TPWRS
2021.110339.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint TrXiv::412.6980 2014.
Thomas Kleinert, Martine Labbé, Ivana Ljubić, and Martin Schmidt. A survey on mixed-integer programming techniques in bilevel optimization. EURO Journal Computational Optimization, 9:100007, 2021.
Jiale Li, Zhenbo Liu, and Xuefei Wang. Public charging station localization and route planning of electric vehicles considering the operational strategy: A bi-level optimizing approach. Sustainable Cities and Society, 87:104153, 2022.
Shiyu Liang and Srikant. Why deep neural networks for function approximation? In International Conference on Learning Representations, 2017.
Ivana Ljubić and Eduardo Moreno. Outer approximation and submodular cuts for maximum capture facility location problems with random utilitics. European Journal of Operational Research, 266 (1):46-56, 2018
Daniel Zhuoyu Long, Jin Qi, and Aiqi Zhang. Supermodularity in two-stage distributionally robust optimization. MANAGEMENT SCIENCE, 2023. doi: 10. 10.1287/mmss:2023.47777
Garth P McCormick. Computability of global solutions to factorable nonconvex programs: Part iconvex undercstimating problems. Mathematical programming, 10(1):147-175, 1976.
Ioana Molan and Martin Schmidt. Using neural networks solve linear bilevel problems with unknown lower level. Technical report, Tech. rep. url: http://www. optimizationonline org/DB HTMBHHHMMM/2222/7/777777 html, 2022.
Nima Nasiri, Ahmad Sadeghi Yazdankhah, Mohammad Amin Mirzaci, Abdolah Loni, Behnam Mohammadi-Ivalle Kazem Zare, and Mousa Marzband. A bi-level market-clearing for coordi- systems in presence of energy storage technologies. Sustainable Cities and Society, 63:102439, 2020.
nated regional-local
George L. Nemhauser and Laurence A Wolsey. Maximizing submodular set functions: formulations and analysis of algorithms. Studies on Graphs and Discrete Programming, 11:279-301, 1981.
Mingyao Qi, Ruiwei Jiang, and Siqian Shen. Sequential competitive facility location: Exact and approximate algorithms. To appear in Operations Research, https: /ttp://////....111112117777 opre. 2022 2022.
Feng Qiu, Shabbir Ahmed, Santanu S Dey, and Laurence A Wolsey. Covering linear programming with violations. INFORMS Journal on Computing, 26(3):531-546. 2014.
Ted Ralphs, Tahernajad, S DeNegre, M Güzelsoy, and A Hassanzadeh. Bilevel integer optimiza- tion: Theory and algorithms. In International Symposium on Mathematical Programming, pp. 184, 2015.
Maria João Santos, Eduardo Curcio, Pedro Amorim, Margarida Carvalho, and Alexandra Marques. A bilevel approach for the collaborative transportation planning problem. International Journal of Production Economics, 233:108004: 2021.
Ryo Sato, Mirai Tanaka, and Akiko Takeda. A gradient method for multilevel optimiza- tion. M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 7522-7533. Cur- ran Associates, Inc., 2021. URL ht htt.p//pw......oo///////eeeeeeo.
files/

---

Published as a conference paper at ICLR 2024
Thiago Serra, Christian Tjandraatmadja, and Srikumar Ramalingam. Bounding and counting linear regions of deep neural networks. In Proceedings of the 35th International Conference on Machine Learning. pp. 4558-4566, 2018.
David Simchi-Levi, Xin Chen, Julien Bramel, al. The logic of logistics. Theory, algorithms, and applications for logistics and supply chain management, 2005.
J Cole Smith and Yongjia Song. A survey of network interdiction models and algorithms. European Journal of Operational Research, 283(3):797-811, 2020.
Sahar Tahernejad, Ted Ralphs, and Scott DeNegre. A branch-and-cut algorithm for mixed integer bilevel linear optimization problems and its implementation. Mathematical Programming Computation, 12(4):529-55) -081229596--11 2020.
Haoxiang Wang, Han Zhao, and Bo Li. Bridging multi-task learning and meta-learning: Towards efficient training and effective adaptation. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pp. 10991-1100) 2021.
M Hosein Zare, Juan S Borrero, Bo Zeng, and Oleg A Prokopyev. A note on linearized reformu- lations for a class of bilevel linear integer problems. Annals of Operations Research, 272(1):
99-117, 2019.
Bo Zeng and Yu An. Solving bilevel mixed integer program by reformulations and decomposition. Optimization online, pp. 1-34, 2014.
Bo Zhou, Jiakun Fang, Xiaomeng Ai, Shichang Cui, Wei Yao, Zhe Chen, and Jinyu Wen. Storage right-based hybrid discrete-time and continuous-timg flexibility trading between energy storage station and renewable power plants. IEEE Transactions on Sustainable Energy, (4(01))465-488) 2023.
12

---

Published as a conference paper at ICLR 2024
A PROOFS
A.1 PROOF OF PROPOSITION 1
Any arbitrary o: {0, 1) R can be rewritten as
d(x) (1- - -1--)
+ 1 - -3))3) ...
(11-1))a2(1 (00-))))))011-01011101100....
+ + - mz) ... ... (1,1,,)),,,00....0))
(1-21)(1--)) 000
+...
+ 22222'3--0 ". 1n 1.....1....
Σ ((2) =E(0,1)
(1 24) + (2z, 1)x:
This shows that & a polynomial of (at most) degree n. Then, & admits the following representation:
<(x)
II 2
Σ by II
Σ
Ck.
II x x
iデj
ifj.
where b, are coefficients of the polynomial. Since E {0, 10,1)" we can rewrite
II r max S
Σ
x. - (n (n (51 o Σ σ; (n - |s| -S1
ES
igS
for subsets [n]. It follows that
((x)= Σ 24 (n (1-11) + Σ byo
Ix--(n- (n 2
i/1
Σ (A.))(m) E
Ckjo
Σ
- (n - 3) +...
There are 2n terms in this expression and the first-order terms and the zeroth-order term do need activation. Therefore, (x)) can be represented using a neural network with an architecture as in Figure 2 with K 1 and 2n n 1 neurons (activations). This proves claim (i).
For GNN shown in Figure d(x) admits closed-form expression in (6). Let 2k R for all & e [K]. then the number neurons in the hidden layers Nur = Ek=1 ny and the total number of parameters Npm GNN trains is
K +1 mm + (nx + 1 + m)
N pm = (n +
k=2
T1 + (n+ + 1)
n +
K
7g + + + (n + 1) +
Then, with fixed Nar. maximize Npm by adjusting the GNN architecture boils down optimizing Σkz 7- + n over (integer) variables and n2, nx. To this end, we state the following technical lemma, whose proof is relegated to Appendix A.1.1.
13

---

Published as a conference paper at ICLR 2024
Lemma 5. Consider the following (non-convex) quadratic program with continuous decision vari- ables:
K
max
7 + 7
(15a)
8.1.
Nn
(15b)
k=1
Vn E [K]. Then, an optimal solution this program ng 0 for all k E [K 2), nk-1 (NNr 1)/2, and
Wk 0
(15c)
a = (Nar + 1)/2. Additionally, the corresponding optimal value is (Nur + 1)2//4
By Lemma 5, the optimal solution to formulation (15) remains the same if Nor is odd and we restrict each ng to take integer values only, because (NNr - 1)/22, (Nvr + 1)/2 E N.
In case that Nur is even, consider the (integer) solution nik = 0 for all A E [K 1K2 -x Nnr/22 and nı / The objective value of this solution is (N² + 2N) that is, 1/4 less than the optimal value of (15). Since all coefficients of the objective function (15a) the objective value of any (integer) solution must be integer-valued, too. It follows that there does not exist any other (integer) solution, whose objective value is strictly larger than (N2 + 2N /4. Therefore, the solution n& = 0 for all & E [K - 2]. fk-1 = Nur/2 /2, and Nx- = Nnr/2 is optimal when N is even. This proves claim (ii).
integer- valued,
To fit the N, sample-label pairs exactly, we need the number of parameters Npm of the GNN to be least large No, i... Npm Ng. But claim (ii) implies that, when choosing nk-1 7kK Nnr/2,
1 [ + 2Nrr) + (n 1) + +) N..
N pm
1
follows that
Nur 2 ((2n + 2) + 4N, + 1 (2n + 3),
proving claim (iii).
A.1.1
PROOF LEMMA
For each local optimal solution n to formulation (15), the KKT condition states that there exist Lagrangian multipliers R and À R, such that
P&T = 0 Vk E [K],
(16a)
Pk 7 + n1++1 0
Vk [K],
(16b)
where no := and n == 1. For example, for the solution n with nt = 0 for all k E [K - 2].
k (Nnr 1)/2, and ni (Nnr + 1)/2, we pair it with the Lagrangian multipliers X
(Nnr + 1)/2 and
(N nr + +1)
if E [K -
1
if K
Wk E [K].
if -
Since (n*, X*) satisfies (164)-(16b) as well as (15b)-115c, n is local optimal solution to (15) with objective value (Nnr + 1)2/4. In what follows, we shall show that any other local optimal solution to (15) does not achieve strictly larger objective value, establishing the global optimality of n²
To this end, we multiply both sides of (16b) by ny and then sum them up for all & [K] to yield
K
K K
Hgng
n
n 1 + nK = 0
km
2
K
1
+ nK
(Nnrd + -nk),
2
14

---

Published as a conference paper at ICLR 2024
which follows from (16a) and (15b). This rewrites the quadratie program (15)
1 ] (Nard +n)
max n PER
2 s.t.
(16a)-(16b).
For any local optimal solution together with its corresponding Lagrangian multipliers (n, 2) sat- isfying (15b)-(15e) and (16b), we discuss the following three cases.
1. If μK 0, then (16b) with K implies that
= p- PK = 7k-11 + = = #k-1 + 1.
We denote m ;= k=1 nh with 0 - m S Nur. The following two sub-cases show that n cannot be better than
(a) If n 0, then À == 1 and nk Nur Exsl ne Nor m. It follows that the obejetive value of n is
(N 1)²
(N + +gk) N 5mh S Nar
that is, smaller than that of n.. This finishes the proof.
(b) If nx-1 > O, then μK- = 0 by (16a). It follows from (16b) with k = K that
nK-2+Rk --C---------- -K-1-- X
nK-1 -
But 1K + TKK = N - m by definition of m. This implies that
RK-1 = 2 (Nur - m nk-2 and nk 2
(N m 7x
Then, the objective value of n
(Nard + +nk)=
2 Nn +
2"" n
2
N + (nx-2+ +2 22-m)aar + nK- 2 m)
1 N + 22N + +(1 (1- (m2mm
< (N + 1)²,
where the first inequality because, with fixed m, setting nk-2 = m maximizes the objective value. Hence, the objective value of n is smaller than that of n², finishing the proof.
2. If pk and pK- 0, then nıK 0 by (16a) and (16b) implies that, with k = K-1,
0 = μK-1 - + nk-2 + ngK = -X + 72--2
that is, A nk-2. The following two sub-cases show that n cannot be better than n'.
(a) If 4K-2 > 0, then = n 0 by (16a). Hence, the objective value of n is (N + nx)/2 smaller than that of n²" This finishes the proof.
(b) pK-2 then (16b) with is K 2 implies that
0 = μK + nk 3 + カ + = пи + nk-1 + A
that is, A nk 7k But (15b) implies that
Nur
Nur 2 nk nk + MK-1 21 = d S
2
It follows that the objective value of n is
N
(Nur + 1)2
1 (N nx) <
4
4
that is, smaller than that of n². This finishes the proof.
15

---

Published as a conference paper at ICLR 2024
3. If Mx > 0 and μK-1 O, then nx R--1 0 by (16a). L [K 2] denote the largest index such that n > 0. Then, ,4, = by (16a). We discuss the following three sub-cases about j to show that cannot be better than n².
(a) j then in O and n2 0 by definition of j. It follows from (16b) with that
0 n2 + p1
0.
Hence, the objective value of n (NarA + nx) /2 = 0, smaller than that of n.. This finishes the proof.
(b) If j 2 2 and n > O, then Pj-1 = Mj = by definition j. It follows from (16b)
with & = j and le j that
Mj-1 - n3 =0 X = - +ij. M - - =0 =
À == == Ri+1+
But (15b) implies that
N
Nnr 2 nz
2
follows that the objective value m
J [NurA nk) S
N2
(N + 1)²
2
that is, smaller than that of n². This finishes the proof.
(c) j 2 2 and nj-1 = 0, then (16b) with k = implies that
= μ + カ5ー n5+1 =
= A ==
where the second equality follows from the definition of j. Hence, the objective value of n is (Nar + n2 O, smaller than that of n². This finishes the proof.
A.2 PROOF OF PROPOSITION
Consider functions / R2, g R" R, and their composite [g(x)). It can be shown that [[(())) is supermodular in x if one of the following is satisfied (see, e.g., Proposition 2.2.5(c) in Simchi-Levi al. (2005)):
1. f is increasing and convex, and g is increasing and supermodular;
2. f is decreasing and convex, and g is increasing and submodular.
Here, function g is called increasing if g(x) S g(x') whenever x - x', decreasing if -g is increasing, and submodular if -g is supermodular.
We proceed to prove by induction. For each A K, suppose that 2 is increasing and supermod- ular in x. Then, W - + by + Dạ is increasing and supermodular in x because Wa 2 0 and D& 2 O. It follows that 2ą is increasing and supermodular in because is increasing and convex. In particular, žх increasing and supermodular in x.
Finally, ò is supermodular in x by (6c) because WK-1- 2 0 and D Dx+12 supermodular (actually linear) in This finishes the proof.
16

---

Published as a conference paper at ICLR 2024
A.3 PROOF OF PROPOSITION 4.
Define x ;= 1 - for all i E [n]] and my == min [0.1)" 6(2). Then, any arbitrary o: (0, 10" 7 can be rewritten
6(x)
Σ
((2)
200 + (1 2x)
2E (
72
2500 + (1 - 2:)oo(] - (n - 1), of
Σ
(>)) max
=E(0,1)*
Σ
((z) o
Zi30is + - 2:)aa] (n 11 11
=E(0,1)
Σ (Ø(z) mo)o Σ =E(0,1)
2e04 + (1 =))ad] (n 11)
m +
where the second equality because all 24 and x in the product are binary, the third equality is by the definition of the activation function, and the last equality is because if and only
x. Then, the coefficients of all T, and x are nonnegative (because z, 0 and 2 0), and the coefficients of all activation functions are 6(2) ms 2 0. It follows that ((x) can be represented using an ISNN with an ar- hidden layer, which consists of 2" neurons and nonnegative coefficients only. This proves claim (i).
(E= [,,0, + (1 - =) (n-1))
chitecture as in Figure with K
For an ISNN as shown in Figure 2, ē(²) admits closed form expression as in (6). Since W1:+11 +1 and D2.K are nonnegative, the expression 2k,5, the ith entry of 2k, in ISNN can be rewritten as
Zhii
Wk.tf&-11 + bR. + Dr
Σ Wu.jo(Wa-11JF-22 +
+b b + DR- DR-15T) + ba.e + Du.
Σ WoWWEiWAL-JjR4-22 Wa.tbu-15 + Wk.ss Vx- + bu. + Dr
Σ @WW + 1/-1.5 + D2-1 + byr + Dr,
where WK,i denotes the entry of Wt row i and column j. W denotes the jth row W and
W W = Wk 2 O, bk-1=== Wk.sbba1.j Dz-1. Wk.ij Dk-15 2
That is, the same zk,i can be represented by an ISNN with Wa.y = for all i, ie., We is an all-one matrix. In light of this, the ISNN admits the following equivalent expression:
21 (W bb)
24. = o(1=-1 + by + Daā),, k == 2..,K K
0 1 Tzx +bx+1 + DK+1I,
17

---

Published as a conference paper at ICLR 2024
where 1 denotes an all-one vector all-one matrix. Hence, the total number of parameters Npm to be trained in the ISNN
K
Npm
Toi (2n +
+
Σ
The (2n + 1) (2n +
km
Σ nk +1
(2n + 1)
-NNr + 1) + 1),
where 24 E IRT for all [K] and so the number of neurons in the ISNN is Nur = Eker nu. Therefore, for fixed Nar, Npm is independent of W and K. This proves claim (ii).
To fit the N₂ sample-label pairs exactly, we need the number of parameters Npm of the ISNN to be at least large N,, i.e.
Npm- (Nur + +)(22n + 1) 2 N..
It follows that
N₂ 1,
Nur 2
2n + 1
proving claim (ii).
A.4 APPROXIMATION ERROR OF Ø
We evaluate the error of using the obtained from neural networks (GNN ISNN) to approximate the true (but unknown) value function c, defined as ||ø -ollo max lØ(x) Φ(x)). The following proposition provides an upper bound of the approximation error, as function of the lower-level problem and the neural network.
Proposition 6. Consider mixed-integer and linear lower-level problem, i.e,, g(y,x) == g and (xx) == ty E KT- 2., Gy Tx + k.), where R G E Raxm T E R8xx and h E 189 Then, holds that
- 2llco < (Cllalla + L)d,
where is a constant depending only on n and G. L is a constant depending only on the neural network coefficients {W K+1 and {Dx
K+1 and d is the Hausdorff distance between the sets of
sampled and unsampled points in X, i.e.,
d: max min | 21|2.
ZEX ZE
Pick any x' E X. First, using the representation (6) of the neural network, have || - zill2 & [Willl(x x'), and
Proof:
| <L1l2 S |[alle - - 2
|Dallla || x' 2'112
for all k == 2.... K. A mathematical induction on A: yields that
lō(x) - 3(0')1 I T|| x' I|2,
where L = -I K+1 [ + E
K+1 | K+1 [Wll)))
Second, since d(x) represents the optimal value of mixed-integer linear program parameterized by x, by Cook et al. (1986) we have
(x(x) ((x') s Cillella
+
where Cr, C2 constants that depend only and G. But since x, xf E [0,1)" - o(x') 0 when æ = x' and || 22112 2 1 when x o s'. It follows that
(x(x)) x(x')1 C lll1|2 2 || o'lla Callella ||x x'112
C ||او| ||ac x lll.
18

---

Published as a conference paper at ICLR 2024
where C C1 C2.
Third, pick a E urgmin.eol || - 2112. Then,
(x(x) [())1 = 10(x)) ((2) + (2) (=) {())
S (0(z) =0(2)++++++++ (0(2) - 0(x)
S C 119|22 || ---= L || 2||2
(Cl1g112 + D)|l= - =212,
where the first equality is because 2 E f2. Since æ is arbitrary, taking the maximum over all E X of both sides of the inequality finishes the proof.
B EFFICIENT CALCULATION FOR ISNN
Although Section 3 provides framework to solve general BPs through value function approxi- mation, there are total Nur neurons in the neural network, we need to introduce Nur auxiliary binary variables 6, Nur auxiliary continuous variables z', and 4N, big-M constraints (see (11)) to represent o(ž), which weakens its scalability.
replace (12c) with linear inequalities without introducing auxiliary binary variables or big-M M constraints as in (11). To this end, define set [2n] = {1, ..... 2n and the indicator set of binary vector E as S(ū) {1
Thanks to the supermodularity of ((t) promised by Proposition 3,
2n 2 = 1}. addition, we define function yp 2 E such that ((((()) ((5)
for all ā E (0, 1 Then, the following proposition recasts the (nonlinear) constraint g(v, 2 ((ū) as linear inequalities.
Proposition 7. (Adapted from Theorem 6 of Nemhauser & Wolsey (1981)). If (()) supermodular in I, then for all ā E {0, 1 g(-) 2 x(ã) if and only if
g(-) 2 ((S) Σ
p[[2m] [k], - 3x)
KES
(17)
((S,k)īK
AS S [2m]
\S
for) where p(S,k) : ((S) all S C [2m] and l e [2n]AS..
As compared to the reformulation (11), (17) does not introduce any new auxiliary variables. Yet, it involves an exponential number of linear inequalities which can nevertheless increase the formula- tion size. Fortunately, the inequalities (17) can be easily separated, that is, given fixed & and other variables. in polynomial time one certify that g(-) 2 o(i) find an violated inequality (17) with respect to some S. This is because a worst-case S on the right-hand side of (17), i.c.,
5 E argmax (((S) - Σ
(((2m) ((A), tx)
SC(2m]
KES
Σ
P(S, k)
admits a closed-form solution S.t (see Ljubié and Moreno (2018) and al. (2022))- Consequently, there is need to incorporate the exponentially many inequalities (17) up front, and we only need to incorporate the violated ones on-the-fly, e.g., in a branch-and-boun algorithm for solving the BP.
C INSTANCE GENERATION
Following the instance generation rules of (Tahernejad et al., 2020), all instances are generated in the following forms:
19

---

Published as a conference paper at ICLR 2024
min cxx + diy
s.t. A s b1
x E {0, 1
y E arg max d2
s.t. A2x B2y S b2
0 < S
y E RdTh (or y 2
where è is 1 vector, d1 and d2 X 1 vectors. The constraint number is set to be the same as the variable number of each level, and hence, A is a n x n matrix, b is a n x 1 vector, A2 is x n matrix, B2 is X matrix, and b2 is m vector. For instances with linear lower-level problem, y is continuous variables; for instances with a mixed-integer linear lower-level problem, y is integer variables. In our numerical experiments, we set m 20, y 1, and generate instances with n = 10, 20, 30, 40, 50, 60. The coefficients are randomly generated in a range given in the following table, where 8 200/ (m + n) and d1 d2.
Table I: Range of randomly generated coefficients.
COEFFICIENT
RANGE
COEFFICIENT
RANGE
C, di,,ds
50, 50]
A2
106, 106
A
20, 20
B2
-0, ò
b
30, 130
by
10, 110
D SUPPLEMENTARY RESULTS
D.1 OPTIMALITY CUT IN ILLUSTRATIVE EXAMPLE
profec
2
a5
1.5
(a) Exact.
Frocile riginn ofrcaud prokins Feyba prin
(yyleulty exper Faesibe pibli
15
15
65
45
155
85
1.5
12
(b) GNN.
(c) ISNN.
Figure 8: Optimality cut.
20

---

Published as a conference paper at ICLR 2024
D.2 DETAILS IN GENERATED INSTANCES
45
Nirraai Nivabo
35
2
20
GNN
ISNN
GNN
ISVN
(a) n1 10.
(b) 11 20.
40
-0
35
40
25
30
20
15
GNN
GNN
(c) 30.
(d) 40.
Maraios Nirriiss
Nicmiia
35
40
-5
25
15
(c)
50.
60.
Figure Distribution of objective difference in different instances with linear lower level.
1800
Sampling
Training
Solving
1500
1200
600
300
a=10
m=22
am330
n=40
n=50
n=60
a-10
n=
m=30
m=40
a=50
n=60
n=10
n=20
a=30
n=40
a>50
n=60
Niteration=1
Niteration=2
Niteration=3
(a) GNN.
Sampling Training Solving
1200
Thme
900
600
300
n=10
n=20
m=30
n=40
n=50
n=60
n=10
n=20
n=30
n=40
n=50
n=60
n=10
n=20
n=30
n=40
n=50
n=60
Niteration=1
Niteration22
Niteration=3
(b) ISNN.
Figure 10: Itemized time in different instances with linear lower level.
21

---

Published as a conference paper at ICLR 2024
Nicralia Steraia 22 Nimaise-))
Nicralia Nicaise-2 Skanika
14
0.
-5
GNN
ISNN
GNN
ISNN
10.
20.
Nivration 1 Moratia 21 Niwatce))
Nacsikon -1 Nieniier 2 Nicviiss
09
20
10
-20
GNN
ISNN
30.
40.
Nicraiia Nivraies 2 Nioviss
Nivration -1 Niwatise Niwation))
13
14
GNN
ISNN
GNN
ISNN
50.
60.
Figure 11: Distribution of objective difference different instances with mixed-integer linear lower level.
2200
Sampling Training Solving
1800
1600
1400
1200
800
600
400
200
n=10
n=20
n=30
n=40
n=50
n=60
n=10
m=20
n=30
n=4
n=50
m=00
n=10
n=20
n=30
n=40
n=50
n=60
Niteration=1
Niteration=2
Niteration=3
(a) GNN.
2000
Sampling
Training Solving
1800
1400
1200
1000
800
400
200
n=10
n=20
n=30
n=40
n=50
n=60
n=10
n=20
n=30
n=40
n=50
n=60
n-10
n-20
n=30
n=40
n-30
n=60
Niteration:2
(b) ISNN.
Figure 12: Itemized time in different instances with mixed-ine nixed-integer linear lower level.
22

---

Published as a conference paper at ICLR 2024
D.3 COMPARISON OF SAMPLING METHODS
We compare the efficiency of three sampling methods, including random sampling, Latin hypercube sampling (LHS), and the proposed enhanced sampling. We fixed sampling time of 30 minutes and adopt the generated MILP instances for numerical experiments. We report the number of ob- tained samples Figure 13. We note that when 10. the maximum possible number of samples is 1024, for which enumeration can be adopted to find all feasible z. Yet when n 2 20, enumeration becomes impractical and sampling is required.
From Figure 13, we see that the performance of random sampling is unstable and strongly depends on the lower-level problem. When the lower-level problem has a relatively large feasible region, such the instance 20 or n 60, it is easy for random sampling to find feasible for which sampling efficiency is high. Yet if the lower-level problem has relatively small feasible region, 30
such the instance and can not find enough samples for training. In addition, because we do not have probability infor- mation, LHS has similar performance random sampling most instances. However, through the proposed enhanced sampling, we guarantee to find a feasible æ each sampling, and from Figure 13, the number of obtained samples almost linearly decreasing with n. It is reasonable because larger n causes longer time consumption for each sampling and thus fewer samples under fixed sampling time. Therefore, avoid the unstable performance of naive random sampling other non-trivial sampling, we adopt the proposed enhanced sampling for all instances.
40, the sampling efficiency of random sampling gets quite low
30
Random
LHS
Proposed
Semple 10
25
20
15
10
20
40
60
Figure 13: Neural network architecture.
D.4 HIGHER DIMENSIONAL INSTANCES
We generate some larger instances with n up to 120. We consider a MILP lower level and conduct experiments with similar setup Section 4.2.1. Figure 14 shows the numerical results.
From Figure 14(a), the average objective difference tends to increase as n gets larger. When > 100, the average objective difference gets higher than 5%. The results are reasonable be- cause more samples are required to approximate (xx) when gets larger yet our sample number is fixed. From Figure 14(b), the computational time increases
gets larger, and when > 100, the computational time is longer than the set time limit of MiBS. Similar to Figure 12, sampling con- sumes the most computational time, for which more advanced sampling required in future works.
MIDDS
GNN
ISNN
27200
10
5
3600
0
1800
GNN
-SNNNN
100
(a) Average objective difference.
(b) Computational time.
Figure 14: Instances with a MILP lower level.
23


Reviews:
[{"review_id": "gctoahEluE", "reviewer": "Reviewer_skLM", "summary": "This paper employs neural networks to help solve the bilevel programs with binary tenders. Specifically, they adopt neural networks to approximate the optimal value of the lower-level problem as a function of the binary tenders. In order to train the neural networks, the authors also introduce an enhanced sampling method to generate high-quality samples. Lastly, numerical experiments were conducted to demonstrate the performance of these approaches.", "questions": "Major questions:\n1. In Figure 5 and Figure 6, the computational time of MiBS increases very fast with the increase of n. My question is, even though the MiBS solver might take a long time to reach optimality, have you considered to set a time limit, and compare the relative error of the best found solution given by the solver at time limit with your approaches?\n2. In your Conclusion section you mentioned: \"we demonstrated that the enhanced sampling helps reduce average relative error\". However, the whole point of your sampling method is simply to get enough data points for training the neural network. In order to claim that your proposed enhanced sampling can bring some extra benefits, you need to at least compare with other non-trivial sampling methods.\n3. In Conclusion section: \"The computational time of using ...... is significantly shorter than that of MiBS ......\" I admit that this is true, but your methods also do not reach optimality, for fair comparison you either need to have enough samples for exactly learning the value function (so that your approach can also produce a true optimal solution), or you need to compare the best feasible solution found by MiBS within a given time limit. \n\nMinor question:\n1. I suggest to give an exact definition for \"binary tender\". As far as I know, I do not think this is a well-known term in the community.\n2. Still about the definition of \"binary tender\". On Page 2 \"we assume that the entries of x appearing in the lower-level formulation are binary-valued\". However, in Appendix C, about the instance generation for the experiments, I notice that the binary constraint on x is enforced in the upper-level instead of the lower-level formulation.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "This paper gives a good attempt of incorporating ML (especially neural networks) methods to facilitate the solution of traditional mathematical optimization problems, which in my opinion is an area that deserves more attention from the community. Overall, this paper is clearly written, easy to understand, and the theoretical results in Section 3 are very neat. I also find it to be very impressive that I cannot find any typo throughout this paper and the propositions also appear to be of their own independent interest.", "weaknesses": "Main concern:\n1. The lack of ablation study, especially on the enhanced sampling part. For instance, why do you want to solve the quadratic programming problem (5) to get the samples? I understand that matrix Q is selected to be PSD is for the polynomial-solvability, but what is the main reason of solving the quadratic program in the first place? If we replaced this enhanced sampling with some other more naive sampling methods, how would it affect the experiment results?\n2. Limitation of the experiment setup and analysis: \n   (a) Instance dimension up to 60 is too small.\n   (b) The selection of the lower level problem is LP and MILP, which both have linear lower level objective function. At least you could have tried some simple nonlinear functions like quadratic function.\n   (c) The experimental results do not support the conclusion well. I list some of the points in the next Questions section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tE3PTzvlVB", "reviewer": "Reviewer_nqbY", "summary": "This work studies the Bilevel programs (BPs)  with discrete decision variables. A neural network is trained to approximate the optimal value of the lower-level problem, as a function of the binary tender. Then a single-level reformulation of the BP through a mixed-integer representation of the value function can be obtained. Moreover, an enhanced sampling method is proposed for high-dimensional BPs.", "questions": "I have some questions about the complexity, the comparison with the existing approximation-based methods in bilevel optimization, and the theoretical guarantee of the proposed method. Please see the Weakness.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. This work proposes an approximation-based method for Bilevel programs (BPs)  with discrete decision variables, which is interesting.\n\n2. an input supermodular neural network (ISNN) is proposed, which ensures a supermodular mapping from input to output.\n\n3. an enhanced sampling method is proposed for solving high-dimensional BPs.", "weaknesses": "1. The author should conduct some complexity analysis, such as time complexity [1], to show the effectiveness of the proposed method.\n\n2. This work employs neural networks to learn and approximate the value function $\\phi(x)$. However, training the neural networks is more computationally complex than directly approximating the lower-level optimization problem [2]. What is the advantage of the proposed method over the polyhedral approximation in [2]?\n\n3. Since one of the key contributions in this work is to employ neural networks to learn and approximate the value function $\\phi(x)$, I suggest the authors clearly discuss the existing approximation-based methods (i.e., which approximate the lower-level optimization problems and reformulate the bilevel optimization problems as single-level optimization problems, for instance, polyhedral approximation) in bilevel optimization since I can't find any discussion about the approximation-based methods in bilevel optimization.\n\n4. Training the neural networks to approximate the value function $\\phi(x)$ may introduce some variance which will lead to the solution of the resulting single-level problem far away from the original bilevel optimization problems. Can you provide a more theoretical guarantee for the proposed method?\n\n[1] A Gradient Method for Multilevel Optimization. NeurIPS 2021.\n\n[2] Asynchronous Distributed Bilevel Optimization, ICLR 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "x107YoWXbJ", "reviewer": "Reviewer_YfM6", "summary": "The authors explore the use of machine learning techniques to solve complex bilevel programs with binary variables. They introduce a sampling algorithm for obtaining high-quality sample data for train neural networks. Two neural networks ( general neural networks (GNN) and input supermodular neural network  (ISNN) are proposed to estimate the lower-level value function based on upper-level decisions. These estimates are used to enhance the optimization process by adding optimality cuts to a single-level problem.", "questions": "1) What is the stopping criterion for Algorithm 2 so that the generated point $(x,y)$ is guaranteed to be feasible to Problem 1? That is, $x,y$ is satisfied $x \\in X(y), y \\in Y(x), g(y,x) \\ge \\max_{z \\in Y(x)}  g(z,x)$ \n\n2) What is the definition of ``relative error” in Figures 4 and 7? How do you define an error when an optimal solution cannot be found, for example, for large-scale problems with n = 60?\n\n3) Why relative errors can be negative (as in Fig. 7)?\n\n4) Why we need to care of supermodularity of an ISNN? I notice that ISNN requires more neurons compared to GNN (i.e., the lower bound value for $N_{nr}$ in Eq. (8) vs Eq. (7))\n\n5) Do we have to solve $\\max_{y \\in Y(\\hat{x})} g(y, \\hat{x})$ to optimality in Algorithm 1?\n\n6) For feasibility purpose, why do we need the second term $h^T x$ in Eq. (5)?\n\n7) I don’t fully understand the notation ``GNN representability” in Prop 1 (ii). This statement “we measure the GNN representability by the number of parameters it trains” seems to be not mathematically sound. \n\n8) The architecture details of GNN of ISNN should be reported.\n\n9) The name of section 3.3 seems to be incomplete. \n\n10) I don’t see any part of Figures to support the claims “achieving relative error less than 5%,”, “the relative error of instances with n = 30, 50, 60 is larger than 15%.”", "limitations": "", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The authors propose a machine learning-based algorithm for solving bilevel programs with binary variables in the outer level problem. \n\n- The authors are able to provide some theoretical analysis for proposed neural networks with binary inputs. \n\n- The paper is easy to follow.", "weaknesses": "1) The numerical experiments lack sufficient evidence to demonstrate the advantages of the proposed algorithm. It appears that the algorithm can only achieve optimality for test instances with n = 10. The authors solely conduct comparisons with the baseline method MiBS based on computational times, without assessing solution quality, such as the optimality gap.\n\n2) The test instances are only from the authors’ synthetic generated ones. Some public datasets like MIBLP-XU and IBLP-FIS in Tahernejad et al., 2020 should be used.   \n\n3) This method should only be suitable when functions $f$ and $g$ are linear or MILP-representable because one needs to solve Eq (5) and (12) to optimality.\n\n4) The cost to solve a single instance of Problem (5) is very expensive because it is at least as hard as a mixed-integer quadratic program (MIQP).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)