ING TWO-STAGE STOCHASTIC INTEGER PROGRAMS
Year: 2022

Abstract:
Many practical combinatorial optimization problems under uncertainty can be modeled as stochastic integer programs (SIPs), which are extremely challenging to solve due to the high complexity. To solve two-stage SIPs efficiently, we propose a conditional variational autoencoder (CVAE) based method to learn scenario representation for a class of SIP instances. Specifically, we design a graph convolutional network based encoder to embed each scenario with the deterministic part of its instance (i.e. context) into a low-dimensional latent space, from which a decoder reconstructs the scenario from its latent representation conditioned on the context. Such a design effectively captures the dependencies of the scenarios on their corresponding instances. We apply the trained encoder to two tasks in typical SIP solving, i.e. scenario reduction and objective prediction. Experiments on two graph-based SIPs show that the learned representation significantly boosts the solving performance to attain high-quality solutions in short computational time, and generalizes fairly well to problems of larger sizes or with more scenarios.

heading: 1 INTRODUCTION
text: Stochastic integer programs (SIPs) are a class of combinatorial optimization problems (COPs) with uncertain parameters. They focus on finding a best decision for a practical task to optimize both the deterministic cost and the expectation of stochastic costs on a set of potential scenarios. Since the uncertain elements commonly exist in reality, SIPs are widely applied and studied in many fields such as transportation, inventory management, energy optimization, and so on (Jammeli et al., 2019; Dillon et al., 2017; Bruno et al., 2016; Wallace & Ziemba, 2005). However, optimally solving SIP is intractable especially when many scenarios are involved, due to the fact that: 1) most integer programs themselves are already NP-hard; 2) considerable (even continuous) scenarios significantly increase the computational complexity (Klein, 2021). In reality, solving moderate-sized problems may require prohibitively long time. Thus efficient and high-quality approximate solutions for SIPs are often practically pursued.
The advances in deep (reinforcement) learning has aroused extensive studies on solving optimization problems with neural networks. The related works have shown promising results in tackling some classic COPs, e.g., vehicle routing, job shop scheduling, combinatorial auction, etc. (Kool et al., 2019; Zhang et al., 2020; Lee et al., 2020; Kwon et al., 2020; Gasse et al., 2019; Chen & Tian, 2019; Ichter et al., 2018; Hottung et al., 2020; Wu et al., 2021). However, most of the existing deep learning based methods only focus on deterministic problems, which could hinder their practical applications in the uncertain situation. On the other hand, a few early attempts introduce deep learning to solve SIPs (Larsen et al., 2018; Abbasi et al., 2020; Joe & Lau, 2020; Bengio et al., 2020; Nair et al., 2018). However, they are generally designed for only specific problems and also lack the ability of generalizing across different problem sizes.
In this paper, we focus on graph-based SIPs, which is a general class of problems with broad applications. We propose to learn the deep latent representation of scenarios in SIPs, which can be
∗Equal contribution. †Corresponding author (zhiguangcao@outlook.com).
applied to various downstream tasks in solving SIPs, e.g., scenario reduction and objective prediction. Based on the basic two-stage form of SIPs, we leverage conditional variational autoencoder (CVAE) to learn latent and continuous representations of stochastic scenarios while considering the deterministic part of a SIP instance, which we term the context. Specifically, we employ a graph convolutional network (GCN) to learn embeddings of each deterministic integer program, which models the first-stage problem with a single scenario. The learned embeddings are mapped into continuous representations, and then decoded back into corresponding scenarios given the embedding of the context. This design ensures that the learned representations of scenarios are correctly linked to the corresponding instances.
The resulting latent space well represents scenarios in a class of SIP instances, and we demonstrate its power with two tasks commonly used in solving SIPs. First, we leverage it for scenario reduction by finding representative scenarios via an off-the-shelf clustering algorithm, which decreases the problem complexity and meanwhile attains high-quality approximate solutions. Second, we further extend the encoder with an additional sub-network to predict objective values for each scenario. Due to the generalization in the latent continuous space, we readily deploy the prediction task in semisupervised manner with a small number of target objective values. In doing so, the latent space is augmented with the predicted objective value and hence yields more representative scenarios, which further shrinks the approximation gap to the optimal solutions of the original SIPs.
In summary, this paper aims to learn scenario representation for solving graph-based two-stage SIPs. The main contributions are as follows:
• We firstly present a deep generative framework to learn representations for scenarios in graphbased two-stage SIPs. It considers dependencies of the scenarios on the deterministic context of each instance, and hence well differentiates scenarios in the latent space.
• We apply the learned representations to scenario reduction and objective prediction for SIPs. The representative scenarios derived from clustering of the representations can be used to efficiently compute high-quality solutions for SIPs. Moreover, we further shrink the approximation error by predicting objective values of scenarios via semi-supervised learning.
• We evaluate the proposed method on two SIP problems, i.e., the network design problem and facility location problem. Experimental results verify the superiority of the learned scenario representation for solving SIPs. Our method significantly outperforms existing deep learning based baselines. Notably, it also generalizes well to problems with larger sizes or more scenarios.

