XOR-Sampling for Network Design with Correlated Stochastic Events

1 Introduction

Many problems such as influence maximization [Kempe et al., 2003], river network design [Neeson et al., 2015] and predisaster preparation [Peeta et al., 2010] can be formulated as stochastic network design problems, in which the edges of the network are present or absent stochastically. The goal is to take protection actions, which reinforce the connectivity of a subset of edges with a limited budget, to maximize the connectivity of the network under stochastic events. For example, in the road network design problem, as shown in Figure 1, edges can be present or absent from the network stochastically when a natural disaster strikes. The problem then is to invest smartly under a limited budget to strengthen a few key road segments, so that ambulance centers can still remain connected to population centers under the attacks of multiple natural disasters.

∗ equal contribution.

Despite the worst-case PP-hard complexity to evaluate the connectivity values under a fixed protection plan [Dyer and Stougie, 2006], approximate algorithms based on Sample Average Approximation (SAA) have been developed to find nearly optimal solutions [Sheldon et al., 2010; Wu et al., 2015], which convert the original stochastic optimization problem into a deterministic one, by optimizing the connectivity objectives over a fixed set of sample scenarios.

The success of previous approaches depends on the core assumption that the state of each edge (e.g., present or absent) is independent of those of other edges. Under this assumption, high-quality samples can be drawn with simple schemes such as flipping a biased coin for each edge. Also, it has been shown that often a small number of samples (e.g., 10) is sufficient to find a good policy [Sheldon et al., 2010], despite the fact that a much larger number of samples is needed for theoretical guarantees [Kleywegt et al., 2002].

Nevertheless, the states of edges are often correlated in practice. For example, in the road network design problem shown in Figure 1, edges tend to fail concurrently when a regional event (e.g., a natural disaster such as an earthquake or a snow storm) disables the connectivity of all edges in a given area. In this paper, we study this more realistic setting where the states of multiple edges can be correlated. We propose to use a general probabilistic model - Markov Random Field (MRF) - to model the joint distribution of edge states, which can be built with domain knowledge or learned from data.

We make several contributions: (1) We propose a new stochastic network design problem framework, taking into account correlated stochastic events based on MRF; (2) We extend the SAA framework with a novel sampling scheme to convert the original stochastic network design problem into a deterministic network design problem, which is formulated as a Mixed Integer Program (MIP) and solved by CPLEX 12.6. As we will show, also mentioned in paper [Kumar et al., 2012], the MIP can only scale to a small number of samples (e.g., 30 samples). Consequently, a small but good quality set of samples is critical to produce high-quality solutions; (3) Another key contribution of this work is the coupling of XOR-sampling [Gomes et al., 2007; Ermon et al., 2013] with SAA; (4) We show that XORsampling coupled with SAA outperforms the standard Gibbs sampler when coupled with SAA. As shown in paper [Ermon et al., 2013], the XOR sampler can produce MRF samples with better quality but often has worse scalability than the Gibbs sampler. However, in our setting, the bottleneck of the runtime is in solving the MIP, not the time required by the XOR sampler. This is because the MIP solver can only handle relatively small samples (therefore they have to be carefully selected). Given the relatively small sample sizes, the time required by the XOR sampler to produce such samples is negligible. (5) We test our algorithms on the flood preparation problem with real-world road networks [Wu et al., 2016]. The goal is to maintain the connectivity of road networks for emergency medical services (EMS). The results show that SAA with XOR sampling produces policies that are never inferior in quality compared to SAA with Gibbs. In other words, the policies produced by the SAA with XOR sampling consistently have similar or better quality than SAA with Gibbs sampling. Furthermore, in some cases, a small number of samples (e.g., 10) from the XOR sampler, is enough for the SAA to converge. Also, the solution quality with XOR sampling is consistently good, not very sensitive to the variance of samples used in SAA. In comparison, the solution quality produced by a Gibbs sampler has a higher variance when the number of samples used in SAA is small. These observations suggest that, on large-scale networks, the SAA method can still produce a good policy if we use high-quality samples such as those from XOR sampling.

2 Problem Statement

One important application of the stochastic network design problem is road network design [Peeta et al., 2010]. Natural disasters such as earthquakes and floods break road segments and destroy the connectivity of the whole road network. The paper [Wu et al., 2016] introduces the problem of flood preparation for EMS as shown in Fig. 1. In this problem, during floods, roads may be washed out and paths connecting ambulance centers to patients are cut off. To maintain good connectivity of the road network in case of an emergency, spatially critical roads should be identified and protected beforehand, at certain costs. The question is how to find these spatially critical roads, given a budget limit, such that most of the patients can be reached when floods happen. This problem is an example of the stochastic network design problem.

