12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015Network Reliability Analysis for Cluster Connectivity UsingAdaBoostRaphael E. SternGraduate Student, Dept. of Civil and Environmental Engineering, University of Illinoisat Urbana Champaign, Urbana, IL, USAJunho SongAssociate Professor, Dept. of Civil Engineering, Seoul National University, Seoul, KoreaDaniel B. WorkAssistant Professor, Dept. of Civil and Environmental Engineering and CoordinatedScience Laboratory, University of Illinois at Urbana Champaign, Urbana, IL, USAABSTRACT: In the aftermath of a natural disaster, knowledge of the connectivity of different regions ofinfrastructure networks is crucial to aid decision makers. For large-scale networks it can be extremelytime-consuming to obtain a converged estimate by performing a large number of Monte Carlo simulationsto compute the network failure probability. To reduce computational requirements, this work develops asurrogate model using an AdaBoost classifier for predicting probabilities of disconnections between nodeclusters in lifeline infrastructure networks. The proposed approach uses spectral clustering to partitionthe network, and it estimates the connectivity of these clusters using an AdaBoost classifier. Numericalexperiments on a California gas distribution network demonstrate that using the surrogate model to de-termine cluster connectivity introduces less than five percent error and is two orders of magnitude fasterthan methods using an exact network model to estimate the probability of network failure through MonteCarlo simulations.1. INTRODUCTIONFollowing a natural disaster, fast response is criticalto minimize the loss of life and damage to infras-tructure. In lifeline networks such as gas distribu-tion networks, power grids, and water pipelines, theability to rapidly assess the connectivity of compo-nents within the network is vital for prompt disasterrelief (Bruneau et al., 2003; Boin and McConnell,2007). This assessment requires methods to quicklyand accurately estimate the probability of networkfailure given the individual component failure prob-abilities conditioned on the event.A variety of system reliability analysis (SRA)methods have been introduced to identify the prob-ability that a network will remain functional in theaftermath of an event. One set of approaches, forexample those discussed by Rausand and Høyland(2004), aim to exactly compute the network fail-ure probability from individual component failureprobabilities. Lim and Song (2012) introduced amethod that intelligently enumerates all possiblefailure combinations by preferentially identifyingdisjoint cut sets and link sets to calculate the prob-ability of network disconnection, while Reed et al.(2009) and Vugrin et al. (2010) proposed to analyzeinfrastructure resilience to natural disasters moregenerally.While these methods are useful for small net-works when precise failure probabilities are knownin advance, they have several limitations for appli-cation to large, infrastructure-sized networks. Asthe number of nodes in the network increases, ex-ploring all possible node failure combinations thatlead to network failure can quickly become com-putationally intractable. Furthermore, before theearthquake event, one can compute the failure prob-112th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015abilities of components for each possible earth-quake scenario, or compute the aggregate failureprobability for all possible earthquakes representedby a probabilistic seismic hazard model. Both ofthese approaches have inherent uncertainty sincethey require conditioning on a specific earthquakeevent. Knowledge of the precise component failureprobability for a specific event are not known untilafter the event occurs, which means that for resultswith the least uncertainty the analysis may need tobe performed in the immediate aftermath of a dis-aster, when time is critical.Because of the above limitations, there is re-cent interest to develop approximate methods tocompute network failure probabilities, or boundson such probabilities (Der Kiureghian and Song,2008; Song and Der Kiureghian, 2003). For exam-ple, Monte Carlo simulations (MCS) are frequentlyused to estimate the probability of network failurein SRA (Papadrakakis et al., 1996; Ditlevsen andMadsen, 1996). However, the Monte Carlo tech-niques still require one to determine the connectiv-ity of nodes in a network for each realization of thenetwork, which can be very time consuming if reli-able estimates are desired.To further speed up the calculation time, Sternet al. (2014) proposed the idea of using a surrogatemodel to determine the connectivity of two nodepairs, instead of testing the connectivity on the ac-tual infrastructure network. The surrogate model isa simplified model that approximates the actual net-work model of interest and on which calculationscan be performed much more quickly. Surrogatemodels are commonly used to simplify computa-tions in fields as diverse as structural optimization,waste water modeling, and supply chain manage-ment (Jansson et al., 2003; Meirlaen et al., 2001;Wan et al., 2005).Stern et al. (2014) used MCS to approximate theconnectivity of a single source and terminal node,using an AdaBoost (Freund and Schapire, 1995)classifier as the surrogate model. This model istrained on a large labeled dataset of network statesin an offline pre-processing stage, and the samplesare chosen independently of any specific node fail-ure probabilities. As a result, the trained modelcan be used in the immediate aftermath of an eventwhen the node failure probabilities are realized, byperforming an MCS on the surrogate model. Theresults in Stern et al. (2014) indicate the approxi-mate model can estimate the probability of discon-nection of a source and terminal node to within 3%under worst-case conditions, and is six times fasterthan using a shortest path method to determine net-work failure in the MCS.Due to the hierarchical structure and regional di-visions of emergency management authorities, de-cisions regarding disaster response often are madeon a regional basis. Therefore, for large-scale riskassessment on infrastructure networks it is impor-tant to determine the probability of disconnectionof these large regions of the network, or node clus-ters. In this work, we extend the AdaBoost surro-gate model SRA framework to the problem of de-termining the probability of clusters of nodes be-coming disconnected after an event.The main contribution of this work is as fol-lows. By decomposing the network into denselyconnected clusters and training a surrogate model todetermine inter-cluster connectivity, we show thatit is possible to quickly estimate the probability ofcluster disconnection using MCS.The remainder of the article is organized as fol-lows. In Section 2, the AdaBoost algorithm and anefficient sampling technique to generate balancedtraining data are reviewed. In Section 3, the clus-ter connectivity surrogate model is presented. Sec-tion 4 demonstrates the application of the devel-oped methods in a numerical example of the Cal-ifornia gas distribution network. In Section 5 con-clusions and future work are discussed.2. BACKGROUND ON ADABOOSTSURROGATE MODELS FOR SRAIn this section, the main techniques used in thesource–terminal network connectivity problem (Raiand Aggarwal, 1978) are reviewed. The generalmethodology of constructing a surrogate model us-ing AdaBoost and training the model on a labeleddataset generated via guided random walk samplingare discussed. These tools will also be leveraged inthe cluster connectivity problem used in this work.212th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 20152.1. Estimating source–terminal disconnectionsThe infrastructure network is modeled as a graphG = (V ,E ), where the number of nodes n = |V |is the cardinality of the vertex (node) set V , andE ⊆{(i, j) : i, j ∈ V } is the set of edges. A networkstate is defined as x ∈ {0,1}n, where each elementxi encodes the state of the corresponding networkcomponent asxi ={0 if node i has failed,1 if node i remains intact.(1)Similarly, the network status ys,t with respect tosource node s and terminal node t is defined as fol-lows:ys,t ={1 if s and t are connected,−1 otherwise.(2)Suppose, following an event E, the failure prob-ability Pr(xi = 0|E) of each component is known.In order to estimate network failure probabilityPr(ys,t = −1|E), we generate network realizationsx(k), where k denotes the sample index, and com-pute ys,t(k) = f (x(k)) for each sample k. The oper-ator f that determines the connectivity of the sourceand terminal node is a shortest path algorithm. Ifa finite cost path exists between s and t, then thenodes are connected.Given K realizations of the network state, the net-work failure probability is approximated as:Pr(ys,t =−1|E)≈1KK∑k=1I ( f (x(k)) =−1) (3)where I(·) is the indicator function that takes thevalue 1 if the event occurs, and 0 otherwise. Eval-uating f for large K can be computationally expen-sive. To improve the computational performance,Stern et al. (2014) proposed the use of a fast but ap-proximate surrogate model f˜ , which approximatesthe value of ys,t as y˜s,t = f˜ (x).2.2. Generating data to train the surrogate modelWhen considering the binary status of componentsin the network, there exist a factorial number ofpossible network states (combination of failed andnon-failed nodes). Therefore, it is necessary tosample a subset of these data points to train and testthe surrogate model f˜ . Furthermore, due to the in-herent robustness of infrastructure networks, thereis a bias toward connected networks when samplinglikely network states. For best classification per-formance, it is important to have a dataset that isbalanced between examples of failed and not failednetworks (Japkowicz and Stephen, 2002).In order to generate a balanced training dataset,data points are sampled according to a guidedrandom walk sampling (GRWS) method (Sternet al., 2014). The sampling method, based on theMetropolis-Hastings algorithm (Metropolis et al.,1953), seeks to identify the boundary between net-work failure and non-failure cases on the basis ofthe number of nodes that have failed. A completedescription of the guided random walk samplingtechnique can be found in Stern et al. (2014).The network status ys,t must be determined foreach data point in the training data using the exactmodel f , which is still computationally expensive.The key benefit is this training data needs to be gen-erated once offline, and future queries can use thesurrogate model after it has been trained.2.3. AdaBoost classificationThe surrogate model for source–terminal connec-tivity is constructed using the AdaBoost machinelearning classifier. AdaBoost is chosen because ithas been shown to learn complex structure withlimited training data (Schapire et al., 1998). Ad-aBoost is a binary classification algorithm that com-bines multiple weighted weak classifiers that, inaggregate, produce a more sophisticated classifier.This general approach is called boosting in the ma-chine learning community (Freund and Schapire,1999). The final classifier is given as:f˜ (x) = sign(T∑t=1αtht(x)), (4)which is a linear combination of T weak classifiersht with weights αt .Each weak classifier ht is a one dimensional clas-sification rule. AdaBoost sequentially computeseach weak classifier, indexed by t, and determinesits contribution αt to the overall classifier according312th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015to:αt =12ln(1− εtεt)> 0, (5)where εt is the classification error of the t th weakclassifier, defined subsequently.The weak classifier ht is determined by minimiz-ing the classification error of the m training datapoints {(x(1),y(1)), . . . ,(x(m),y(m))} as:εt =1mm∑k=1[Wt(k) · (ht(x(k)) 6= y(k))]. (6)To compute the classification error, the datapoints are weighted according to the weights{Wt(1), · · · ,Wt(m)}.For the first classifier, these weights are initial-ized as W1(k) = 1/m. For subsequent classifiers,the new weight for data point x(k) is computedfrom the previous data weight Wt−1(k), the previousclassifier ht−1 operating on x(k), and the classifierweight αt−1 as follows:Wt(k) =Wt−1(k)Zt−1exp(−αt−1y(k)ht−1 (x(k))) , (7)where Zt−1 is a normalization constant.The AdaBoost classifier in Stern et al. (2014) istrained on data collected using the guided randomwalk sampling method, and applied to classify anode pair in the network as connected or discon-nected. Cross validation results (not shown) pro-duced similar error on both the training and testdata, indicate the model is not being over-fitted.3. SURROGATE MODELS FOR CLUSTERCONNECTIVITYIn this work, we extend the AdaBoost surrogatemodel for use in an approximate method to esti-mate the probability of disconnection of two clus-ters of nodes. To construct the method, the infras-tructure network is first partitioned into a numberof densely connected clusters of nodes. Once theclusters have been determined, the connectivity ofa cluster–cluster pair can be determined by gener-ating a surrogate model for inter-cluster connectiv-ity. Owing to the good performance on the source–terminal connectivity problem, we again use an Ad-aBoost classifier, where one classifier is trained foreach cluster–cluster pair. The details of this proce-dure are described next.3.1. Graph clusteringWe assume a common structure for large infrastruc-ture networks. Specifically, we assume there area number of densely connected components (e.g,in an urban infrastructure grid), which are looselyconnected to other densely connected components(e.g., in another urban area) via long-range links.These network structures lend themselves well tocluster connectivity analysis. Clustering of infras-tructure networks has also been proposed to under-stand network hierarchy (Gómez et al., 2013) andto facilitate analysis of large networks on multiplescales (Lim et al., 2015).Spectral graph clustering (SGC) (von Luxburg,2007) is used to partition graphs into densely con-nected clusters. To split a graph into two clustersLaplacian L is computed asL = D−A, (8)where D is a diagonal matrix containing the degreeof each node on the diagonal, and A is the adja-cency matrix of the graph. Clusters are determinedby finding the second smallest eigenvalue of L andassigning node labels, qi according toqi ={0 if vi < 0,1 otherwise,(9)where vi is the ith element of the eigenvector as-sociated with the second smallest eigenvalue of L.A cluster is the set of all nodes that have the samenode label qi.For an arbitrary number of clusters, k, the first keigenvectors of L are clustered using k-means clus-tering and used to determine the cluster assignment(von Luxburg, 2007).3.2. Training the surrogate modelSimilar to the methodology in Stern et al. (2014),we generate training data using the guide randomwalk sampling method, with the distinction thatconnectivity is defined with respect to cluster i andj as opposed to the source and terminal nodes in(2). Since each surrogate model only describes412th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015Figure 1: California gas distribution network with fourclusters.the connectivity of two specific clusters within thenetwork, we construct(nc2)models describing eachpossible cluster-cluster pair, where nc is the totalnumber of clusters.Once the surrogate models are built, the modelenters the application phase. When an event oc-curs, fragility models are used to estimate the fail-ure probability of each network component, andnetwork state samples are drawn from these compo-nent failure distributions. Each sample is then runthrough the AdaBoost classifier to estimate if theclusters in question are disconnected for the givennetwork realization. Finally, the samples are usedcollectively to generate an estimate on the proba-bility of cluster disconnection given the event.4. NUMERICAL EXAMPLE4.1. California gas distribution networkWe use the California gas distribution network asa numerical example to test the proposed method.The network topology, obtained from Lim et al.(2015), is shown in Figure 1. The network consistsof 244 components (70 nodes and 87 bi-directionaledges). Each node represents a substation in thegas distribution network, and the edges representgas pipelines.We consider a case where all nodes in the net-work have the same uncorrelated failure probabil-ity. Without loss of generality, only nodes are al-lowed to fail in this network. To relax this assump-tion, one can place a node at the midpoint of eachedge and assign to it the corresponding edge failureprobability. We divide the network into four nodeclusters and consider all possible cluster combina-tions. A surrogate model is trained for each pair ofnode clusters in the network. Each model is trainedusing AdaBoost, and built using m = 5,000 trainingpoints and a maximum of T = 200 weak classifiers.We evaluate the trained models by estimating theprobability that the two clusters become discon-nected using the corresponding AdaBoost classi-fier. All nodes are assumed to have a uniform fail-ure probability p f = 0.15, since it was found toyield the results with the highest error in Stern et al.(2014). A uniform component failure probability isassumed for simplicity, but this can be generalizedto non-uniform component failure without modi-fying the methodology, though performance mayvary. The inter-cluster disconnection probability isestimated via MCS with K = 5,000 samples. Wecompare this to the estimate computed using Dijk-stra’s shortest path algorithm to check connectivitybetween clusters (Dijkstra, 1956).In addition to the failure probability estimate, wealso compute the coefficient of variation (COV) onthe true probability of network failure based on theestimate, δ pˆ f , computed by:δ pˆ f =s√K pˆ f(10)where K is the total number of MCS, s is the samplestandard deviation of network failure probabilityestimates, and pˆ f is the estimated network failureprobability. The time required to estimate the net-work failure probability using the AdaBoost classi-fier is recorded and compared with the time to esti-mate the same probability using Dijkstra’s shortestpath method. The algorithms are implemented inMatlab on a quad core 2.6 GHz MacBook Pro, andare available online at Stern (2015).4.2. ResultsThe node clusters for the California gas distribu-tion network are shown in Figure 1. This choice ofclusters divides the network into four densely con-nected regions of approximately equal size.512th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015(a) Disconnection probability estimate ofclusters 1 and 2 convergence, and COV.(b) Disconnection probability estimate ofclusters 1 and 3 convergence, and COV.(c) Disconnection probability estimate ofclusters 1 and 4 convergence, and COV.(d) Disconnection probability estimate ofclusters 2 and 3 convergence, and COV.(e) Disconnection probability estimate ofclusters 2 and 4 convergence, and COV.(f) Disconnection probability estimate ofclusters 3 and 4 convergence, and COV.Figure 2: Convergence of disconnection probability estimate of cluster pairs.612th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015The results for each cluster pair are summarizedin Table 1, which gives the estimation error, as wellas the computational time for the MCS estimate us-ing the AdaBoost surrogate model (TML), and usingthe shortest path method (TSP). The error is com-puted as the difference between the MCS estimatesusing the exact shortest path method and the sur-rogate model. Overall, the error introduced due tothe surrogate model is less than five percent acrossthe various clusters, while the computation time isapproximately two orders of magnitude faster. It isnoted that the surrogate model is observed to under-predict failure in almost all cases. However, moretesting is needed to identify if this is a general re-sult or specific to the experiments presented in thisarticle.An interesting result is the perfect performancefor predicting connectivity between clusters 2 and3. This is due in part to the network structure,since there is only one edge between the two clus-ters. As a result, a model is learned that only re-quires one weak classifier for accurate classifica-tion. Thus, TML is reduced compared to the othersurrogate models. Figure 2d shows the convergenceof the the estimated probability of disconnection forboth MCS approaches, and is an example of thebest-case performance of the surrogate model.Figures 2a through 2f show the general trend thatclusters that are further apart have a higher prob-ability of disconnection than those that are adja-cent. Furthermore, the clusters that are further apartexhibit a greater difference between the surrogatemodel estimate and shortest path estimate. For ex-ample, Figure 2c shows that the two estimates haveconverged on slightly different values. The greatestconvergence gains occur in the first 1,000 simula-tions.5. CONCLUSIONS AND FUTURE WORKThis work demonstrates how to use machine learn-ing methods to construct a surrogate model for esti-mating cluster connectivity via MCS. The approachimproves the efficiency significantly and introducesrelatively small (less than five percent) errors to thefailure probability estimate.Several extensions to this work are currently un-der exploration. First, the influence of correlatedCluster Error (%) TML (s) TSP (s)1 - 2 0.72 110.6 1.480×1041 - 3 2.82 109.3 2.060×1041 - 4 4.52 112.2 1.504×1042 - 3 0.00 1.691 3.727×1042 - 4 0.82 109.4 2.067×1043 - 4 0.74 112.9 2.074×104Table 1: Summary of results for all six cluster-to-cluster models with component failure probabilityp f = 0.15. Computation time for machine learningmethod TML and shortest path method TSP provided.component failures on the surrogate model accu-racy are being investigated, since this may occurin many practical applications. Second, we are alsointerested in developing a machine-learning basedregression approach to predict the probability ofnetwork failure in terms of network flow quanti-ties. Furthermore, we are exploring the possibil-ity of combining all cluster-cluster surrogate mod-els into one super-model that would incorporate allthe individual surrogate models.6. ACKNOWLEDGMENTSThe authors would like to thank the US NationalScience Foundation for funding under grant num-ber CMMI 1031318. Any opinions, findings andconclusions or recommendations expressed in thismaterial are those of the authors and do not nec-essarily reflect the views of the National ScienceFoundation. The second author, Junho Song wouldlike to acknowledge support by a grant entitledDevelopment of cutting edge technologies for themulti-faceted representation of design earthquakeground motions based on analyses of accelera-tion records [NEMA-NH-2013-71] provided by theNatural Hazard Mitigation Research Group, Na-tional Emergency Management Agency of Korea.7. REFERENCESBoin, A. and McConnell, A. (2007). “Preparing for crit-ical infrastructure breakdowns: The limits of crisismanagement and the need for resilience.” Journal ofContingencies and Crisis Management, 15(1), 50–59.Bruneau, M., Chang, S. E., Eguchi, R. T., Lee, G. C.,O’Rourke, T. D., Reinhorn, A., Shinozuka, M., Tier-712th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015ney, K., Wallace, W. A., and von Winterfeldt, D.(2003). “A framework to quantitatively assess and en-hance the seismic resilience of communities.” Earth-quake Spectra, 19(4), 733–752.Der Kiureghian, A. and Song, J. (2008). “Multi-scale re-liability analysis and updating of complex systems byuse of linear programming.” Reliability Engineering& System Safety, 93(2), 288–297.Dijkstra, E. (1956). “A note on two problems in con-nexion with graphs.” Numerische Mathematik, 1(1),269–271.Ditlevsen, O. and Madsen, H. (1996). Structural Relia-bility Methods. John Wiley & Sons, Chichester, UK.Freund, Y. and Schapire, R. (1995). “A desicion-theoretic generalization of on-line learning and an ap-plication to boosting.” Computational Learning The-ory, P. Vitányi, ed., Vol. 904 of Lecture Notes in Com-puter Science, Springer Berlin Heidelberg, 23–37.Freund, Y. and Schapire, R. E. (1999). “A short intro-duction to boosting.” Journal of Japanese Society forArtificial Intelligence, 14(5), 771–780.Gómez, C., Sanchez-Silva, M., Dueñas-Osorio, L.,and Rosowsky, D. (2013). “Hierarchical infrastruc-ture network representation methods for risk-baseddecision-making.” Structure and Infrastructure Engi-neering, 9(3), 260–274.Jansson, T., Nilsson, L., and Redhe, M. (2003). “Us-ing surrogate models and response surfaces in struc-tural optimization – with application to crashworthi-ness design and sheet metal forming.” Structural andMultidisciplinary Optimization, 25(2), 129–140.Japkowicz, N. and Stephen, S. (2002). “The class imbal-ance problem: A systematic study.” Intelligent DataAnalysis, 6(5), 429–450.Lim, H.-W. and Song, J. (2012). “Efficient risk assess-ment of lifeline networks under spatially correlatedground motions using selective recursive decompo-sition algorithm.” Earthquake Engineering & Struc-tural Dynamics, 41(13), 1861–1882.Lim, H.-W., Song, J., and Nolan, K. (2015). “Seis-mic reliability assessment of lifeline networks usingcluster-based multi-scale approach.” Earthquake En-gineering Structural Dynamics, 44(3), 355–369.Meirlaen, J., Huyghebaert, B., Sforzi, F., Benedetti, L.,and Vanrolleghem, P. (2001). “Fast, simultaneoussimulation of the integrated urban wastewater systemusing mechanistic surrogate models.” Water Scienceand Technology, 43(7), 301–307.Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N.,Teller, A. H., and Teller, E. (1953). “Equation of statecalculations by fast computing machines.” The Jour-nal of Chemical Physics, 21(6), 1087–1092.Papadrakakis, M., Papadopoulos, V., and Lagaros, N. D.(1996). “Structural reliability analysis of elastic-plastic structures using neural networks and montecarlo simulations.” Computer methods in applied me-chanics and engineering, 136(1-2), 145–163.Rai, S. and Aggarwal, K. K. (1978). “An efficientmethod for reliability evaluation of a general net-work.” IEEE Transactions on Reliability, 27(3), 206– 211.Rausand, M. and Høyland, A. (2004). System Reliabil-ity Theory: Models, Statistical Methods, and Appli-cations. John Wiley & Sons, Chichester, UK.Reed, D., Kapur, K., and Christie, R. (2009). “Method-ology for assessing the resilience of networked infras-tructure.” Systems Journal, IEEE, 3(2), 174–180.Schapire, R. E., Freund, Y., Bartlett, P., and Lee, W. S.(1998). “Boosting the margin: a new explanation forthe effectiveness of voting methods.” The Annals ofStatistics, 26(5), 1651–1686.Song, J. and Der Kiureghian, A. (2003). “Bounds onsystem reliability by linear programming.” Journal ofEngineering Mechanics, 129(6), 627–637.Stern, R. (2015). “California gasdistribution network source code,.Stern, R. E., Song, J., and Work, D. B. (2014). “Machinelearning-based surrogate models for network reliabil-ity analysis.” 17th IFIP WG 7.5 July 3-7.von Luxburg, U. (2007). “A tutorial on spectral cluster-ing.” Statistics and Computing, 17(4).Vugrin, E., Warren, D., Ehlen, M., and Camphouse, R.(2010). “A framework for assessing the resilienceof infrastructure and economic systems.” Sustain-able and Resilient Critical Infrastructure Systems, K.Gopalakrishnan and S. Peeta, eds., Springer BerlinHeidelberg, 77–116.Wan, X., Pekny, J. F., and Reklaitis, G. V. (2005).“Simulation-based optimization with surrogate mod-els—application to supply chain management.” Com-puters & Chemical Engineering, 29(6), 1317–1328.8