heading: 2 RELATED WORK
text: Learning based methods for SIPs. Recent works of deep learning for COPs mostly address deterministic problems, with scarce attention to uncertainty that commonly exist in practice. Though a few methods apply deep learning to improve approximate solutions for two-stage SIPs, they are mostly designed for some specific problems. Nair et al. (2018) train a bit-flipping policy with reinforcement learning to improve solutions iteratively, but it is limited to binary decision variables and assumes no constraints in the first-stage problem. Larsen et al. (2018) predict properties of optimal solutions to integer programs (i.e. the first-stage problem with respective single scenario) in a supervised manner. This method is only applicable when operational solutions are not needed. In contrast, Abbasi et al. (2020) predict values of decision variables for blood transshipment problem with fixed-sized instances. To reduce the complexity of SIPs, Bengio et al. (2020) aim to yield one scenario for the first-stage problem, so that the derived solution achieves a small objective value when evaluated with the original set of scenarios. This method assumes complete recourse, with any generated scenario being feasible for the second-stage problem. All the above methods attempt to reduce the expected cost of stochastic scenarios but their applications are limited by distinct assumptions or fixed problem sizes. In contrast, our method can be used for more general two-stage SIPs, and generalizes well to varying problem sizes or scenario quantities.
Scenario reduction for SIPs. In the SIP formulation, a set of scenarios is used to model stochastic events on top of the first-stage problem. It introduces extra complexity which may exponentially increase with its cardinality, owing to the combinatorial search space (Dyer & Stougie, 2006). Since the direct optimization with massive scenarios is intractable, the scenario reduction is often conducted to find a subset of scenarios that well replaces the original set with a small approximation
error. Existing reduction methods mainly branches into distribution-oriented and problem-oriented paradigms. The former ones pursue a group of scenarios that closely estimate the probability distribution under the scenario space (Henrion & Römisch, 2017). The distance between the estimated and actual probability are often measured by Wasserstein distance (Rujeerapaiboon et al., 2018), polyhedral discrepancy (Henrion et al., 2008) or difference of moments (Høyland et al., 2003). These methods neglect the influence of the context on the objective function, e.g., similar scenario distributions may cause different expected costs in disparate SIP instances. Thus, there is a recent trend towards the study of problem-oriented methods (Henrion & Römisch, 2018; Fairbrother et al., 2019; Keutchayan et al., 2020; 2021), which attempts to involve problem-specific properties in scenario reduction. For example, Keutchayan et al. (2021) propose to find K scenarios with their objective values separately approximating the expected objective values in the K scenario subsets. Most the above methods are mainly for theoretical proof on stochastic programs (without the restriction of integers). In contrast, we firstly propose a deep learning method for SIPs to learn scenario representation, which can be used for both scenario reduction and objective prediction.

heading: 3 PRELIMINARIES

heading: 3.1 TWO-STAGE STOCHASTIC INTEGER PROGRAMS
text: As an optimization problem, SIP is usually characterized by both the uncertain parameters that potentially follow certain probability distributions and the discrete solution space due to the integer restrictions. It is commonly described by a two-stage formulation as below:
min x
µ⊤x+ Eω[Q(x, ω)] (1)
s.t. Ax ≤ b, x ∈ Rn1−p1 × Zp1 , (2)
where Q(x, ω) := miny{q⊤ω y|Wωy ≤ hω − Tωx; y ∈ Rn2−p2 × Zp2}. In particular, Eq. (1) and (2) prescribe the first-stage problem with x ∈ Rn1 being the decision variable, where n1 ≥ p1; Q(x, ω) prescribes the second-stage problem, with y ∈ Rn2 being the decision variable, where n2 ≥ p2. We refer to the group of static parameters µ ∈ Rn1 , A ∈ Rm1×n1 and b ∈ Rm1 in the first stage as the context, and assume that the group of uncertain parameters qω ∈ Rn2 , Wω ∈ Rm2×n2 , Tω ∈ Rm2×n1 and hω ∈ Rm2 follows a distribution P . Here we focus on SIPs defined on graphs, which are a family of problems with practical applications in many domains such as networks, transportation and scheduling (Rahmaniani et al., 2018; An & Lo, 2016).
The primary method to solve the above SIP is sample average approximation (Kleywegt et al., 2002), which converts it into a mixed integer program (MIP) by Monte Carlo simulation and in turn optimizes the following objective:
O(x) := min x
µ⊤x+ 1
N ∑N i=1 Q(x, ωi), (3)
where {ωi}Ni=1 is a set of scenarios from P , i.e., an independently and identically distributed (i.i.d.) random sample of N realizations of uncertain parameters. Typically, a large scenario set is often required to make the distribution of scenarios P̃ well fit P , which may cause intractable MIP for a given solver. In this paper, we aim to find a small number of informative representatives in {ωi}Ni=1 by learning representations of scenarios, which is supposed to considerably boost the computation efficiency with tolerable approximation errors.

heading: 3.2 CONDITIONAL VARIATIONAL AUTOENCODER
text: As an unsupervised generative model, CVAE is developed on top of VAE, and it further controls the data generation process conditioned on additional random variables (Sohn et al., 2015). These conditional variables could be either class labels or certain properties of the data with specific distributions, which are engaged in the input to both the encoder and decoder. Optimized with the stochastic gradient variational bayes (SGVB) framework (Kingma & Welling, 2014), the evidence lower bound objective (ELBO) of CVAE on the marginal likelihood for the input data is defined as:
log pθ(X, c) ≥ ELBO(X, c) = Eqϕ(z|X,c)[log pθ(X|z, c)]− KL[qϕ(z|X, c)∥p(z|c)], (4)
where X , c, z denote input, latent and conditional variables, respectively; both the encoder qϕ and decoder pθ are typically parameterized by deep neural networks. Rather than simple distributions of conditional variables in most existing works, in this paper we exploit CVAE to learn scenario representation in SIPs, with conditional variables derived from continuous parameters in the context. For more details of CVAE, we refer readers to (Kingma et al., 2014; Feng et al., 2021).

heading: 4 METHODOLOGY
text: Given a class of SIP instances {Xm}Mm=1 with parameters drawn from a distribution D, we regard each instance as a 2-tuple Xm = (Dm, {ωim}Ni=1), where Dm denotes the context (i.e. the group of static parameters) in the m-th instance, and ωim denotes the i-th scenario (i.e. the i-th realization of uncertain parameters) in the m-th instance. We aim to learn latent representations (variables) {zim}Ni=1 of scenarios {ωim}Ni=1 under the context Dm in each instance.