A stochastic network design problem is defined on a directed graph G = {V,E} with a source node s. Each node is associated with a weight w(v) (e.g., population density). Each edge e ∈ E is associated with a binary random variable θe that describes the state of the edge during disasters. θe = 1 means that the edge e is present or not destroyed. θe = 0 means that e is absent or destroyed. Let θ = {θe, e ∈ E} denote the vector of random variables of all edges, each random variable taking value in {0, 1}|E| and Pr(θ) denote their joint probability distribution. Pr(θ) may be constructed from domain knowledge or learned from realworld data. Note that we use the bold symbol θ for random variables and the unbold symbol θ for their values. For a given θ, the indicator function I(s ⇝ v|θ) indicates whether node v is reachable from the source s in the network by a path of all present edges. If node v is reachable, it equals to 1. Otherwise, it equals to 0. The expected weighted number of nodes that are reachable from the source s is written as

To be general, one protection action can protect one or multiple edges. If an edge is protected, it is always present. Formally, we are given a collection of M edge sets E = {E1,E2, ...,EM} where each Es ⊆ E represents an action to protect all edges in Es and is associated with a cost cs. A policy π is a subset of E that indicates all edges being protected. Protecting an edge from u to v can be equivalently viewed as adding an edge (u, v) that is always present in the network. That is, if (u, v) is protected by π, u is connected to v no matter θ(u,v) = 1 or 0, which makes the indicator function depend on π as I(s ⇝ v|θ, π). In other words, protection actions only change the network structure but not Pr(θ), the probability of the stochastic events such as natural disasters that is independent of the network structure. The decision making problem is to find the policy π, subject to a budget limit B, to maximize the expected weighted number of nodes that are reachable from s, that is

This problem definition is general and subsumes many special cases. For example, in the real-world problem in our experiments, only a subset of edges are present or absent stochastically while the remaining edges are known to be present. For this special case, we can set the protection costs of those edges that are always present to zero so they are protected automatically in the optimal solution.

Multiple Source Case In problem (2), it is assumed that there is only one source. It can be extended to multiple sources. For example, there are several ambulance centers, each of which is represented by a source. The house of a patient is encouraged to be reachable from multiple centers in case that ambulences in one center are busy and can't respond. In this new problem formulation, instead of having a single source s, we are given a set of sources S. Connecting to more sources increases the objective value accordingly. For each source in S, we count the weighted number of nodes that are reachable from the source separately. The objective is to maximize the expected total weighted number of nodes that are reachable from each source in S, that is,

In this case, each node in S issues a separate network flow to each other node. Later, we will introduce the algorithm to solve single source network design problems and it can be extended to solve multiple source problems.

Model Pr(θ) We propose to use the Markov random field (MRF) to model the joint probability distribution Pr(θ) of all variables. MRF is a general model for the joint distribution of multiple correlated random variables. Using notions in MRF, we define the probability of edge states Pr(θ) as:

where θα is a subset of variables in θ that the function ψα depends on. ψα : {θ}α → R+ is a potential function, which maps every assignment of variables in {θ}α to a nonnegative real value. I is an index set. Z is a normalization constant, which ensures that the probability adds up to one: Z = ∑θ∈{0,1}|E| ∏α∈I ψα({θ}α). In this paper, we will also use ˜ Pr(θ) to represent the unnormalized probability density:

A potential function ψα({θ}α) defines the correlation between all variables in the subset {θ}α. As a simple example, suppose {θ}α contains the subset of variables for all network edges within a small spatial region. All roads in set {θ}α may be either all present or all absent, depending on whether a regional event (e.g., an earthquake) happens or not. In this case, ψα({θ}α) will take high values when all variables are 1 or 0 simultaneously. The structure of the MRF or the set I can be built from domain knowledge and potential functions can be learned from real-world data. The focus of this paper is not on how to construct the MRF but is how we solve problem in Equation 3 in general when Pr(θ) is given in the form shown in Equation 4.

Theorem 1. For given policy π, evaluating the objective of problem (3) is #P complete.