heading: 4.1 CVAE FOR SCENARIO REPRESENTATION
text: Our CVAE comprises an encoder for the inference process and a decoder for the generation process, which are parameterized by deep neural networks qϕ and pθ, respectively. In the generation process, the decoder approximates the posterior distribution of uncertain parameters in scenarios, given the latent and conditional variables, such that:
pθ(ω i m|zim, cm) = f(ωim; zim, cm, θ); qϕ(cm|Dm) = h(cm;Dm, ϕ), (5)
where f and h denote likelihood functions (e.g. Gaussian or multinomial distributions) that are instantiated by pθ and qϕ, respectively. Note that we do not take the raw context Dm as conditional variables for the decoder since they are continuous and high-dimensional. It may not only intensify the computational complexity but also overwhelm the latent variable which is typically low-dimensional. In contrast, we leverage the graph neural network (GNN) in the encoder to derive a low-dimensional conditional variable cm from Dm. Further by preserving the marginal independence between latent and conditional variables, the decoder will more effectively govern the generation of ωim according to the context Dm (Kingma et al., 2014).
In the inference process, the encoder approximates the posterior distribution of latent representations (variables) of scenarios given (ωim, Dm), which is expressed as:
qϕ(z i m|ωim, Dm) = N (zim|µϕ(ωim, Dm), σ2ϕ(ωim, Dm)), (6)
where we assume a Gaussian distribution of scenario representations with µϕ and σϕ being parameterized by neural networks, respectively. As aforementioned, we exploit a GNN in the encoder to embed the continuous and high-dimensional ωim and Dm, which is also used by the decoder to embed the sole Dm. This parameter sharing design contributes to a fast learning of embeddings for deriving the latent and conditional variables, i.e., zim and cm. To learn disentangled representations for these two types of variables, we force the conditional independence as:
qϕ(z i m, cm|ωim, Dm) = qϕ(zim|ωim, Dm)qϕ(cm|Dm). (7)
In summary, our CVAE hinges on two important steps, i.e., 1) we first use the encoder to attain scenario representations with respect to the context in each SIP instance; 2) we then use the decoder to reconstruct the scenarios from the latent space according to the context.
Semi-supervised CVAE. With the inferred scenario representations, we can directly conduct various downstream tasks for solving SIPs, such as scenario reduction using an off-the-shelf clustering algorithm. However, it may ignore the relationship between scenarios and objective values, which is nontrivial since even similar scenarios would deliver different solutions. Therefore, we also extend CVAE to predict the objective function in a semi-supervised manner. Given that the scenarios are embedded into a continuous space by the encoder, the prediction is supposed to be well generalized (Gómez-Bombarelli et al., 2018).
Most CVAE models for semi-supervised learning predict discrete targets and regard the ground truth, e.g., class labels, as conditional variables. In contrast, we predict a continuous property of the data (i.e. the objective value) and directly approximate it through the encoder. To this end, we use
CPLEX solver to collect optimal objective values {Y im}Ni=1 of the resulting MIP problems, which are defined by respective scenarios with the context, for only a small volume (1% ·M ) of instances in {Xm}Mm=1. We refer to the underlying joint distribution of (ωim, Y im, Dm) as DY . Then we infer the objective values with a sub-network, which processes the GNN embedding of (ωim, Dm) in the encoder, such that the objective function σ is parameterized as:
rψ(Y i m|him) = σ(Y im;him, ψ); qϕ(him|ωim, Dm) = g(him;ωim, Dm, ϕ), (8)
where him denotes the GNN embedding that is derived by qϕ.

heading: 4.2 NEURAL NETWORKS
text: A major challenge to parameterize CVAE for SIPs is how to derive an effective embedding of the context, since it is high-dimensional and continuous. Moreover, it is also involved in the input to both the encoder and decoder, and significantly influences the eventual scenario representation. A direct use of multi-layer perceptron (MLP) to process static parameters in Dm may fail to leverage the problem structure and cannot generalize across instances of different sizes. To overcome these limits, we exploit a GNN (as aforementioned) to embed the context Dm for attaining conditional variables in the decoder. Meanwhile, we apply the same GNN to embed each scenario with its context (ωim, Dm) for attaining latent representations in the encoder. To this end, we first describe the context with scenarios on a graph.
SIP graphs. We define a complete graph G = (V, E), where V = {v1, · · · , vn} denotes nodes with features V ∈ Rn×dv ; E = {ejk|vj , vk ∈ V} denotes edges with features E ∈ Rn×n×de . We formally represent the features on the j-th node as vjωi = [s j ωi ; d
j ] ([;] means concatenation). In specific, sjωi ∈ ω i m denotes a realization of uncertain parameters on the j-th node and d
j ∈ Dm denotes static parameters on the j-th node. Similarly, the features on edge ejk is represented as ejkωi = [s jk ωi ; d jk], where sjkωi and d jk denote the parameters from the i-th scenario and the context, respectively. Different from existing works that only consider varying parameters on nodes, e.g., Nazari et al. (2018), our graph representation is more generic and could be applied to a broader class of SIPs. Particularly, two examples of the graph representation on the network design problem and facility location problem are introduced in Section 4.4.
On top of the SIP graphs, we can exploit GNN to derive embeddings for either the context Dm or varying scenarios with the context, i.e., (ωim, Dm), i = {1, . . . , N}. In our method, we adopt the graph convolutional network (GCN), which is a commonly used GNN variant (Kipf & Welling, 2017; Yao et al., 2019). The architecture of our GCN is similar to the one in Joshi et al. (2019), which was leveraged to compute an adjacency matrix for solving the travelling salesman problem. Specifically, we extend it for graph-based SIP problems to attain graph embeddings for Dm and (ωim, Dm), denoted as hm and h i m, respectively. Detailed structure of the GCN is described in Appendix A. Next, we will revisit our CVAE given the graph embeddings above.
Encoder. The encoder processes the graph embedding him with two separate linear projections to compute 2-dimensional vectors of means and standard deviations for latent variables, respectively, i.e., µϕ(ωim, Dm) and σϕ(ω i m, Dm) in Eq. (6). Also, we linearly project the graph embedding hm into a 2-dimensional vector, i.e., cm in Eq. (5).
Decoder. We concatenate cm with the latent variables zim and pass them through an MLP to reconstruct the scenario. The MLP is structured by two hidden layers with ReLU activation functions, which are of 128 and 256 dimensions, respectively.
Semi-supervised learning. To echo the semi-supervised CVAE in Section 4.1, we leverage another MLP as a sub-network to process him. It comprises a 512-dimensional hidden layer with a ReLU activation function, and outputs a single value to estimate the objective value Y im.
The basic structure of neural networks is illustrated in Figure 1, where ω̂im and Ŷ i m are the reconstruction and predicted objective value for the i-th scenario, respectively. Note that we pad uncertain parameters ωim with 0 to embed the sole context (in the lower grey square). In this way, the GCN could be shared to attain both graph embeddings hm and him.
4.3 TRAINING & INFERENCE Algorithm 1 Training procedure Input: encoder pθ , decoder qϕ, num. of epochs E , batch size B . 1: Initialize parameters θ(0) and ϕ(0) 2: for epoch = 1, 2, . . . , E do 3: Shuffle the data inD andDY ; 4: for t = 1, 2, . . . , ⌊ |D|B ⌋ do 5: Retrieve a batch of instances fromD; 6: Get a random scenario in each instance; 7: Compute ELBO in Eq. (9); 8: θ, ϕ← Adam(θ, ϕ,∇L(θ, ϕ)); 9: if objective prediction is trained then 10: Retrieve a batch of instances fromDY ; 11: Compute MSE of objectives in Eq. (9); 12: ϕ, ψ ← Adam(ϕ, ψ,∇L(ϕ, ψ)); 13: end if 14: end for 15: end for The training objective is to minimize three loss functions. The first one is defined by the mean square error (MSE) between the reconstructed scenarios from decoder and the original ones. It is essentially employed in CVAE to maximize the marginal likelihood of the data. The second loss is defined by Kullback–Leibler divergence from the latent variables to the prior Gaussian distribution p(z) = N (z|0, I), which is used for regularization. The third loss is defined by the MSE between the predicted objective values and their ground truth. Formally, the objective for training our CVAE is expressed as:
L(θ, ϕ, ψ) = ED[−ELBO(Xm)]− α · EDY [log rψ(Y im|him)], (9)
where ELBO(Xm) = Eqϕ(zim|ωim,Dm)[log pθ(ω i m|zim, cm)]−β ·KL[qϕ(zim|ωim, Dm)∥p(z)] involves the first MSE and KL divergence; α and β are hyperparameters to balance the three losses. Following the typical training paradigm for VAE, we jointly optimize θ and ϕ via reparameterization trick and Monte Carlo approximation (Kingma & Welling, 2014). The training procedure is summarized in Algorithm 1 where we use the Adam optimizer to update the parameters (Kingma & Ba, 2015).
As shown in Algorithm 1, we randomly pick one scenario for each instance so that we can process it with its context in a batch. The resulting generalization across both the contexts and scenarios of different instances will allow for a fast training of neural networks. However, in the inference we employ the trained encoder to attain the latent representations, by processing all scenarios from a set of instances in parallel. It is realized by duplicating the context in each instance for all its scenarios. In doing so, the computational efficiency will be considerably improved for solving the instance with a large cardinality of scenarios, and thus a large amount of instances. For each instance, the derived scenario representations from the encoder are then used for clustering to attain a few centers that correspond to representatives among original scenarios, so we end with reduced scenarios.
4.4 APPLICATIONS
We apply our method to learn scenario representation in two typical SIP problems, i.e., the network design problem (NDP) and facility location problem (FLP) (Keutchayan et al., 2021). We represent them as the SIP graphs, whose features on nodes and edges including both static and uncertain parameters are listed in Table 1. For FLP, the parameters tokenized by F. and C. mean that they are peculiar to facility and customer nodes, respectively. More details about the two problems and their parameters can be found in Appendix B and C.

heading: 5 EXPERIMENTS
text: We evaluate our method on NDP and FLP to demonstrate its effectiveness. On the one hand, we train the CVAE to only learn representations for scenarios (Line 5∼8 in Algorithm 1), and the trained encoder can attain representations in each instance, which are directly used for scenario reduction by finding representatives via an off-the-shelf clustering algorithm. On the other hand, we train the semi-supervised CVAE to simultaneously learn scenario representations and predict objective values (Line 5∼13). After the training, the neural network can attain representations and objective values for each scenario and we concatenate them for clustering to find more informative scenarios as representatives. We refer to the above two paradigms as CVAE-SIP and CVAE-SIPA, respectively. In the following experiments, we compare them with other learning based methods, and assess their generalization performance across problems with varying sizes and scenario cardinalities.

heading: 5.1 SETUP
text: Instance generation. Our instances are generated following the way in (Keutchayan et al., 2021). Specifically, we generate NDP instances with 14 nodes including 2 source nodes, 2 target nodes and 10 transition nodes. For FLP, we designate 10 facility nodes and 20 customer nodes. We generate 200 scenarios for each instance in both problems. To evaluate the generalization performance, we generate larger instances by doubling and quadrupling the number of transition nodes for NDP, respectively. We do the same to both facility and customer nodes for FLP. Besides, we also double and quadruple the scenario cardinality in both problems. More details are provided in Appendix D.
Training. We collect 12, 800 instances (of normal size) for NDP and FLP, respectively, and normalize each type of parameters into [0, 1] across instances. Regarding CVAE-SIP, we train our model without the sub-network (used for the prediction task) by solely optimizing the ELBO loss in Eq. (9). Regarding CVAE-SIP, we train neural networks for learning representation and objective prediction simultaneously, as shown in Algorithm 1. We set hyperparameters β = 0.005 and α = 100 in Eq. (9). For both problems, CVAE-SIP and CVAE-SIPA are trained with 100 and 400 epochs, respectively. The batch size and learning rate are set to 128 and 10−4. We run the training on a server using a single GeForce RTX 2080 Ti GPU.
Testing. Pertaining to CVAE-SIP, we employ theK-medoids algorithm (Park & Jun, 2009) to obtain representative scenarios based on the latent representations. We keep its default settings and only change the number of clusters (i.e. K). Given the attained K representatives, we can derive the approximation solution by reformulating Eq. (3) as:
Õ(x) := min x µ⊤x+ Nk N ∑K k=1 Q(x, ωk), (10)
where ωk means the representative from the k-th cluster and the number of intra-cluster population equals to Nk. Accordingly, the approximation error is defined as the gap1 between objective values of the approximate solution derived from Eq. (10) and the optimal solution from Eq. (3). Pertaining to CVAE-SIPA, we infer both the latent representations and objective values of scenarios, which are concatenated as augmented representations to attain the representatives. To reduce the influence by the potential concurrent processes, we conduct the evaluation on a different desktop computer, which uses a single Nvidia Quadro P4000 GPU with a Xeon W-2133 CPU@3.60 GHz.