Proof. This is a consequence of the #P-completeness of the graph reliability problem [Dyer and Stougie, 2006], which is to compute the probability that two given vertices is connected if the edges are present or absent independently. Our problem simplifies to the graph reliability problem when all potential functions are singular.

Figure 2: A compact MIP encoding of problem (5).

Other Applications The stochastic network design problem (2) can be used to model other network optimization problems in a more realisitc way by allowing correlation between random events. Here are two examples.

Barrier removal problem [Wu et al., 2014a; 2017] is an important problem in computational sustainability where the goal is to maximize the expected amount of habitat that one randomly selected group of fish can reach. Previous problem formulation assumes that the probability that a group of fish is capable of passing a barrier is independent of the probabilities of passing other barriers. However, in practice, the probability for a group of fish to pass a barrier mainly depends on the structure of the barrier and the capability of the fish. Therefore, it is more realistic that a fish is able to pass multiple structurally similar barriers or pass none of them. Our algorithm introduced later can be applied to solving the barrier removal problem with correlated passage probabilities as well.

Influence maximization problem [Kempe et al., 2003] is a well-known problem in social computing, in which people who are infected by an influence will continue to infect their friend stochastically. The independent cascade model assumes that each person will infect another person with some probability independent from infections from other sources. Nevertheless, the infections coming frommultiple friends can be correlated, since friends from the same social group tend to post messages of similar content.

3 Our Method

In this chapter, we introduce the algorithm to solve the single source network design problems (2) with Pr(θ) defined by (4). The algorithm can be extended to solve the multiple source network problem (3) with a little modification to the MIP formulation. We use the Sample Average Approximation [Kleywegt et al., 2002] (SAA) method to solve problem (2). The basic idea is to draw N samples {θ1 , ...., θN } of θ from the distribution Pr(θ) and convert the stochastic optimization problem (2) into the following deterministic optimization problem:

Theorem 2. Let π∗ be the optimal policy of problem (5) with N samples and πOPT be the optimal policy of problem (2). The number of samples for π∗ to converge to πOPT is Θ(|E|) [Kleywegt et al., 2002].

Problem (5) can be formulated as a MIP. Using the maxflow encoding from one source to multiple targets, we introduce a more compact MIP formulation than the one in paper [Wu et al., 2015] which uses O(|V |) times more variables (|V | is the number of nodes). This new compact formulation is shown in Fig. 2. The basic idea is as follows. For each sample θi , and for each node v, we define a binary variable ziv , which represents whether node v is reachable from s (= 1) in the ith sample or not (= 0). Variables ziv can be relaxed to be continous variable as shown in Fig. 2 because the value of ziv is either 0 or 1 in the optimal solution. The objective function in Equation 5 can be written as:

We then defines an indicator variable yEs for the action to protect all edges in Es. Protecting edges in Es costs cs, therefore we have the budget constraint (5) in Fig. 2.

Next, we enforce flow constraints (2,3,4 in Figure 2) to guarantee that zik = 1 if and only if node k is reachable from s in the ith sample. The idea is to inject ∑k̸=s zik unit of flow into the source s, and to pull out zik unit of flow (zik is either 1 or 0, depending on whether node k is reachable from s in the ith sample or not) at all nodes other than the source s. Here we use xie as flow variables to represent the flow on edge e in the i-th sample. Constraint (2) is used to enforce a flow of this kind. If a node k is reachable from s (zik = 1), one unit of flow will be consumed by node k. If k is not reachable from s (zik = 0), k will not consume any flow. Constraint (4) restricts the amount of flow on each edge between 0 and |V | - 1, which must be true since the total amount of flow is ∑k̸=s zik and we have ∑k̸=s zik ≤ |V | - 1. Constraint (3) ensures that no flow is allowed on edge e unless e is present in the i-th sample (θe = 1) or e is protected by at least one action Es that covers e and the corresponding yEs = 1.

SAA is sensitive to the number of samples N used to approximate the objective function. Although in general a large N is preferred in order to reduce the variance, the number of variables in the MIP encoding in Figure 2 scales linearly with respect to N, which prevents us from using a large number of samples (e.g., millions of samples) in practice. Under the independent edge sampling case, despite theoretical results [Kleywegt et al., 2002] claim that O(|E|) number of samples are required to obtain probabilistic guarantees, in practice, people found that a small number of samples (e.g., 10) are sufficient to find a good policy [Sheldon et al., 2010].