heading: 5.2 COMPARISON ANALYSIS
text: We compare with four baselines2, i.e., 1) CPLEX IBM (2017), with default settings to compute optimal solutions. It is also used to solve the approximation problem with representative scenarios. 2)K-medoids, the same clustering algorithm as used in our method. However, we apply it to directly cluster the original scenarios rather than the learned low-dimensional representations. 3) ScenarioM, a supervised method with hand-crafted features to predict a scenario for solving a variant of FLP (Bengio et al., 2020). 4) Solution-M3, a supervised method with raw features to predict the first-stage decision variables for solving blood transshipment problem (Abbasi et al., 2020).
1In specific, we denote the optimal solutions in Eq. (10) and Eq. (3) by x̃ and x∗ respectively. The gap is defined as |O(x̃)−O(x∗)|/(|O(x∗)|+ ε), where we assume minimization problems and ε = 10−10.
2A short comparison with methods in operations research is also given in Appendix I. 3The training details of Scenario-M and Solution-M are summarized in Appendix E.
We evaluate all methods on 200 instances of normal size for NDP and FLP, respectively. We set K = 5, 10, 20 in K-medoids for both problems so that we could access the performance against the number of representatives. Note that the computation by CPLEX, Scenario-M and Solution-M is independent of K. The results are displayed in Table 2, where we record the average of objective values, approximation errors and computational time over all instances, respectively. As shown, both CVAE-SIP and CVAE-SIPA consistently outperform K-medoids in all cases. Though K-medoids can also improve the solution quality as K grows, our methods efficiently attain solutions with much smaller errors. It implies that we achieve more informative representatives with the learned low-dimensional representations, which help deliver high-quality solutions. On the other hand, our methods significantly outperform Scenario-M and Solution-M for the three respective values of K. This suggests that the representatives found by our methods are able to approximate the original problem more accurately than the approximation with the single scenario or solution. Furthermore, we observe that CVAE-SIPA is generally superior to CVAE-SIP, which indicates that the predicted objective values could effectively enrich the scenario representations for achieving more informative representatives and thus better solutions. Last but not least, all the learning based methods consume fairly short time compared to the sole CPLEX.

heading: 5.3 GENERALIZATION ACROSS PROBLEM SIZES
text: A desirable generalization to larger problems is necessary for learning based methods. To verify such ability of our methods, we apply the previously trained networks to directly infer another 50 larger instances (of two levels) for both NDP and FLP, respectively. The instances are generated following the operation of doubling and quadrupling described in Section 5.1. The baselines are also directly applied to these larger instances except Solution-M since its output dimension (i.e. number of decision variables) is fixed only for the problem of normal size. Similarly, Scenario-M is not applicable to FLP where the scenario dimension (i.e. number of customer nodes) changes with problem sizes. For CPLEX, we set a 1h time limit on FLP (120; 200). In Table 3, the upper and lower half contains the results on doubling and quadrupling instances, respectively. As shown, CVAE-SIP is superior to K-medoids and Scenario-M in most cases except NDP (44; 200) with K=5. On the other hand, CVAE-SIPA generally achieves the smallest errors but performs inferior to CVAE-SIP on NDP (44; 200) with K=10 and FLP (120; 200) with K=5. This might stem from the degenerated accuracy of objective prediction, which induces solutions with relatively large penalties on some instances.
Notably, the computational time by CPLEX rises drastically on larger problems, in comparison with those in Table 2. Such low efficiency might not be acceptable in practice, especially where real-time decisions are required for handling a number of problem instances simultaneously. In contrast, our methods deliver solutions with small penalties in short runtime. In Appendix F, we draw curves of the average error against the number of representative scenarios used in FLP problems.

heading: 5.4 GENERALIZATION ACROSS SCENARIO CARDINALITIES
text: From a practical view, the cardinality of scenarios in SIP problems might be changed due to certain reasons, e.g., new observations are added or outdated ones are removed (Issac & Campbell, 2017). Therefore, a learning based method for SIPs should be able to handle this varying dimension via generalization. To this end, we evaluate our methods and baselines on 50 instances of NDP and FLP with doubling and quadrupling scenarios, respectively. The results in Table 4 show that our methods can still achieve the smallest errors across most settings, except NDP (14; 400) with K=5, where Solution-M performs better. It implies that our methods are less sensitive to the scenario cardinality. Since our methods process all scenarios in parallel, they solve instances with scenarios up to 800 in about 2 seconds, which is significantly faster than CPLEX. Additionally, we also test with scenarios generated from varying distributions and different dependencies on the context. The details are provided in Appendix G.

heading: 5.5 OBJECTIVE PREDICTION
text: We further demonstrate the efficacy of the semi-supervised learning for objective prediction. The trained model is directly adopted to infer objective values of scenarios in NDP and FLP instances. The results are gathered in Figure 3 and 4, which show our prediction is fairly accurate for both problems even when generalizing to larger problem sizes or scenario cardinalities. It well verifies the effectiveness of our semi-supervised design, which only needs to solve a small volume (1%) of training instances. In addition, the sample efficiency is also enhanced through the generalization over the latent space. More details can be found in Appendix H.

heading: 6 CONCLUSIONS AND FUTURE WORK
text: In this paper, we present a deep learning method to solve graph-based two-stage SIPs. We leverage CVAE to learn latent representations of stochastic scenarios conditioned on the context of a SIP instance. With the learned representations, we conduct two downstream tasks for solving SIP, i.e., scenario reduction and objective prediction. The evaluation on two classic SIP problems has shown that our method is able to find less yet informative representatives for improving both solution quality and computation efficiency. Additional studies also verify that our method can be well generalized to larger problem sizes, larger cardinalities of scenarios and different distributions of scenarios. In the future, we will attempt to automate the decision on the optimal number of representatives for scenario reduction, and extend our method towards solving general two-stage SIPs.

heading: ACKNOWLEDGEMENTS
text: This research was conducted at Singtel Cognitive and Artificial Intelligence Lab for Enterprises (SCALE@NTU), which is a collaboration between Singapore Telecommunications Limited (Singtel) and Nanyang Technological University (NTU) that is supported by A*STAR under its Industry Alignment Fund (LOA Award number: I1701E0013). Wen Song is supported by the National Natural Science Foundation of China under Grant 62102228, and the Shandong Provincial Natural Science Foundation under Grant ZR2021QF063.

heading: A GRAPH CONVOLUTIONAL NETWORK
text: Given a complete SIP graph as defined in Section 4.2, we first process the features of nodes and edges by separate linear transformations as below:
V0 = VW0v + bv; E 0 = EW0e + be, (11)
where V ∈ Rn×dv and E ∈ Rn×n×de are the raw features (i.e. parameters in Table 1) on nodes and edges, respectively; W0v ∈ Rdv×128 and W0e ∈ Rde×128 are trainable weights with corresponding bias vectors bv ∈ R1×128 and be ∈ R1×128. Note that we also use different linear transformations for heterogeneous nodes or edges. For example, the location nodes of facilities and customers are processed separately.
Next, we advance the node and edge embeddings by the graph convolution layer for L iterations, and each iteration has its own parameters, which are expressed as:
vl+1j = v l j + PReLU
( BN ( vljW l 1 + ejk∈V∑ k ηljk ⊙ vljWl2 )) , ηljk = σ(eljk)∑ejk∈V k σ(eljk) + ϵ
el+1jk = e l jk + PReLU
( BN ( eljkW l 3 + v l jW l 4 + v l kW l 5 )) , l = 0, . . . , L− 1
(12)
where Wl1, . . . ,W l 5 ∈ R128×128; ⊙ means the element-wise product; σ is the sigmoid activation function; ϵ is a small value to avoid numerical error. Initially, the node embedding v0j is the j-th row of matrix V0 and edge embedding e0jk is the element at the j-th row and k-th column of matrix E0. The above graph layer is similar to the one in Joshi et al. (2019) except that we use parametric rectified linear unit with its learnable parameter being 0.25.
After L (L=10) iterations of the above convolution, the node and edge embeddings are represented as {vLj }nj=1 and {eLjk}nj,k=1, respectively. We attain the graph embedding by the mean pooling over all the node embeddings. We neglect edge embeddings considering that their information has already been involved in node embeddings, after multiple convolutions.
Parameter sharing The graph embedding network based on Eq. (11) and (12) is able to process the sole context Dm with its static parameters on the graph. Also it can process the context with any scenario in the instance, by featurizing the graph with both static parameters in Dm and a scenario ωim. To foster the learning efficiency, we use a single network to process these two featurizations. It is achieved by padding the uncertain parameters with zeros when we only process the context. We represent the resulting graph embeddings by hm and him, which are further used to produce conditional and latent variables and meanwhile predict objective values, as stated in Section 4.2.

heading: B NETWORK DESIGN PROBLEM
text: Given a complete directed graph G = (V, E), each edge ejk ∈ E points from vj to vk where vj , vk ∈ V . We denote the set of source and target nodes as S and T respectively. The features on an edge include its capacity qjk, opening cost ojk and transportation cost hjkc for one unit of commodity c ∈ C. Scenarios The uncertainty exists in the demands on each node. In specific, the parameter divc denotes the demand of commodity c on node v, in the i-th scenario. A penalty fs is given to a source node if its commodities are not fully transported in a scenario due to the capacity limit of opened edges.
Decision variables The decision variable in the first stage is denoted by xjk. xjk = 1 means that the edge ejk is open to transport commodities, otherwise xjk = 0. The decision variable in the second stage is denoted by yijkc, which means the units of commodity c transported along the edge ejk in the i-th scenario.
The objective of NDP is to decide the edges to be opened, which minimizes the total opening cost and expected transportation cost (including the penalty). The SIP formulation is expressed as:
min ∑ ejk∈E ojkxjk + 1 N N∑ i=1 (∑ c∈C ∑ ejk∈E tjkch i jkc + ∑ s∈S fszsc )
(13)
s.t. ejk∈E∑ vj=v yijkc − ejk∈E∑ vk=v yijkc = d i vc, ∀(v, c, i) ∈ V × C × {1, . . . , N} (14)
ejk∈E∑ vj=t yijkc = ejk∈E∑ vk=s yijkc = 0, ∀(s, t, c, i) ∈ S × T × C × {1, . . . , N} (15) disc ≤ ejk∈E∑ vj=s
yijkc +Mzsc, ∀(s, c, i) ∈ S × C × {1, . . . , N} (16)∑ c∈C yijkc ≤ qjkxjk, ∀(vj , vk, i) ∈ V × V × {1, . . . , N} (17)
zsc, xjk ∈ {0, 1}; yijkc ∈ [0,∞). (18) where M is a large value which ensures that the source node will not be penalized if all its commodities are transported.