We call this small set of samples as good samples or wellrepresentative samples. That is, the number of samples is relatively small but these samples can be used to find good policies. To make SAA work, it is crucial to obtain a smallsize (N small) but well-representative set of samples, so that not only the quality of the sample approximation is good, but also the size of the sample set is small enough for MIP to scale. It is unclear and potentially more difficult to obtain good samples when edges are correlated. In the following sections, we introduce two candidate methods to attempt to obtain these so-called good samples {θ1 , ..., θN }.

3.1 Gibbs Sampling

Gibbs sampling is a special case of the Metropolis-Hastings algorithm to obtain a sequence of samples which approximates a given distribution Pr(θ). The high-level idea is to build a Markov chain whose transition probability is

It has been proved that the stationary distribution of the Markov chain converges to the target probability distribution Pr(θ). Notice that we only need the unormalized density ˜ Pr(.), rather than Pr(.) to compute Pr(θi|θ1, .., θi-1, θi+1, .., θ|E|). This is the crucial fact that makes the transition of Gibbs sampling efficient in practice.

3.2 XOR Sampling

Sampling with XOR constraints [Gomes et al., 2007; 2006; Ermon et al., 2013; Chakraborty et al., 2013] was developed, thanks to the rapid progress in solving NP-complete problems in the last decade. Despite the intractability in the worst case, modern MIP and Satisfiability Testing (SAT) solvers are capable of exploiting the problem structure, and can scale up to millions of variables [Pipatsrisawat and Darwiche, 2010; Sontag et al., 2008; Riedel, 2008; Poon and Domingos, 2006; Gogate and Dechter, 2011]. XOR sampling transforms the problem of sampling, which is an #P-complete problem, into solving satisfiability problems with additional XOR constraints.

With XOR sampling, the resulting sampling distribution can be proved to be arbitrarly close to uniform in the unweighted case [Gomes et al., 2007; Chakraborty et al., 2013], and to be within a constant factor of the true probability distribution in the weighted case [Ermon et al., 2013]. This is a significant improvement compared with MCMC based techniques, which often require an exponential amount of samples to reach its stationary distribution on many problems with intricate combinatorial structure.

XOR sampling is a special fit to our SAA framework. The rationale behind XOR sampling is to trade computational efficiency with provable near-optimal quality in the samples obtained. XOR sampling is more expensive than Gibbs sampling because it has to solve combinatorial problems at each sampling step. The recent progress in constraint programming makes it possible to solve large-scale NP-complete problems efficiently, therefore the benefit often outweights the cost. This is especially true in our SAA framework. First, the MIP encoding for SAA prevents us from optimizing over millions of examples simultaneously. Therefore, the samples we use in the optimization has to be small, but representative enough for SAA to balance over all likely scenarios. Second, since the samples are obtained prior to the SAA optimization, we can afford to spend some time on obtaining good samples. We can even run several sampling programs in parallel. In experiments, we found that the time spent at the sampling step was negligible compared with solving the MIP.

General Idea To illustrate the high-level idea of XOR sampling, let us first consider the unweighted case where the unnormalized density ˜ Pr(θ) is either 1 or 0. Suppose there are 2l different θ assignments which make ˜ Pr(θ) = 1. We call all these assignments "satisfying assignments". An ideal sampler draws samples from the set of satisfying assignments ∆ = {θ : ˜ Pr(θ) = 1} uniformly at random; i.e., each member in ∆ has 2-l probability to be selected.

At a high level, XOR sampling obtains near-uniform samples by finding one satisfying assignment subject to additional randomly generated XOR constraints. To be specific, we keep adding XOR constraints to the original problem of finding satisfying assignments. We can prove that in expectation, each newly added XOR constraint rules out about half of the satisfying assignments at random. Therefore, if we start with 2l satisfying assignments in ∆, after adding l XOR constraints, we will be left with only one satisfying assignment in expectation, by which time we return this assignment as our first sample. Because we can prove that the assignments are ruled out at random in each step, we can guarantee that the last assignment left must be a randomly chosen one from the original set ∆. Please see [Gomes et al., 2007] for details.

For the weighted case, our goal is to guarantee that the probability of sampling one θ is propotional to the unnormalized density ˜ Pr(θ). We implement a horizontal slice technique to transform a weighted problem into an unweighted one, taking a 2-approximation. For the easiness of illustration, consider the simple case where P0 = minθ ˜ Pr(θ), Pk = maxθ ˜ Pr(θ) and Pk = 2k P0. Let δ = (δ0, . . . , δk-1) ∈ {0, 1}k be a binary vector of length k. We sample (θ, δ) uniformly at random from the following set ∆w using the afromentioned unweighted XOR sampling:

If we sample (θ, δ) uniformly at random from ∆w and then only return θ, it can be proved that the probability of sampling θ is propotional to P02i-1 when ˜ Pr(θ) is sandwiched between P02i-1 and P02i . We refer the readers to [Ermon et al., 2013] for details.

4 RelatedWork

Our work is related to an active research field on reliability of physical communication networks under natural disasters and malicious coordinated attacks [Rahnamay-Naeini et al., 2011; Neumayer and Modiano, 2010]. In those papers, probabilistic models are built based on spatial point processes to evaluate the reliability of networks with geographically correlated failures. However, the problem of optimizing the reliability of networks by taking protection actions is not discussed and the proposed models are too complex to be used to define solvable decision-making problems for network protection.

Sheldon et al. (2010) formulate the problem of maximizing the spread of cascades as stochastic network design and develop an algorithm combining SAA and MIP. In their model, all variables are assumed to be independently distributed so samples of each variable can be drawn independently. Their MIP works only when the network is a DAG and has only one source while our MIP works for general directed graphs and an arbitrary number of sources. Wu et al. (2014b) solve the river network design problem with a single source and multiple sources by rounded dynamic programming algorithms [Wu et al., 2014a; 2014b]. However, their algorithms only work for tree-structured networks, and they also assume that stochastic events are all independent. Recently, several fast approximate algorithms based on the primal-dual technique have been developed to solve several deterministic network design problems similar to our problem (5) [Wu et al., 2015; 2016], but they do not guarantee to produce the optimal solutions. It is an interesting future work to see how we can combine their algorithms and XOR sampling techniques to tradeoff the optimality with efficiency on large-scale networks. Xue et al. (2017) use an optimization algorithm by adding XOR constraints to solve a related landscape connectivity problem. Their optimization problem takes into account of spatial capture re-capture model, which is in a different setting from ours. More importantly, we use XOR sampling within the SAA framework, which is different from their approach. 0

5 Experiments

We test our algorithm on a real-world problem, the so-called Flood Preparation problem for the emergency medical services (EMS), on road networks introduced in paper [Wu et al., 2016]. Edges of the graph represent road segments. Nodes represent either the intersections of roads or locations where EMS are requested through phone calls, which are sampled based on the patient population density. Some edges or road segments are above roadway stream crossing structures, such as culverts, which are vulnerable to floods and earthquakes. The failure of a crossing structure will disrupt the road segment above it and cause the road to become impassable. In our study area, only a subset of edges are associated with crossings, each of which defines a binary variable θe. The value of θe indicates whether the crossing under e fails (= 0) or not (= 1). A protection action can replace a crossing with a sturdy type such that it will not fail during floods. The goal is to decide which crossings (or which edges out of the subset of edges above crossings) to protect beforehand so patients remain connected to ambulance centers even during floods, that is, we will decide for each edge above a crossing whether to protect it or not. To measure the connectivity during floods, we use the objective of the multiple source case introduced in section 2. Each source represents an ambulance center. Since each center usually has a limited number of ambulances available (sometimes only one ambulance available), a location reachable from multiple centers has a lower risk of no ambulance responding during floods than a location reachable from only one center. The problem can also be treated as a single source problem in which it is sufficient to connect a node to only one center or being able to connect to more than one center doesn't further increase the objective. While our methodology applies to both cases, we only show the results for the multiple source case.

We use a similar way as in paper [Neumayer and Modiano, 2010; Rahnamay-Naeini et al., 2011] to synthetically generate failure correlation between these crossing structures. The idea is as follows. Each crossing may fail with a small chance independently and multiple spatially nearby crossings tend to fail concurrently if affected by a natural disaster. To generate disaster locations, instead of using spatial point processes as in paper [Neumayer and Modiano, 2010] due to the lack of landscape features, we uniformly generate a small number of disaster centers and define their effect regions with randomly generated radiuses. All variables within a region are correlated with different levels of strengths (strong or weak).

Small Networks We test our algorithm on a small network of 502 edges and 20 crossings extracted from the road networks, which consists of 20 binary random variables. 20 is small enough to compute the exact objective value of Equation (3) for any given policies so we can compare different policies exactly. We test on a different number of sources.