heading: C FACILITY LOCATION PROBLEM
text: The facility location problem comprises F and C locations for facilities and customers, respectively. Each facility location is featured by an opening cost of , f ∈ {1, . . . , F} and a facility has q units of resources that can be provided to customers. We assume that the number of facilities to be opened is not larger than v. For each customer c ∈ {1, . . . , C}, it consumes qcf units of resources if it is served by the facility at the location f , and the induced service cost is denoted by scf .
Scenarios The uncertainty exists in the presence of customers. In specific, the parameter hic refers to the presence of the customer c in the i-th scenario, i.e., hic = 1 means that the customer c is present to be served in scenario i, otherwise hic = 0. A penalty bf is given to the facility f in a scenario if one additional unit of resource is needed for it to meet the demand.
Decision variables The decision variable in the first stage is denoted by xf . xf = 1 means that a facility is open at location f to provide services, otherwise xf = 0. The second-stage decision variables include yicf and z i f . In specific, y i cf = 1 means that the customer c is served by facility at location f in the i-th scenario, otherwise yicf = 0. z i f means the quantity of additional resources required at facility location f to meet the total demand of its customers, in the i-th scenario.
The objective is to decide the locations to build the facilities, so as to minimize the total opening cost and the expected service cost including the potential penalty. The two-stage SIP formulation of the facility location problem is expressed as:
min ∑ f∈F ofxf + 1 N N∑ i=1 (∑ c∈C ∑ f∈F scfy i cf + ∑ f∈F bfz i f ) (19)
s.t. ∑ f∈F
xf ≤ v, (20)∑ c∈C qcfy i cf ≤ qxf + zif , ∀(f, i) ∈ F × {1, . . . , N} (21) zif ≤Mxf , ∀(f, i) ∈ F × {1, . . . , N} (22)∑ f∈F yicf = h i c, ∀(c, i) ∈ C × {1, . . . , N} (23)
xf , y i cf ∈ {0, 1}; zif ∈ [0,∞), (24)
where M is a large value that constrains zif to 0 if xf = 0, so as to ensure that the facility f will not be penalized if it is not open.

heading: D DETAILS OF INSTANCE GENERATION
text: We generate NDP instances with 14 nodes, in which the commodities on two source nodes are transported to another two target nodes to meet their demands. The remaining nodes are only for transition without demands. We set two types of commodities and their quantities are uniformly sampled from [5, 15] in each scenario. For each edge, we uniformly sample the opening cost, shipping cost and capacity from [3, 11], [5, 11] and [10, 41], respectively. A penalty 1000 will be triggered for each source node if its commodities cannot be fully transported due to the capacity limit of opened edges. For FLP, we consider 10 facility nodes and 20 customer nodes, with their coordinates uniformly sampled from the unit square [0, 1]× [0, 1]. The opening cost and capacity for each facility are uniformly sampled from [40, 81] and [30, 60]. We generate binary values through Bernoulli sampling for customer nodes to indicate their presence in a scenario, with probabilities uniformly sampled from [0.4, 0.6] for each node. A penalty 1000 is set for each facility node if it cannot meet one unit of demands. We generate 200 scenarios for each instance of both problems. To test the generalization performance, we generate larger instances by doubling and quadrupling the number of transition nodes for NDP, respectively. We do the same to both facility and customer nodes for FLP. Regarding the scenario, we also double and quadruple its cardinality in both problems.

heading: E TRAINING DETAILS FOR THE TWO BASELINES
text: Scenario-M This method aims to generate one scenario, which can result in a good approximate solution. It is not easy since for each instance, sufficient scenarios need to be evaluated to attain the target scenario for supervised learning. To avoid much computational burden, Bengio et al. (2020) compute the average of scenarios as the initial one and then iteratively improve it to derive a similar solution as the optimal one. However, the improvement heuristic used cannot guarantee that the attained one exists in the original set of scenarios. To tackle this issue, we adapt the method to track the improvement along existing scenarios. Specifically, we find the 20 scenarios nearest to the average and compute the MIPs defined by each of them with the context. We determine the target scenario by picking the one with the most similar solution to the optimal one. This method could be used for both NDP and FLP.
We conduct the same featurization as described in Appendix A.4 in (Bengio et al., 2020), and use both linear regression (LR) and artificial neural network (ANN) with MSE loss. In addition, we also use the same GCN in our method to directly process raw parameters in SIP instances. Specifically, we only input the context with the average of scenarios to GCN, and the graph embedding is linearly projected into output to estimate scenarios. The reported results in the main paper are delivered by GCN, which we find outperforms LR and ANN in most cases.
Solution-M Following (Abbasi et al., 2020), we directly predict the first-stage solution in both NDP and FLP. Despite the machine learning algorithms listed in the paper, we use the same GCN in our method for training. It is more powerful to process the raw parameters under the graph representation. In our cases, the raw parameters are the static parameters in the context for both NDP and FLP. We take them as the input into GCN and the output dimension is same as the number of first-stage decision variables. If there is no static parameters on nodes or edges, we use 0 for them to run GCN as usual.
For both supervised methods, we use 12,800 instances (of normal size) in the training with the number of epochs 100, batch size 64 and learning rate 0.001. We test them with diverse instances as described in corresponding experiments in the main paper.

heading: F ADDITIONAL RESULTS ON LARGE FLP
text: We evaluate our methods on FLP of large sizes to show the trend of the average error as the number of representative scenarios grows. We follow the experimental setting in Section 5.3 and apply the trained model to solve 50 instances from FLP (60; 200) and FLP (120; 200). We set K from 1 to
20 and thus attain solutions with varying numbers of representative scenarios. The curve of average approximation errors is plotted in Figure 2. We observe that both CVAE-SIP and CVAE-SIPA persistently reduce errors towards 0 as K grows. It means that our methods are consistently attaining more effective representatives from scenarios, which help deliver better solutions. Meanwhile, CVAE-SIP and CVAE-SIPA outperform K-medoids with different K on both problems, especially on FLP (120; 200). It again verifies the merit of the learned representations for solving SIP problems. The comparison between CVAE-SIP and CVAE-SIPA implies that though CVAE-SIPA is superior to CVAE-SIP at the early stage of the curve, i.e., K=1,. . .,3, they are on par with each other when more representatives are used eventually.