For each problem instance, we vary the sample size from 10 to 180 and budget size from 10% to 40% of the total cost to protect all crossings. For sample size N, we produce 10 policies, each computed with N samples independently drawn from the Gibbs and XOR sampler. The objective value of each policy is computed with an exact solver. The mean and variance over 10 policy values are reported. The results of four "representative" instances are shown in Fig. 3, and the results of other experimental instances show similar trends. In most of the sample sizes and testing instances, the mean values produced by SAA with the XOR sampler is as good or better than those obtained by SAA with the Gibbs sampler. Moreover, the policy produced by the XOR sampler is more stable in the sense that the variance of policy values is smaller than the variance produced by the Gibbs sampler, as shown in Fig. 3(a-c). The results by the Gibbs sampler are comparable to the results by the XOR sampler in some instances. In Fig. 3(d), both methods produced almost identical means and variances. In Fig. 3(a), the XOR sampler gives smaller means with 10 samples but quickly outperforms the Gibbs sampler with more samples. In Fig. 3(a)(c), for the 10% budget cases, the Gibbs sampler produces larger means but appears very unstable. In most of these results, as the number of samples increases, the variance by the XOR sampler quickly reduces to 0 while the variance by the Gibbs sampler remains large (e.g., in Fig. 3(b) and (c)).

Large Networks We test on a large road network of 1, 562 edges and 81 binary random variables where the exact policy value is not available. Fig. 4 shows that the MIP quickly becomes unscalable as the sample size and the number of sources increase. CPLEX takes around 2.5 hours to solve the MIP of 30 samples and 20 sources and cannot finish within 4 hours for 50 samples and 20 sources. Therefore, to work with 20 or more sources or networks of larger sizes, it is important that good policies can be obtained by a small number of samples (e.g., 30 samples for 20 sources). For 10 source case, it is interesting that the runtime decreases a little bit after a certain number of samples for some budget sizes.

Also, for both samplers, we do experiments to find out the number of samples that are sufficient for the sample average method to obtain a reasonably good estimation of policy value. To test whether a sequence of values converges, we borrow the idea from Cauchy sequence test to see if the difference of two subsequent values converges to zero. The results are shown in Fig. 5. We create a sequence of sample sizes {N1,N2, ...} varying from 10 to 5000. In Fig. 5, the value reported on y-axis at x = Ni equals the difference between the average of Ni policy values and the average of Ni-1 policy values, each policy value estimated by an independent drawn sample of θ. If the difference goes to zero, it implies the sequence converges. As one can see from the figure, the value difference for both the XOR and the Gibbs sampling narrows down to a small interval while the value difference for the XOR sampler converges to a much smaller range near 0 and much faster than for the Gibbs sampler. The results also imply that the empirical mean of 5000 samples can produce good estimations of the policy value.

At last, we compare the solution qualities of our algorithm on large networks. Two representative results are shown in Fig. 6. Again, the results of instances that are not shown in the figure have similar trends. Since we can not calculate the exact value of a policy (e.g., an expectation) on the large network, we estimate the expectation using 5000 empirical values of the policy, each calculated by an independent sample. As demonstrated, 5000 samples are sufficient to obtain a good estimation. In the figure, the XOR sampler outperforms the Gibbs sampler with better and stable policy values, in a similar way as for the small networks. In conclusion, the XOR sampler produces solutions with higher qualities and is less sensitive to sample variance.

6 Conclusion

Previous approaches for solving stochastic network design problems rely on a restricted assumption that the states of edges are independently distributed. An SAA based algorithm can produce high-quality solutions by taking advantage of this restricted assumption, although they only scaling up to a small number of samples. In this paper, we relax the independence assumption and propose a more general problem framework by allowing correlated stochastic events to affect the states of edges in networks. Novel algorithms have been developed based on the sample average method and various samplers. We show that our algorithm can produce high-quality solutions with a small number of samples if these samples are obtained from a XOR sampling technique. Testing on a real-world flood preparation problem for EMS, empirical results show that SAA with the XOR sampler produces better and more stable solutions compared with SAA with a Gibbs sampler.

Acknowledgements

This research was supported by the National Science Foundation under Grant No. CNS-0832782, CNF-1522054, CNS1059284, ARO grant W911-NF-14-1-0498 and the Atkinson Center for a Sustainable Future.