heading: G GENERALIZATION ACROSS DISTRIBUTIONS AND DEPENDENCIES
text: We test the generalization performance of our method on NDP with scenarios sampling from the Normal distribution and Binomial distribution. For the former, we set the expectation and standard deviation to 10 and 2, respectively. For the latter, we set the number of runs and the probability to 20 and 0.5. The results are summarized in the upper half of Table 5. We observe that our methods consistently outperform other learning based approaches with scenarios from both distributions. When K=10, both CVAE-SIP and CVAE-SIPA achieve near-optimal solutions with slight approximation errors. On the other hand, it can be found that the results delivered by CVAE-SIPA are generally superior or comparable to those by CVAE-SIP, suggesting a better cross-distribution generalization.
In our training, we do not assume explicit dependencies between the context and scenarios, following the settings in (Keutchayan et al., 2020; Nair et al., 2018; Abbasi et al., 2020). However, we would like to show that the trained networks can well generalize to the dependent settings. In specific, we construct in NDP two dependencies between the demand scenarios and costs in contexts as follows:
1) D0: for each instance, the lower bound of Uniform demands is determined by the smaller value of the average opening and shipping cost in the context; 2) D1: the lower and upper bounds of Uniform demands are determined by the smaller and larger value of the average opening and shipping cost, respectively. We generate 50 instances for each dependency. From the lower half of Table 5, we observe that our methods outperform learning based baselines except the cases with K=3, where CVAE-SIP loses its advantage over a few baselines. Generally, it indicates that our methods possess fairly good generalization abilities to the setting where scenarios depend on the contexts.

heading: H RESULTS ON OBJECTIVE PREDICTION
text: We first evaluate our trained model for the prediction task in three FLP instances, which are randomly generated from FLP (30; 200), FLP (120; 200) and FLP (30; 800), respectively. For each instance, we draw a scatter plot of the 2-dimensional latent representations for scenarios and represent the normalized objective values using different colors. The results are shown in the upper half
of Figure 4, and the corresponding ground truth (i.e. optimal solution values attained by CPLEX) are displayed in the lower half. It is clear that the prediction is fairly accurate with almost the same patterns as the target values in the ground truth, even for the larger problem size or scenario cardinality. It indicates that our design for semi-supervised learning is effective and the trained model exhibits a desirable generalization ability.
We further evaluate the objective prediction on NDP. Specifically, the trained model is used to predict objective values for the scenarios from three instances, which are randomly generated from NDP (14; 200), NDP (44; 200) and NDP (14; 800), respectively. From the results in Figure 3, we observe that our method can also predict the objective values accurately for NDP. Especially, the pattern of predicted values for NDP (14; 800) is almost indistinguishable from the one by ground truth. It reveals that our design for the semi-supervised learning is effective across distinct SIPs.

heading: I COMPARISON WITH OR METHODS
text: In this section, we compare our methods with two mature methods from operations research (OR), i.e., 1) Scen-red (Heitsch & Römisch, 2003), which seeks a subset of scenarios that is closest to original scenarios in terms of Fortet-Mourier probability metric, and 2) CSSC (Keutchayan et al., 2020), which solves a MILP to derive representative scenarios, each of which estimates the average objective of some other scenarios. They are on behalf of the distribution-oriented and problemoriented paradigm, respectively. The tested instances are the same as the ones in Section 5.2 and 5.3. From the results in Table 6, we observe that our methods are comparable to or better than the baselines and especially on NDP. While CVAE-SIP and CVAE-SIPA are inferior to CSSC on NDP with K=5, they attain the best approximate solutions when K=10 and 20. Notably, the errors of our methods are close to 0 on both normal and large instances when K=20. For FLP, while our methods are inferior to CSSC and Scen-red on FLP when K=5 and 10, CVAE-SIPA surpasses CSSC on FLP (30; 200) and gains the smallest objective value and error on FLP (120; 200), when K=20. Since our runtime is shorter than CSSC and similar to Scen-red, the results again verify that our methods are effective to reduce scenarios for good approximations.

heading: J ABLATION STUDY
text: We further conduct an ablation study on NDP to verify the effectiveness of the learned scenario representations, and meanwhile investigate the mixed use of these representations with the original scenarios. In specific, we enrich the original scenario set for clustering in the following ways: 1) CVAE-SIPA+Ori.: we concatenate the predicted objective values with their original scenarios (rather than the learned representations in CVAE-SIPA); 2) CVAE-SIP+Ori.: we concatenate the learned representations with their original scenarios (rather than merely use the learned representations in CVAE-SIP). These two kinds of mixed representations are used for clustering to find representative scenarios, where we set K=5, 10, 20. The tested instances are the same as the ones in Section 5.2 and 5.3. As shown in Table 6, both CVAE-SIP and CVAE-SIPA are generally superior to their
counterparts. The advantage of CVAE-SIPA over CVAE-SIPA+Ori. indicates that with the same predicted objective values, our learned representations are effective to attain better representatives than the original scenarios. Moreover, CVAE-SIP is superior to CVAE-SIP+Ori. probably because the learned representations are overwhelmed by the original scenarios to some extent. It is also interesting to identify an advantage that CVAE-SIP > K-medoids > CVAE-SIP+Ori, which means the pure use of learned representations or original scenarios is better than the mixed use of both.

Review summaries:
This paper studies the problem of generating representative scenarios for two-stage stochastic integer programmings in which the parameters could be either static (referred to as context) or stochastic (forming a space of scenarios). The proposed method leverages conditional variational autoencoder to learn the distribution of scenarios conditioned on the context, where the encoder computes the latent variables (based on the context and scenarios) while the decoder recovers the scenarios. In particular, the functions are parametrized by graph neural networks and multilayer perceptron. Experiments on two concrete problems are performed. 

The paper aims to use deep-learning towards solving generic stochastic integer programs. The proposed approach involves using a conditional VAE to model the distribution over scenarios, and leveraging the scenario's latent space embeddings to enable better representative scenario selection: the embeddings are used as inputs to a standard clustering algorithm that returns representatives.  The paper demonstrates that using VAE embeddings improves performance when compared to using the same clustering algorithm on the raw scenario representations. It further demonstrates that using embeddings concatenated with a predicted objective generally improves performance over using the embeddings alone. Finally, the paper shows improvement over two recent approaches.

Solving stochastic integer programs in practice can be difficult for large problems with many scenarios (i.e., sample discretizations of the uncertain distribution), motivating relevant problems such as scenario reduction and objective prediction. This paper explores these problems by using a graph conditional VAE to learn low-dimensional representations of scenarios and perform downstream tasks in this space. The results show improvements over existing scenario reduction and objective prediction methods over the original high-dimensional space.