

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK By JULIE ANNA BRIGHT Bachelor of Science in Industrial Engineering and Management Oklahoma State University Stillwater, OK 74075 2009 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE December 2011 COPYRIGHT c By JULIE ANNA BRIGHT December 2011 ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK Thesis Approved: Balabhaskar Balasundaram Thesis Advisor Manjunath Kamath Ricki Ingalls Sheryl Tucker Dean of the Graduate College iii TABLE OF CONTENTS Chapter Page 1 Introduction 1 1.1 Organization of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 The Deterministic Shortest Path Problem . . . . . . . . . . . . . . . . 2 1.2.1 Formal Definition . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Classical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.3 Linear Programming Formulation . . . . . . . . . . . . . . . . 4 1.3 Stochastic Shortest Path Problems . . . . . . . . . . . . . . . . . . . . . 5 1.4 Research Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 Research Tasks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Review of Modeling Approaches 8 2.1 Stochastic Programming Models . . . . . . . . . . . . . . . . . . . . . 8 2.2 Robust Optimization Models . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Other Modeling Approaches . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Models Based on Conditional ValueatRisk . . . . . . . . . . . . . . . 11 2.5 Formal Definition of CVaR . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.6 Using CVaR in Optimization . . . . . . . . . . . . . . . . . . . . . . . . 15 3 Model Development and Analysis 17 3.1 Modeling Arc Failures . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Formulation of the CVARSPP . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.1 Additional Parameters . . . . . . . . . . . . . . . . . . . . . . . 18 iv 3.2.2 CVARSPP Optimization Model . . . . . . . . . . . . . . . . . . 18 3.3 Loss Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.1 Path Reliability . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.2 Number of Arc Failures . . . . . . . . . . . . . . . . . . . . . . 22 3.3.3 Number of Detours . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 Comparison of Loss Functions . . . . . . . . . . . . . . . . . . . . . . 29 3.4.1 Magnitude of Losses . . . . . . . . . . . . . . . . . . . . . . . . 29 3.4.2 Dependence of CVaR on b . . . . . . . . . . . . . . . . . . . . . 30 4 Computational Experiments 31 4.1 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2.1 Group 1: Complete Enumeration . . . . . . . . . . . . . . . . . 33 4.2.2 Group 2: Branch Networks . . . . . . . . . . . . . . . . . . . . 34 4.2.3 Group 3: Uniform Random Networks . . . . . . . . . . . . . . 35 4.3 Comparison of Loss Functions . . . . . . . . . . . . . . . . . . . . . . 36 5 Conclusion and Future Work 38 Bibliography 39 A Optimization Models 44 A.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 A.2 Decision Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 A.3 Path Reliability Loss Function . . . . . . . . . . . . . . . . . . . . . . . 45 A.4 Arc Failure Loss Function . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.5 Detours Loss Function . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 B Database Schema 48 B.1 Database Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 v B.2 Database Schema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 C Objective Function Graphs 55 D Test Network Diagrams 70 D.1 Small Network: exp0002 . . . . . . . . . . . . . . . . . . . . . . . . . . 70 D.2 Branch Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 D.3 Uniform Random Networks . . . . . . . . . . . . . . . . . . . . . . . . 72 vi LIST OF TABLES Table Page 3.1 Summary of Loss Functions (LF) . . . . . . . . . . . . . . . . . . . . . 29 4.1 Summary of Test Bed . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 vii LIST OF FIGURES Figure Page 2.1 An example of VaR and CVaR . . . . . . . . . . . . . . . . . . . . . . . 12 3.1 Path Reliability Loss Function: The loss for this path is 1. . . . . . . . 21 3.2 Arc Failures Loss Function: The loss for this path is 3. . . . . . . . . . 22 3.3 Detours Loss Function: The loss for this path is 2. . . . . . . . . . . . 24 4.1 Network exp0002 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Branch Network Structure . . . . . . . . . . . . . . . . . . . . . . . . . 35 B.1 Legend: Crow’s Foot Notation . . . . . . . . . . . . . . . . . . . . . . 53 B.2 Database Schema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 C.1 exp0002: Level sets of the optimal path length . . . . . . . . . . . . . 56 C.2 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 57 C.3 branch02: Level sets of the optimal path length . . . . . . . . . . . . . 58 C.4 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 59 C.5 branch03: Level sets of the optimal path length . . . . . . . . . . . . . 60 C.6 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 61 C.7 branch04: Level sets of the optimal path length . . . . . . . . . . . . . 62 C.8 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 63 C.9 gnp01: Level sets of the optimal path length . . . . . . . . . . . . . . . 64 C.10 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 65 C.11 gnp02: Level sets of the optimal path length . . . . . . . . . . . . . . . 66 C.12 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 67 viii C.13 gnp03: Level sets of the optimal path length . . . . . . . . . . . . . . . 68 C.14 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 69 D.1 Network exp0002 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 D.2 Network branch02 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 D.3 Network branch03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 D.4 Network branch04 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 D.5 Network gnp01 — GNP(25, 0.1) . . . . . . . . . . . . . . . . . . . . . . 72 D.6 Network gnp02 — GNP(25, 0.25) . . . . . . . . . . . . . . . . . . . . . 73 D.7 Network gnp03 — GNP(25, 0.4) . . . . . . . . . . . . . . . . . . . . . . 73 ix CHAPTER 1 Introduction to Shortest Path Problems Under Uncertainty The shortest path problem (SPP) is to find the path through a network with the lowest “cost” — whether that cost is measured in miles, money, minutes, or any other metric. This problem has been extensively studied over the past 50 years (Pollack and Wiebenson, 1960; Dreyfus, 1969; Deo and Pang, 1984). When the structure and costs of the network are known, then very fast and powerful tools exist to find a shortest path, such as Dijkstra’s algorithm (Dijkstra, 1959). In the real world, though, networks are afflicted with failures, and costs can fluctuate randomly. This limits the applicability of classical, deterministic SPP algorithms. Often the definition of a “shortest” path needs to be adjusted. Different mathematical models and algorithms are needed to handle uncertainty. The family of shortest path problems under uncertainty (modeled by probability distributions) are collectively called stochastic shortest path problems (SSPPs). This thesis investigates one type of SSPP. In this SSPP, the arcs of a network can fail randomly and independently with some probability. This model is useful for roads that close due to accidents or severe weather or for PERT charts where activities may not get done. The effect of arc failures on the overall network and the shortest path will be investigated in order to create measures to quantify the loss of an arc. Secondly, a quantitative measure of downside risk called Conditional ValueatRisk (CVaR) will be used to find short and robust paths in the network (Rockafellar and Uryasev, 2002). 1 1.1 Organization of this Thesis Chapter 1 discusses the classical or deterministic shortest path problem (SPP), including a formal definition and a review of the most common algorithms, before providing some motivating examples to show how uncertainty in graphs affects the SPP. This chapter also introduces the research problem. Chapter 2 reviews the modeling approaches in the literature to show what kinds of uncertainty have been captured and what kinds of uncertainty remain relatively lessstudied. This chapter also introduces the concept of CVaR. Chapter 3 formulates the CVaRConstrained Stochastic Shortest Path Problem (CVARSPP) and develops three loss functions to model loss due to probabilistic arc failures to use in the CVARSPP. This chapter also analyzes the loss functions, deriving their distributions for uniform arc failures and some implications for their behavior in CVARSPP models. Then, Chapter 4 presents some exploratory computational experiments that illustrate the behavior of the CVARSPP models with the loss functions, leading to some preliminary guidance on setting parameters. Finally, Chapter 5 presents the conclusions of this thesis and some directions for further research. 1.2 The Deterministic Shortest Path Problem 1.2.1 Formal Definition Consider the network G = (N, A), consisting of N, a set of nodes, and A N N, a set of directed arcs connecting pairs of nodes. For each arc (ij) 2 A, node i is called the tail and node j is called the head. Each arc also has an associated cost cij 2 R+. For each node i 2 N, its inneighborhood is defined as N(i) = fjj(ji) 2 Ag and its outneighborhood is defined as N+(i) = fjj(ij) 2 Ag. Adirected path P is a sequence of k distinct nodes (i1, i2, . . . , ik), where (ij, ij+1) 2 2 A for all j 2 f1, . . . , k1g. Throughout the remainder of this thesis, directed paths will be referred to simply as “paths”. The edgeinduced subgraph of the path P is denoted by G(P) = (N(P), A(P)). The path P can also be represented by an incidence vector x 2 f0, 1gjAj of arcs in P, which induces the subgraph G(x) = (N(x), A(x)). The cost, also called weight or length, of a path P is cP = å (ij)2A(P) cij. Definition 1.2.1. (Deterministic Shortest Path Problem) Given network G = (N, A) and nodes s, t 2 N, the PointtoPoint Shortest Path Problem is to find a minimum cost path from node s to t in G. The SingleSource Shortest Path Problem is to find a shortest path from a source node s 2 N to each node i 2 Nnfsg. The AllPairs Shortest Path Problem is to find a shortest path from every node i 2 N to every other node j 2 N. 1.2.2 Classical Results for Shortest Path Problem In the absence of negative cost cycles in the network, the SPP is a wellsolved problem, for which several stronglypolynomial algorithms are known (Deo and Pang, 1984). These algorithms assume that the weight of each arc in the network is constant, and network structure (nodes and arcs) are known with certainty. Algorithms for the SPP fall into two main categories: labelsetting and labelcorrecting (Deo and Pang, 1984). Dijkstra’s algorithm (Dijkstra, 1959) is an example of a labelsetting algorithm. The BellmanFord algorithm (Bellman, 1958) is a labelcorrecting algorithm derived from dynamic programming. The labelsetting algorithms require all arc costs to be nonnegative, while the labelcorrecting algorithms only require that the network have no negative cost cycles. The Floyd Warshall algorithm for the AllPairs Shortest Path Problem (Floyd, 1962) is also based on dynamic programming. The SingleSource and PointtoPoint Shortest Path Problems can also be formulated as a linear program and solved with the net 3 work simplex method when the network contains no negative cost cycles (Dantzig, 1957). 1.2.3 Linear Programming Formulation of the Shortest Path Problem The linear programming (LP) formulation of PointtoPoint SPP optimization model is worth presenting on its own since this thesis extends this formulation. The constraint matrix for the LP formulation of the SPP is totally unimodular, so the extreme point optimal solutions are integral (Ahuja et al., 1993). Inputs Network G = (N, A) is a directed network. s 2 N is the source node. t 2 N is the sink node. cij is the cost of arc (ij) 2 A. Decision Variables xij for all (ij) 2 A. At any x 2 f0, 1gjAj, xij = 1 implies that arc (ij) is in the solution and xij = 0 implies that arc (ij) is not. Formulation 1.2.2 (Deterministic Formulation). Minimize å (ij)2A cijxij (1.1) 4 subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (1.2) xij 0 8 (ij) 2 A (1.3) 1.3 Stochastic Shortest Path Problems Uncertainty in the system underlying a network model may be modeled by random variation in network parameters such as arc costs and in the structure of the network, particularly node and arc failures. Once uncertainty is allowed into the network model, there is no longer one welldefined shortest path problem. A family of Stochastic Shortest Path Problems (SSPPs) arise instead, depending on which elements in the network are uncertain and how the uncertainty is modeled. Furthermore, SSPPs generally require the redefinition of the concept of feasible and optimal solutions. One SSPP models the arc cost cij as a random variable. Applications of this model include modeling the travel time for an arc in a congested road network or the task duration in a PERT network (Ball et al., 1995a). For this problem, two ways to redefine the shortest path are in terms of expected length, which gives rise to many stochastic programming approaches, and worstcase length, in which arc costs are set at or close to their maximum possible values, which lies behind many robust optimization approaches (Ball et al., 1995b). Uncertainty about which node is the source or sink can be modeled by randomly assigning a supply of +1 to one node or 1 to another. This would be useful in modeling ambulance routing during an emergency when the driver dis 5 covers which ER has room for the patient enroute. Allowing uncertainty in the structure of the network (e.g. uncertainty in the set of nodes or arcs, or only partial knowledge) gives rise to another set of SSPPs. One source of uncertainty is limited knowledge of the network, such as not knowing the cost or head node of an arc until the tail is reached. This is frequently used to model hideandseek or adversarial search on networks (Papadimitriou, 1991) or for navigating real road networks without a map, such as after a major disaster (Papadimitriou, 1991). Secondly, the nodes may be probabilistic—that is, they can fail. This is used in modeling telecommunications networks or computer networks (Jaillet, 1992). Finally, the arcs may fail with some probability, giving rise to the subject of this thesis. Allowing arcs to fail allows modeling a wide variety of disturbances in networks. For instance, road closure in transportation networks, due to accidents, weather, or other events, can be modeled as an arc failure. Cutting through the cables in fiber optic networks, which happens during construction or fires (Ball et al., 1995a), can also be modeled as arc failure. In a social network, an onagain offagain friendship between two people can be modeled as an arc that fails after an argument and is restored after a reconciliation. It is therefore desirable to find a short, reliable path through a network. The redefinition of the shortest path depends on the trade off between reliability and shortness. 1.4 Research Problem This thesis investigates models of the Stochastic Shortest Path Problem with Probabilistic Arc Failures (SSPPPAF) using Conditional ValueatRisk (CVaR) constraints. CVaR is a risk management tool that can trade off risk aversion and performance, so it is able to handle varying levels of risk tolerance. Moreover, CVaR is based on an explicit loss function, which allows more precision in modeling the effects 6 of arc failures. The concept of CVaR is discussed in more detail in Chapter 2. This thesis will limit itself to arcs that fail randomly and independently of each other with some probability. All other graph parameters are deterministic. 1.5 Research Tasks The objective of this thesis is to develop CVaRbased optimization models for the SSPPPAF. The solutions to these models should be paths that are short and robust under arc failures; specifically, the loss due to arc failure as measured by using CVaR should be low. To achieve this objective, the tasks listed below will be pursued. 1. Model the effects of arc failures using different loss functions. 2. Develop and implement optimization models to find the shortest path subject to CVaR constraints. 3. Investigate the effect of different loss functions and CVaR bound on the feasibility of the optimization model. 4. Gain insight into the difference between models using expected value and CVaR via their behavior with different loss functions on carefullychosen test instances. 7 CHAPTER 2 A Review of Modeling Approaches for the Shortest Path Problem under Uncertainty 2.1 Stochastic Programming Models Stochastic programming is one of the approaches available to model optimization problems under uncertainty. Stochastic programming models typically use expected value to model the effect of uncertainty in the objective function and to allow for recourse actions (Sen and Higle, 1999). As a measure of risk, expected value is risk neutral. It weights gains as high as losses whereas human decisionmakers are often much more riskaverse (Canada et al., 2004), preferring to minimize losses rather than maximize gains. There have been a few stochastic programming models of the SSPPPAF. Croucher (1978) incorporated recourse decisions into a labelcorrecting algorithm for directed acyclic networks using dynamic programming. The paper attempts to model a traveler who has imperfect knowledge of a road network. In particular, the traveler does not know whether an arc is blocked until arriving at its head. The objective of the algorithm is to find a path with the shortest expected length from each node to the sink node given that the desired arc is unreliable. In the algorithm, each node is labeled with the expected shortest path length to the sink. Each arc in the network has a probability pij of being usable. At each node an outgoing arc is selected. If the desired arc is broken, one of the other outgoing arcs of the node is selected at random. When the node has only one outgoing arc, it is not allowed to fail. Each node is labeled with the expected shortest path length to the sink, which 8 is the average of the labels of the outneighbors weighted by their probability of transversal. This algorithm is quick to compute, but it relies on the assumption that at least one outgoing arc at each node will be usable. Andreatta and Romeo (1988) built on Croucher (1978)’s paper to explicitly find optimal recourse for a blocked arc, and they find a path that minimizes the expected length of the path including recourse, which they define as detours around a failed arc. The model uses stochastic dynamic programming to solve the problem. This model suffers from the curse of dimensionality as the number of scenarios and the number of failed arcs per scenario grow. Like the previous paper, this one seeks only to minimize the expected path length and does not give special consideration to the extremes of the distribution. 2.2 Robust Optimization Models Robust optimization is another approach to optimization under uncertainty. Robust optimization models seek to find solutions that are feasible no matter how the underlying uncertainty resolves (Bertsimas et al., 2008). The uncertainty is captured in an uncertainty set, which is a collection of possible future states of the system being modeled. The optimal solution must be feasible for every element in the uncertainty set. The uncertainty set need not contain every possible resolution of the uncertainty, only the subset for which the the optimal solution must remain feasible. Robust optimization tends to favor minimax objectives (Bertsimas et al., 2008). Recently Yu and Yang (1998) showed a version of the robust shortest path problem to be NPcomplete. The paper considers a network with uncertain arc costs, which are captured in a set of scenarios, and develops two robust optimization models for the problem. The first model uses a minimax objective over the scenarios. The second model seeks a path with the smallest possible range of lengths. 9 Although the authors develop a dynamic programming algorithm for networks with bounded numbers of scenarios, the algorithm is pseudopolynomial in the number of scenarios. As the number of scenarios needed to capture the behavior of a large graph grows very quickly, this algorithm is not suitable for large scale uses (Yu and Yang, 1998). 2.3 Other Modeling Approaches A conceptually similar approach to the SSSPPAF is the Minimum Cost Reliability Ratio Path Problem (MCRRPP), a bicriteria optimization model proposed by Ahuja (1988). This problem considers paths from a single source node to a single sink in directed network with arc failures. The objective is to find a path with the optimal ratio of cost to reliability, defined as the probability of no arc failures on the path. The author first develops optimality criteria for a bicriteria optimization problem with the criteria of cost and reliability. The paper develops a pseudopolynomial dynamic programming algorithm with a complexity of O(mnDlogm) to solve the problem, although the computational experiments found that averagecase complexity was much lower. The MCRRPP trades off cost and reliability in a riskneutral manner. By using reliability in the formulation of the MCRRPP, the model implicitly a loss of 100% when any arc on the path fails, and 0% otherwise. This approach is similar to the approach used in this thesis in that it uses arc failures and explicitly trades off cost and reliability. In a related approach, Jaillet (1992) incorporates node failures into a SSPP model. The goal of this model is to find an easily repairable path of minimum expected length given node failures. The paper shows that the SSSP with probabilistic node failures is in general NPcomplete although some restricted types of the problem are solvable in polynomial time. 10 2.4 Models Based on Conditional ValueatRisk CVaR was developed in the late 1990s as an improvement on an earlier financial risk measure called ValueatRisk or VaR (Acerbi and Tasche, 2002). VaR has been very popular, even being written into regulations in several countries (Rockafellar and Uryasev, 2000), but it has a number of drawbacks (Rockafellar and Uryasev, 2000). For most probability distributions, VaR is difficult to compute and is not a coherent risk measure (Artzner et al., 1999). CVaR, on the other hand, is a coherent measure of risk and is easier to calculate (Acerbi and Tasche, 2002), where it is called expected shortfall. More pertinently, Rockafellar and Uryasev presented a way to use CVaR in optimization problems via a portfolio optimization example (Rockafellar and Uryasev, 2002). Intuitively, CVaR is the average of the upper tail of the distribution of losses due to uncertainty, which are quantified in a loss function. The loss function and the threshold for the upper tail are provided by the modeler. Figure 2.1 illustrates VaR and CVaR for a negative binomial distribution with the tail threshold set at the 90th percentile. Currently CVaR is applied in several areas of finance such as insurance and credit risk evaluation (Rockafellar and Uryasev, 2000). In these applications, risk is usually specified in terms of monetary losses (Rockafellar and Uryasev, 2000) although utility can be used as well. In more general contexts, though, CVaR is defined in terms of a realvalued loss function L(x,Y), where Y is a vector of random variables that model uncertainties in the underlying system, and x is a vector of decisions. In portfolio optimization, Y would be future share prices for all companies traded in the stock market and x the portfolio allocation. For each decision x, the loss function L(x,Y) is itself a random variable. Generally gains are modeled as negative losses. In portfolio optimization, the loss is often measured in terms of 11 Figure 2.1: An example of VaR and CVaR 0 5 10 15 20 25 30 35 40 45 50 55 60 0 1 102 2 102 3 102 4 102 5 102 6 102 CVaR pmf 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 VaR cd f return on investment (ROI), where the loss is negative if the ROI is positive. A recent paper (Madadi et al., 2010) formulates a CVaRbased supply chain network model in order to minimize the risk of defective pharmaceuticals in a medical supply chain. Supplier failures and resulting defective pharmaceuticals are the sources of uncertainty in this network. The demand at the customer nodes is fixed. In the transportation network the arcs and their costs are deterministic. The objective function minimizes the CVaR of the loss from defective pharmaceuticals. The constraints trade off the cost of inspection versus the cost of defective products. They enumerate all possible scenarios and assign them a probability. This was computationally feasible because the number of suppliers is likely to be small compared to the number of customers. One strength of their model is that it allows for a state of partial failure. This paper uses CVaR for managing productionside risks. In contrast to this thesis, modeling uncertainty in the arcs of the underlying 12 transportation network is not within the scope of Madadi et al. (2010)’s paper. 2.5 Formal Definition of CVaR Definition 2.5.1 (Rockafellar and Uryasev (2002)). Let x 2 X Rn be a vector of decision variables and Y : W ! Rm be a vector of random variables defined on the sample space W. Then the loss function L(x,Y) is a function of the decisions x and the random variables Y. The random variable L(x,Y) has a family of distributions, one for each x 2 X. The sample space of L(x,Y) is denoted by WL(x). Given a decision x, denote the probability density function or the probability mass function of L(x,Y) by fL(`; x,Y) where ` 2 WL(x) and the cumulative distribution function by FL(`; x,Y). Note 2.5.2 (Finite Loss Functions). For each loss function discussed in this thesis, EfL(x,Y)g < ¥ for all x 2 X. The assumption that the uncertainties modeled by Y are not affected by the choice of x is important (Rockafellar and Uryasev, 2002). For instance, in a financial setting, this assumption implies that stock prices fluctuate randomly regardless of who owns how much of each stock. Throughout this thesis the distribution of Y, the random vector of arc failures, is assumed to be unaffected by the choice of path x. In addition, it is assumed that arc failures are independent of each other. In particular, for each arc (ij) 2 A, Yij is a Bernoulli random variable, and Y : W ! f0, 1gjAj (see Section 3.1). Modeling of loss functions in the context of the SSPPPAF is discussed extensively in Section 3.3 of this thesis. The VaR(x, b) is the quantile associated with the threshold probability b of the distribution of L(x,Y) and is defined as follows: Definition 2.5.3 (Rockafellar and Uryasev (2002)). VaR(x, b) = min f`jFL(`; x,Y) bg (2.1) 13 Definition 2.5.4 (Rockafellar and Uryasev (2002)). Given a threshold percentile b, the cumulative upper btail distribution, denoted by Fb,L(`; x,Y), is given by: Fb,L(`; x,Y) = 8>>< >>: 0 if ` < VaR(x, b) FL(`;x,Y)b 1b if ` VaR(x, b) (2.2) Definition 2.5.5 (Rockafellar and Uryasev (2002)). For a given threshold percentile b, the CVaR(x, b) of the loss associated with decision x is given by: CVaR(x, b) = E fL(x,Y)jL(x,Y) VaR(x, b)g (2.3) In this thesis the random vector Y represents arc failures, so the distribution of Y is discrete. Hence, the loss functions induced are also discrete distributions. For loss functions with finite discrete distributions, the CVaR of the loss is defined below. Proposition 2.5.6. (Proposition 8, Rockafellar and Uryasev (2002)) Assume that the sample space of Y is finite, so that for each x 2 X the distribution of L(x,Y) is likewise discrete and finite and FL(`; x,Y) is a step function with jumps at ` 2 WL(x). For a fixed x, order ` 2 WL such that `0 < `1 < < `max. Let `k be the unique element of WL such that FL(`k1; x,Y) < b FL(`k; x,Y). (2.4) The VaR(x, b) of the loss is given by VaR(x, b) = `k. (2.5) In addition, the CVaR(x, b) is given by CVaR(x, b) = 1 1 b " (FL(`; x,Y) b) `k + `max å `=`k+1 ` fL(`; x,Y) # . (2.6) 14 2.6 Using CVaR in Optimization The minimization formulas below were developed by Rockafellar and Uryasev (2000, 2002) and rely on the convexity of CVaR. Theorem 2.6.1 (Rockafellar and Uryasev (2002)). Let Fb(x, z) = z + 1 1 b E n [L(x,Y) z]+ o , (2.7) where [z]+ = maxf0, zg. Fb(x, z) is finite and convex as a function of z 2 R and CVaR(x, b) = min z Fb(x, z). (2.8) Denote the set of minimizers of Fb(x, z) by Z. This set is a nonempty, closed interval, reducing to a single point if the minimum z of Fb(x, z) is unique, and VaR(x, b) = minfz 2 Zg. (2.9) In addition, CVaR(x, b) = Fb(VaR(x, b)). (2.10) Since Fb(x, z) is convex, it may be used in optimization problems. Specifically, the next theorem and corollary show that Fb(x, z) can be used to minimize CVaR(x, b) in the objective function or part of the constraints, which is how CVaR is used in the CVARSPP. Theorem 2.6.2 (Rockafellar and Uryasev (2002)). Minimizing CVaR(x, b) over x 2 X is equivalent to minimizing Fb(x, z) over X R since min x2X CVaR(x, b) = min (x,z)2X R Fb(x, z). (2.11) In addition, (x , z ) 2 X R are minimizers of Fb(x, z) (not necessarily unique) if and only if x 2 X is a minimizer of CVaR(x, b) and z is a minimizer of Fb(x , z). 15 Since Fb(x, z) is a convex piecewise linear function of z, it can be used as a constraint or objective function in an optimization problem. This is helpful since the direct formula in Proposition 2.5.6 requires computing VaR(x, b) in advance, which is usually difficult (Rockafellar and Uryasev, 2002). Using Fb(x, z) instead allows the use of CVaR in optimization without having to compute the VaR(x, b) in advance, as the corollary below mentions. The optimal solution found by minimizing CVaR(x, b) via the minimization formula Fb(x, z) will contain information about VaR(x, b) as a byproduct. Corollary 2.6.3 (Rockafellar and Uryasev (2002)). If (x , z ) minimizes Fb(x, z) over X R, then x minimizes CVaR(x, b) over X and CVaR(x , b) = Fb(x , z ) (2.12) and VaR(x , b) z , (2.13) where equality holds if the set of minimizers Z of Fb(x, z) (see Theorem 2.6.1) reduces to a single point. 16 CHAPTER 3 Model Development and Analysis 3.1 Modeling Arc Failures Throughout this thesis the simplifying assumption is made that failures are independent of each other. The probability that arc (ij) 2 A fails is pij where 0 pij < 1. Arcs for which pij = 0 never fail. The random vector of arc failures, Y, is of length jAj. Each component of Y is a random variable Yij. If arc (ij) fails for some ys 2 W, then ys ij = 0. In other words, Yij is a Bernoulli random variable. Yij Bernoulli(1 pij) 8 (ij) 2 A (3.1) 3.2 Formulation of the CVARSPP This model is an extension of the deterministic LP model of the shortest path problem in Section 1.2.3. It simply adds a constraint on CVaR to the LP model. When solved to optimality, a shortest path that satisfies the CVaR constraint will be chosen. Unfortunately, adding a CVaR constraint generally destroys the integrality of the shortest path polytope, so it is necessary to solve the problem as a MILP, greatly increasing solution time. 17 3.2.1 Additional Parameters In addition to the notation used in Formulation 1.2.2, the following notation is used throughout the rest of this thesis. Note that the set of arcs A is random. In the formulation, we denote by A the set of arcs such that pij < 1. C is the maximum acceptable CVaR of loss. b is the percentile for the upper tail (e.g. 0.95, 0.99) S W is the set of scenarios in the model, indexed by s. S may be a random sample from W or the whole sample space if W is small. ps is the probability of scenario s, normalized so that å s2S ps = 1. For equiprobable scenarios, set ps = 1 jSj 8 s 2 S. ys is the realization of the random vector Y corresponding to scenario s 2 S. 3.2.2 CVARSPP Optimization Model Formulation 3.2.1. Minimize å (ij)2A cijxij (3.2) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.3) CVaR(x, b) C (3.4) x 2 f0, 1gjAj (3.5) Constraint 3.4 can be rewritten in terms of the minimization formula Fb(x, z) from Theorem 2.6.1 as in the formulation below. 18 Additional Decision Variables z 2 R is VaR(x, b) at optimality. zs linearizes the excess loss [L(x, ys) z]+ in scenario s 2 S. Formulation 3.2.2. Minimize å (ij)2A cijxij (3.6) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.7) Fb(x, z) C (3.8) x 2 f0, 1gjAj (3.9) z 2 R (3.10) In Constraint 3.8, Fb(x, z) = z + 1 1 b å s2S ps [L(x, ys) z]+ . (3.11) Constraint 3.8 can be linearized using standard techniques and the distribution of L(x,Y) approximated by sampling to produce Formulation 3.2.3. As the loss functions are developed in the next section, the CVARSPP models for each loss function will be based on this formulation. Formulation 3.2.3. Minimize å (ij)2A cijxij (3.12) 19 subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.13) z + 1 1 b å s2S pszs C (3.14) zs L(x, ys) z 8 s 2 S (3.15) zs 0 8 s 2 S (3.16) x 2 f0, 1gjAj (3.17) z 2 R (3.18) 3.3 Loss Functions In this section three different loss function are introduced. Each loss function corresponds to a different effect of arc failures on the path. While the concepts behind these loss functions have appeared in the literature before, explicitly formulating and analyzing them as a part of a CVaR optimization model is the main contribution of this thesis. All the loss functions measure properties of an s t path x. The first loss function is called the Path Reliability Loss Function and measures the reliability of the path. The second function is the Arc Failure Loss Function, which measures the number arc failures. The final loss function is the Detours Loss Function, which counts the number of detours around broken segments of a path. 3.3.1 Path Reliability The first loss function to be considered derives from the basic physical interpretation of the problem. It asks “Does the path break in the current scenario?” In 20 other words, the loss is 0 if no arc fails on the current path and 1 otherwise. This loss function is worth looking at because it is simple and comparable to other SSPP approaches such as Ahuja (1988). This is the loss function implicit in the Most Reliable Path Problem, which finds the path with the lowest probability of failure (Ball et al., 1995b). In Figure 3.1 the path reliability loss function measures a loss of one. Figure 3.1: Path Reliability Loss Function: The loss for this path is 1. Definition 3.3.1 (Path Reliability Loss Function). Lr(x,Y) = max (ij)2A xij(1 Yij) (3.19) Remark 3.3.2. Lr(x,Y) = 0 if no arcs fail and 1 otherwise. Lr(x,Y) Bernoulli(q0) (3.20) where q0 = Õ (ij)2A(x) (1 pij) for all (ij) 2 A. Note that the path reliability loss function is piecewise linear and convex in x, so Formulation 3.2.3 can be further linearized to give the MILP below. Formulation 3.3.3 (Path Reliability Loss Function). Minimize å (ij)2A cijxij (3.21) subject to: 21 å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.22) z + 1 1 b å s2S pszs C (3.23) zs xij(1 ys ij) z 8 s 2 S (3.24) zs 0 8 s 2 S (3.25) x 2 f0, 1gjAj; z 2 R (3.26) 3.3.2 Number of Arc Failures The second loss function counts the number of arcs that fail on the path. For a decisionmaker who is responsible for repairing the path, this loss function gives an estimate of the amount of work to be done. In Figure 3.2, which is the same path as Figure 3.1, the arc failure loss function measures a loss of three. The drawback to this loss function is that not all arc failures, even on the path, may have an equal impact on the path length. For example, in reality it may be very easy to avoid one broken link. Figure 3.2: Arc Failures Loss Function: The loss for this path is 3. Definition 3.3.4 (Arc Failure Loss Function). La(x,Y) = å (ij)2A xij(1 Yij) (3.27) Remark 3.3.5. If pij = p for all (ij) 2 A, then La(x,Y) is the sum of Bernoulli random variables and La(x,Y) Bin(jA(x)j, 1 p) (3.28) 22 Since the arc failure loss function is linear, Formulation 3.2.3 reduces to the MILP below. Formulation 3.3.6 (Arc Failure Loss Function). Minimize å (ij)2A cijxij (3.29) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.30) z + 1 1 b å s2S pszs C (3.31) zs å (ij)2A xij(1 ys ij) z 8 s 2 S (3.32) x 2 f0, 1gjAj; z 2 R jSj + ; z 2 R (3.33) 3.3.3 Number of Detours The final loss function considered in this thesis counts the number of gaps in the current path. A gap is any set of consecutive failed arcs, such as arcs (2, 3) and (3, 4) in Figure 3.3. This loss function is related to ideas from Lagrangian relaxation and derives from the formulation of the problem. Mathematically the loss function measures total violation of flow conservation due to chance in the form of arc failures. This is a useful loss function to examine because it is common to measure the violation of the constraints in stochastic programming, particularly in terms of basic recourse. Figure 3.3 shows that on the same path as Figure 3.1 and Figure 3.2 the detours loss function measures a loss of 2. 23 Figure 3.3: Detours Loss Function: The loss for this path is 2. Definition 3.3.7 (Detours Loss Function). The loss function is Ld(x,Y) = 1 2 å i2N å j2N+(i) xijYij å j2N(i) xjiYji bi (3.34) where bs = 1, bt = 1 and bi = 0 for all i 2 Nnfs, tg. This loss function measures violation of the flow balance constraints or number of detours needed. The loss at some node i 2 N is Ld i (x,Y) = 1 2 å j2N+(i) xijYij å j2N(i) xjiYji bi , (3.35) and Ld(x,Y) = å i2N Ld i (x,Y). (3.36) Since flow conservation will only be violated at the ends of each gap in the path, this loss function ends up counting the number of gaps. The loss function can be interpreted physically as the number of detours necessary to avoid all the gaps. Since one failed arc renders the following arc inaccessible, the failure of that second arc does not cause additional loss. Until all the consecutive failed arcs have all been repaired, that gap in the path must still be detoured around. Lemma 3.3.8. For any incidence vector x of an s t path, Ld(x,Y) = å i2N(x) Ld i (x,Y) Proof. A feasible s t path x corresponds to a graph G(x) = (N(x), A(x)). Nodes i /2 N(x) have no arcs incident at them, and s, t 2 N(x). Hence, 24 Ld i (x,Y) = 1 2 å j2N+(i) xijYij å j2N(i) xjiYji bi 8 i /2 N(x) (3.37) = 1 2 å j2N+(i) 0 Yji å j2N(i) 0 Yij 0 (3.38) = 0 (3.39) So Ld(x,Y) = å i2N(x) Ld i (x,Y) (3.40) Lemma 3.3.9. If x is an incidence vector of an s t path, then Lds (x,Y) = 1 2 jYsi 1j where i 2 N+(s) \ N(x), (3.41) Ldj (x,Y) = 1 2 Yij Yjk 8 j 2 N(x)nfs, tg, (3.42) where i 2 N(j) \ N(x) and k 2 N+(j) \ N(x), and Ldt (x,Y) = 1 2 1 Yjt where j 2 N(t) \ N(x). (3.43) Proof. A feasible s t path x corresponds to a graph G(x) = (N(x), A(x)) such that one arc incident must be incident to s and t and two arcs to all other nodes i 2 N(x)nfs, tg. At s, N(s) = Æ and bs = 1. Also, xsi = 0 if i /2 N(x). Lds (x,Y) = 1 2 å j2N+(s) xsjYsj 0 bs (3.44) = 1 2 j1 Ysk + 0 1j for some k 2 N+(s) (3.45) = 1 2 jYsk 1j (3.46) 25 Similarly, at t, there exists some k 2 N(t) \ N(x) such that, Ld i (x,Y) = 1 2 j0 1 Ykt (1)j (3.47) = 1 2 j1 Yktj (3.48) For j 2 N(x)nfs, tg, let N(j) \ N(x) = fig, and N+(j) \ N(x) = fkg. Ldj (x,Y) = 1 2 å u2N(j) xujYuj å v2N+(j) xjvYjv bj (3.49) = 1 2 1 Yij 1 Yjk 0 (3.50) = 1 2 Yij Yjk (3.51) Corollary 3.3.10. For all j 2 N(x)nfs, tg, Ldj (x,Y) = 1 2 if Yij 6= Yjk; that is, if one arc fails and the other does not. If Yij = Yjk, then Ldj (x,Y) = 0. Lds (x,Y) = 12 if Ysi = 0, and similarly for t. Each gap consisting of consecutive nodes v1 . . . vk 2 N(x) results in a loss of kå i=1 Ldv i (x,Y) = 1. (3.52) Corollary 3.3.11. The sample space for Ld(x,Y) is Wd = f0, 1, . . . , j jN(x)j 2 k g Distribution of Detours Loss Function Theorem 3.3.12. Suppose that all arcs have an equal probability of failure; that is, pij = p for all (ij) 2 A. In that case, the pmf of the detours loss function Ld(x,Y) is given by fL(`; x,Y) = 8>>>>>>>>>>>< >>>>>>>>>>>: qjA(x)j when ` = 0, jA(x)j`+1 å i=` jA(x)ji å j=`1 piqjA(x)ji(i1 `1)(j1 `2)(jA(x)j i j + 1) when ` = 1 . . . j jN(x)j 2 k , and 0 otherwise, (3.53) where q = 1 p. 26 Proof. From Corollary 3.3.11, we know that there are between 0 and j jN(x)j 2 k gaps in the path subgraph G(x) = (N(x), A(x)). Let us then consider the number of ways to distribute the failed arcs into gaps and the surviving arcs into “fragments”. Let Ld(x,Y) = ` 2 Wd be the loss. Suppose ` > 0. Denote by i the number of failed arcs. By Corollary 3.3.10, the failed arcs must be distributed into exactly ` gaps with at least one arc in each gap, which implies that i `. The number of ways to distribute the i failed arcs into ` gaps with at least one arc in each gap is (i1 `1). Denote by j the total number of surviving arcs between the first and last failed arcs. Between each of the ` gaps are ` 1 “fragments” of at least one surviving arc each, which implies that ` 1 j jA(x)j i. The number of ways to distribute the j surviving arcs into ` 1 fragments with at least one arc in each fragment is (j1 `2). The length of the entire path segment from the first to the last failed arc is i + j jA(x)j. This segment may start anywhere in the first jA(x)j i j +1 arcs. The total number of ways to arrange i failed arcs with j jA(x)ji surviving arcs between the failed arcs into ` gaps with ` 1 fragments between them is i 1 ` 1 j 1 ` 2 (jA(x)j i j + 1). (3.54) The probability of i arcs failing and a loss of ` is then P(L(x,Y) = ` andå (ij)2A(x) Yij = i) = piqjA(x)ji jA(x)ji å j=`1 i 1 ` 1 j 1 ` 2 (jA(x)jij+1). (3.55) Hence the total probability of a loss of ` is P(L(x,Y) = `) = jA(x)j`+1 å i=` jA(x)ji å j=`1 piqjA(x)ji i 1 ` 1 j 1 ` 2 (jA(x)j i j + 1). (3.56) 27 The case of L(x,Y) = 0 is a special case. The only possible number of arc failures is i = 0, and all arcs must fall into a single fragment. The probability of ` = 0 is therefore fL(0; x,Y) = qjA(x)j (3.57) As was the case with the path reliability loss function, the detours loss function is piecewise linear and convex. Formulation 3.3.13 (Detours Loss Function). Minimize å (ij)2A cijxij (3.58) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.59) z + 1 1 b å s2S pszs C (3.60) zs å i2N `s i z 8 s 2 S (3.61) `s i å j2N(i) xjiysj i å j2N+(i) xijys ij bi (3.62) `s i bi + å j2N+(i) xijys ij å j2N(i) xjiysj i (3.63) zs 0 8 s 2 S (3.64) x 2 f0, 1gjAj; ` 2 RjSj RjNj; z 2 R (3.65) Table 3.1 summarizes the different loss functions. 28 Table 3.1: Summary of Loss Functions (LF) Name Formulation Sample Space Path Reliability LF Lr(x,Y) = max (ij)2A xij(1 Yij) f0, 1g Arc Failure LF La(x,Y) = å (ij)2A xij(1 Yij) f0, . . . jA(x)jg Detours LF Ld(x,Y) = 12 å i2N j å j2N(i) xjiYji å j2N+(i) xijYij bij f0, . . . j jN(x)j 2 k g 3.4 Comparison of Loss Functions 3.4.1 Magnitude of Losses As shown in Table 3.1, the sample space for the path reliability loss function is a subset of the sample space for the detours loss function, which is in turn a subset of the sample space for the arc failure loss function. Furthermore, the loss as measured by the arc failure loss function will always be greater than or equal to the loss measured by the detours loss function, which will in turn be greater than or equal to the loss measured by the path reliability loss function for the same realization ys 2 W. Theorem 3.4.1. For all feasible paths x and realizations of the random vector of arc failures ys 2 W, Lr(x, ys) Ld(x, ys) La(x, ys). (3.66) Proof. Consider an arbitrary feasible s t path, corresponding to the incidence vector x, and ys, a realization of Y. Let `r = Lr(x, ys), `d = Ld(x, ys), and `a = La(x, ys). If no two arcs fail consecutively, then each arc failure is its own gap in the path, and `d = `a. Otherwise at least two arcs will fail in a row, and `d < `a. Therefore Ld(x, ys) La(x, ys). 29 If `d = 0, no arcs have failed, so `r = 0 and `d = `r. If `d 1, then at least one arc has failed, and `r = 1. When `d 2, therefore, `d > `m. That is, Lr(x, ys) Ld(x, ys). 3.4.2 CVaR as a function of b Before beginning the next chapter, it is useful to examine how the choice of b affects CVaR(x, b). This section gives some guidance that will be used to check the results in the next chapter. Theorem 3.4.2 shows that CVaR(x, b) is piecewise hyperbolic in b but not necessarily convex in b. Theorem 3.4.2. Let L(x,Y) be a loss function with a finite discrete distribution. Given a fixed feasible path x, the corresponding loss distribution is fL(`; x,Y) with sample space WL(x) = f0, 1 . . . `maxg. For each ` 2 WL(x), CVaR(x, b) satisfies: CVaR(x, b) = ` + E n [L(x,Y) `]+ o 1 b (3.67) when P(L(x,Y) < `) < b P(L(x,Y) b). Proof. For each ` 2 WL(x), let b` = P(L(x,Y) `), let b1 = 0, and consider CVaR(x, b) for each b 2 (b`1, b`]. Applying Proposition 2.5.6 gives: CVaR(x, b) = 1 1 b " `(FL(`; x,Y) b) + `max å i=`+1 i fL(i; x,Y) # (3.68) = 1 1 b " `FL(`; x,Y) `b + `max å i=`+1 (i ` + `)fL(i; x,Y) # (3.69) = 1 1 b " `FL(`; x,Y) `b + ` `max å i=`+1 fL(i; x,Y) + `max å i=`+1 (i `)fL(i; x,Y) # (3.70) = 1 1 b " `(1 b) + `max å i=`+1 (i `)fL(i; x,Y) # (3.71) = ` + 1 1 b " `max å i=`+1 (i `)fL(i; x,Y) # (3.72) 30 CHAPTER 4 Computational Experiments The experiments in this chapter are exploratory, meant to illustrate the use and effects of CVaR. The purpose of these computational experiments are to explore the behavior of the three loss functions as C and b vary, to compare the loss functions to each other, and to illustrate the difference between using CVaR and expected value to chose a path though an uncertain network. 4.1 Implementation Details The experiments were performed on two HP xw4600 workstations with a quadcore 2.50 GHz Intel Core 2 CPU and 8 GB of RAM. The operating system on the first computer wasWindows XP and on the otherWindows Vista. The experiments were programmed using Python 2.7 (Van Rossum, 2003) and Gurobi 4.0.2 (Gu et al., 2010). The networks were implemented using the Networkx graph package (Hagberg et al., 2008). The Numpy matrix package (Oliphant, 2006) and the Scipy scientific computing package (Jones et al., 2001) were used for some calculations and for generating scenarios. Data was collected and stored in a Sqlite3 database using the SQLAlchemy objectrelationship manager package (ORM) (Bayer et al., 2011). Finally, results were analyzed and visualized with a mix of Microsoft Access, Microsoft Excel, and Matplotlib (Hunter, 2007), a Python plotting package. 31 4.2 Experiments Test instances were implemented as a directed network with source and sink attributes, a list of integer arc weights, a list of arc failure probabilities, and a set of scenarios. The network edges, arc weights, arc failure probabilities, and scenarios were all stored in the database using Python’s builtin pickle module. Since these experiments are exploratory, smaller instances were chosen. Preliminary experiments had found that instances of several hundred nodes took a long time to solve. Some of the larger and denser networks ran out of memory as well. The final test bed consisted of seven instances divided into three groups. The instances are summarized in Table 4.1. Table 4.1: Summary of Test Bed Name Nodes Arcs Density Scenarios Notes exp0002 6 9 30.0% 512 Complete enumeration branch02 27 30 4.27% 100 Only one arc per branch fails. branch03 27 30 4.27% 100 Less expensive paths are less reliable. branch04 27 30 4.27% 100 pij Uni f (0, 1), cij DU(1, 100) gnp01 25 59 9.83% 100 pij Uni f (0, 1), cij DU(1, 100) gnp02 25 138 23.0% 100 pij Uni f (0, 1), cij DU(1, 100) gnp03 25 243 40.5% 100 pij Uni f (0, 1), cij DU(1, 100) For each instance, the deterministic shortest path (DSP), the most reliable path (MRP, the path with the lowest probability of failure), and the most direct path (MDP, the path with the fewest arcs from source to sink) were calculated and recorded. All the sets of scenarios were generated as random samples of size 100 and treated as equiprobable except exp0002. All scenarios for that instance were 32 enumerated and their exact probability used in optimization. Descriptions of each component of a test instance were also stored in the database (see Appendix B for details). The main experiment mapped the level sets of the optimal cost as a function of b and C for each of the three loss functions. The full CVARSPP formulation for each loss function are in Chapter 3 and Appendix A. The parameter b was varied from 0 to 0.99 in increments of 0.01. The parameter C was varied from 0 to `max in increments of 0.01`max. The optimization status (infeasible, optimal, or unfinished), setup time, solve time, number of nodes, and number of integer feasible solutions found during optimization were recorded for each optimization run. If the run solved to optimality, the objective value, the values of each decision variable, and the arcs, cost, and reliability of the optimal path were recorded as well. For full details of what was recorded, see Appendix B. 4.2.1 Group 1: Complete Enumeration The purpose of the first test network was to enable direct comparison with the analytical results in Section 3.4. The CVARSPP models were still small enough to solve in <1s on average, despite containing all 512 scenarios. Since the network only contained 6 nodes and 9 arcs, it was possible to completely enumerate all st paths. The arc costs and failure probabilities for each path were chosen to force the MDP to be the MRP and have the highest cost and the DSP to be the least reliable. Figure 4.1 contains a diagram of the resulting test network. The lower edges of the level sets of the optimal path length are CVaR(x, b) for each path x that was optimal for the region. The values predicted by Theorem 3.4.2 are plotted over the level sets in Appendix C and agree with the optimal solutions found during the experiment. The level sets are piecewise hyperbolic but not convex in general. The plot for the detours loss function in Figure C.2(b) shows that 33 Figure 4.1: Network exp0002 t 1 2 s 3 4 1, .25 6, .1 2, .2 1, .25 4, .2 2, .25 3, .25 2, .2 2, .2 the lower edge of the MRP contour is flat at `max = 1 for b 0.675. The other two paths represented on the plot have three arcs, so the lower edge of their level sets keep increasing toward `max = 2. 4.2.2 Group 2: Branch Networks The effect of varying the arc failure probabilities on equallength paths was examined in a branch network consisting of a source and sink node connected by five internally nodedisjoint paths of five nodes each, as shown in Figure 4.2. In the first test network, branch02, only one arc was allowed to fail in each branch. As indicated by Theorem 3.4.1, all three loss functions gave the same results for each combination of b and C (Figures C.4(a), (b), and (c)). In the second test network, branch03, the arc costs and arc failure probabilities were set so that more reliable paths also had higher cost (Figure D.2). For instance, each arc on the MRP had cij = 5 and pij = 0.05, each arc on the next most reliable path had cij = 4 and pij = .1, and so on down to the DSP, where each arc had cij = 1 and pij = 0.25. The goal of this setup was to see if each path would be optimal in some region of the (b, C) plane. This was indeed the case under all 34 Figure 4.2: Branch Network Structure 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t three loss functions. In the final branch network, branch04, arc costs and arc failure probabilities were generated randomly for all arcs in order to observe how the CVARSPP models would respond to a wider range of arc reliabilities and path costs. The arc failure probabilities pij were drawn from a Unif(0, 1) distribution and the costs cij from a DU(1, 100) distribution. Since each branch is six arcs long, the expected reliability of each path is low. Even the MRP of branch04 failed in 98 out of 100 scenarios, giving an expected loss of 0.98 under the path failure loss function. When b > 0.98, the minimum choice of C that resulted in a feasible solution to the CVARSPP model under the path reliability loss function is C = 1 as shown in Figure C.8(a). The other two loss functions discriminated among paths, as E Ld(xMRP, ys) = 1.55 and E fLa(xMRP, ys)g = 2.05. This result suggests that the path reliability loss function is better suited for networks with highly reliable arcs or shorter st paths and should potentially be used with smaller b values. 4.2.3 Group 3: Uniform Random Networks The final group of test networks are uniform random networks (Erdös and Rényi, 1959) of 25 nodes with different densities. The arc costs and failure probabilities were randomly generated from pij Unif(0, 1) and cij DU(1, 100). The purpose 35 of this group of test networks was to observe the effect of increasing densities and the behavior of the CVARSPP models on less controlled test networks. As the density increased from 0.0983 for gnp01 to 0.230 for gnp02 to 0.405 for gnp03, the number of paths that were feasible in some region of the (b, C) plane increased as well, from 2 to 4 to 6 (Figures C.10(c), C.12(c), and C.14(c)). Network gnp02 also provided an illustration of the effects of sampling. The path corresponding to the region labeled 53 had two arcs, one of which failed in all 100 scenarios. That arc had a probability of failure of 0.96, so its failure in all scenarios is an artifact of sampling. This resulted in the path being infeasible for C < 1 under the path reliability loss function but being chosen under the detours and arc failure loss functions, as in Figure C.11(d). Under the detours loss function, though, CVaR(x53, b) = 1 for all b 1. If the scenarios had been fully enumerated, the path would be feasible in some narrow region under the path reliability loss function. 4.3 Comparison of Loss Functions To expand on Section 4.2.2 above, another feature of the the path reliability loss function is that, unless the network contains a path with no arcs that can fail, above b = 1 P(YMRP = 0) the instance will be infeasible for all C < 1. The CVARSPP model under the detours loss function has a feasible solution with C < `max for all b 1 in five out of seven networks. Under the arc failure loss function the CVARSPP model has a feasible solution with C < `max for all b 1 for all seven test networks. This suggests that b should be chosen somewhat lower and C somewhat higher within the range (0, `max), when using the path reliability loss function instead of the detours or arc failure loss functions. The detours loss function can use the same values of b as the arc failure loss function for most networks. Only in networks with higher arc failure probabilities does b need to be adjusted 36 slightly downward to preserve feasibility for C < `max. The CVaR limit C, on the other hand, should be chosen lower in the range (0, `max) than under the arc failure loss function, although precise guidance must wait on further computational experiments. Figure C.2 shows how the average optimal path length, average optimal path reliability, average CVaR loss varied as b increased for the test instances using network exp0002. For each loss function, the instances which solved to optimality were grouped by b. Figure C.3(d) records the percentage of instances in the experiment that were feasible for each b and loss function. Figures C.3(a)(c) record the average optimal path reliability, optimal path length, and optimal path CVaR for the instances that solved to optimality. The results all test networks are recorded in Appendix C. 37 CHAPTER 5 Conclusion and Future Work This thesis formulates CVaRbased models for the SSPPPAF. The literature contains models for the SSPPPAF using robust optimization and expected value objectives in stochastic programming (Chapter 2). There has also been some work attempting to trade off robustness and path length (see Ahuja (1988)). Although CVaR has recently been used for supplychain optimization (Madadi et al., 2010), to our knowledge this is the first attempt to use CVaR for a SSPP. The major contributions of this thesis are formulating the CVARSPP and developing three ways of modeling loss. The results found in Chapter 3 develop the loss functions. The three loss functions measure whether the path breaks, the number of gaps, and the number of arc failures with corresponding physical interpretations. Distributions for the each loss function were found for networks with uniform arc failure probabilities. The distribution of the path failure loss function was identified for any path with probabilistic arc failures. Moreover, an analytical characterization of the shape of CVaR(x, b) in terms of b was derived. Finally, it was shown that, for each possible scenario the arc failure loss function is always at least as large as the detours loss function, which is in turn always at least as large as the path reliability loss function. The lower edge of the level sets of the the fullyenumerated test network exp0002 agreed with the values calculated analytically. The computational experiments suggest that the path reliability loss function be used for networks with more reli 38 able paths and the tail threshold b should be set lower than when using the other loss functions. In general, the CVaRbased formulations find more reliable paths than when using expected value. This thesis has focused on modeling. Future research is needed to investigate computational issues. Preliminary experiments have indicated that scaling is likely to prove challenging. A decomposition approach such as the one proposed by KünziBay and Mayer (2006) may be needed to handle larger instances. Another scaling issue is the number of scenarios. Even the small networks in this study showed the effects of sampling. As the network grows, the number of scenarios needed to approximate the actual loss distribution increases rapidly. Some investigation into a better sampling scheme might also help with scaling. Another area for further investigation lies in extending singlestage CVaRbased formulations to twostage recourse models of the SSPPPAF. Finally, more investigation into the behavior of the CVaRbased formulation over a wider and more representative collection of networks is needed. It would be particularly interesting to test the models against some realworld networks, which include their own patterns of arc failures. Relaxing the assumption of independent arc failures and extending the analysis of the loss functions is another avenue to explore. 39 BIBLIOGRAPHY Acerbi, C. and Tasche, D. (2002). On the coherence of expected shortfall. Journal of Banking & Finance, 26(7):1487–1503. Ahuja, R. (1988). Minimum costreliability ratio path problem. Computers & Operations Research, 15(1):83–89. Ahuja, R. K., Magnanti, T. L., and Orlin, J. B. (1993). Network Flows: Theory, Algorithms, and Applications. Prentice Hall, 1 edition. Andreatta, G. and Romeo, L. (1988). Stochastic shortest paths with recourse. Networks, 18(3):193–204. Artzner, P., Delbaen, F., Eber, J.M., and Heath, D. (1999). Coherent Measures of Risk. Mathematical Finance, 9(3):203–228. Ball, M. O., Magnanti, T. L., Monma, C. L., and Nemhauser, G. L., editors (1995a). Network Models, volume 7 of Handbooks in Operations Research and Management Science. Elsevier Science. Ball, M. O., Monma, C. L., and Magnanti, T. L., editors (1995b). Network Routing, volume 8 of Handbooks in Operations Research and Management Science. Elsevier Science Pub Co. Bayer, M., Kirtland, J., de Menten, G., Trier, M., Jenvey, P., Aasma, A., Johnston, P., Ellis, J., and Others (2011). SQLAlchemy 0.7.1. software. Bellman, R. E. (1958). On a Routing Problem. Quarterly of Applied Mathematics, 16:87–90. 40 Bertsimas, D., Brown, D. B., and Caramanis, C. (2008). Theory and applications of robust optimization. Canada, J. R., Sullivan,W. G., Kulonda, D. J., and White, J. A. (2004). Capital Investment Analysis for Engineering and Management. Prentice Hall, 3 edition. Croucher, J. S. (1978). A note on the stochastic shortestroute problem. Naval Research Logistics, 25(4):729–732. Dantzig, G. B. (1957). DiscreteVariable Extremum Problems. Operations Research, 5(2). Deo, N. and Pang, C.Y. (1984). Shortestpath algorithms: Taxonomy and annotation. Networks, 14(2):275–323. Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1(1):269–271. Dreyfus, S. E. (1969). An Appraisal of Some ShortestPath Algorithms. Operations Research, 17(3). Erdös, P. and Rényi, A. (1959). On random graphs, I. Publicationes Mathematicae (Debrecen), 6:290–297. Floyd, R. W. (1962). Algorithm 97: Shortest path. Commun. ACM, 5(6):345+. Gu, Z., Rothberg, E., and Bixby, R. (2010). Gurobi 4.0.2. software. Hagberg, A. A., Schult, D. A., and Swart, P. J. (2008). Exploring network structure, dynamics, and function using NetworkX. In Proceedings of the 7th Python in Science Conference (SciPy2008), pages 11–15, Pasadena, CA USA. Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science and Engineering, 9(3):90–95. 41 Jaillet, P. (1992). Shortest path problems with node failures. Networks, 22(6):589– 605. Jones, E., Oliphant, T., Peterson, P., and Others (2001). SciPy: Open source scientific tools for Python. software. Kendall, K. E. and Kendall, J. E. (2008). Systems analysis and design. Pearson/ Prentice Hall. KünziBay, A. and Mayer, J. (2006). Computational aspects of minimizing conditional valueatrisk. Computational Management Science, 3(1):3–27. Madadi, A. R., Kurz, M. B., Taafe, K., Mason, S., Pohl, E., Root, S., and Sir, M. (2010). Managing Disruptions in Pharmaceutical Supply Chain Networks. In Johnson, A. and Miller, J., editors, Industrial Engineering Research Conference. Oliphant, T. E. (2006). Guide to Numpy. Tregol Publishing, Provo, UT. Papadimitriou, C. (1991). Shortest paths without a map. Theoretical Computer Science, 84(1):127–150. Pollack, M. and Wiebenson, W. (1960). Solutions of the ShortestRoute ProblemA Review. Operations Research, 8(2). Rockafellar, R. and Uryasev, S. (2002). Conditional valueatrisk for general loss distributions. Journal of Banking & Finance, 26(7):1443–1471. Rockafellar, R. T. and Uryasev, S. (2000). Optimization of conditional valueatrisk. Journal of Risk, 2:21–41. Sen, S. and Higle, J. L. (1999). An Introductory Tutorial on Stochastic Linear Programming Models. Interfaces, 29(2). 42 Van Rossum, G. (2003). The Python Language Reference Manual. Network Theory Ltd. Yu, G. and Yang, J. (1998). On the Robust Shortest Path Problem. Computers & Operations Research, 25(6):457–468. 43 APPENDIX A Optimization Models A.1 Parameters These parameters are used in all three optimization models. Graph G = (N, A) is a directed network. S W is the set of scenarios in the model, indexed by s. S may be a random sample from W or the whole sample space if W is small. cij is the cost of arc (ij) 2 A. ys is the realization of the random vector Y corresponding to scenario s 2 S. ps is the probability of scenario s, normalized so that Så s=1 ps = 1. C is the maximum acceptable CVaR loss. b is the percentile for the upper tail (e. g. 0.50, 0.95, 0.99). A.2 Decision Variables xij indicates which arcs are chosen for the shortest path. xij = 1 if arc (ij) is chosen, and 0 otherwise. z 2 R is VaR(x, b) at optimality. zs linearizes the excess loss [L(x, ys) z]+ in scenario s 2 S. `s i is the loss at node i 2 N in scenario s 2 S and linearizes the loss function. 44 A.3 Path Reliability Loss Function Lr(x,Y) = max (ij)2A xij(1 Yij) (A.1) Formulation A.3.1 (Path Reliability Loss Function). Minimize å (ij)2A cijxij (A.2) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (A.3) z + 1 1 b å s2S pszs C (A.4) zs xij(1 ys ij) z 8 s 2 S, (ij) 2 A (A.5) x 2 f0, 1gjAj; z 2 R jSj + ; z 2 R (A.6) 45 A.4 Arc Failure Loss Function La(x,Y) = å (ij)2A xij(1 Yij) (A.7) Formulation A.4.1 (Arc Failure Loss Function). Minimize å (ij)2A cijxij (A.8) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (A.9) z + 1 1 b å s2S pszs C (A.10) zs å (ij)2A xij(1 ys ij) z 8 s 2 S (A.11) x 2 f0, 1gjAj; z 2 R jSj + ; z 2 R (A.12) 46 A.5 Detours Loss Function Ld(x,Y) = å i2N å j2N(i) xjiYji å j2N+(i) xijYij bi (A.13) Formulation A.5.1 (Detours Loss Function). Minimize å (ij)2A cijxij (A.14) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (A.15) z + 1 1 b å s2S pszs C (A.16) zs å i2N `s i z 8 s 2 S (A.17) `s i å j2N(i) xjiysj i å j2N+(i) xijys ij bi 8 s 2 S, i 2 N (A.18) `s i bi + å j2N+(i) xijys ij å j2N(i) xjiysj i 8 s 2 S, i 2 N (A.19) x 2 f0, 1gjAj; z 2 R jSj + ; ` 2 RjSj RjNj; z 2 R (A.20) 47 APPENDIX B Schema of Experiment Database This appendix describes the database used in the computational experiments. The database contains all data generated during the experiments. The database schema is diagrammed in Figure B.2 using crow’s foot notation. Each table in the database and its fields are described below. The database was normalized to the Backus Naur Third Normal Form (Kendall and Kendall, 2008) before being selectively denormalized for ease of use. B.1 Database Tables Graph The table for the directed networks. This does not contain arc weights, arc failure probabilities, or scenarios. name The name of the directed network. descr A sentence or two describing the network to a human reader. n The number of nodes. a The number of arcs. source The source node of the network. sink The sink node of the network. density The density of the network obj The Python object containing the arc list. Weights The table for arc weights (cij). 48 name A short name for the arc weights. descr Asentence or two describing the arc weights or their generation method to a human reader. obj The Python object containing the list of arc weights. length The number of weights in the list of arc weights. graph_id Foreign key to Graphs table. AFPs The table for arc failure probabilities (pij). name A short name for the arc failure probabilities. descr A sentence or two describing the arc failure probabilities and their generation method to a human reader. obj The Python object containing the arc failure probabilities. length The number of arc failure probabilities in the list. graph_id Foreign key to Graphs table. Scenarios The table for scenarios (ys and ps for all s 2 S) name A short name for this set of scenarios. descr Asentence or two describing the scenarios and their generation method to a human reader. s The number of scenarios. scen_array The Python object containing an jAj s Numpy matrix of the scenarios. Each row is ys, the incidence vector of arc failures. ps The Python object containing a Numpy vector of scenario probabilities. Each element is corresponds to ps. empafp The Python object containing a Numpy vector of jAj entries. Each entry is the percent of scenarios in which that arc fails. 49 afp_id Foreign key to AFPs table. TestInstance The table for instances. This is a full network scenarios and with weights and failure probabilities attached to each arc. name A short name for this instance. graph_id Foreign key to Graphs table. wts_id Foreign key to Weights table. scen_id Foreign key to Scenarios table. Path The table for information about paths of interest through a network. rel The apriori reliability (i.e. Õ (ij)2A(x) (1 pij)). cost The cost of the path (i.e. å (ij)2A(x) cij). arcs Number of arcs in this path. type Classification—DSP, MRP,MDP, or an optimal solution to some CVARSPP instance. inst_id Foreign key to the TestInstance table. sol_id Optional: Foreign key to the Solution table. PathArcs The table that stores the arcs in a path. head The head of an arc in a path. tail The tail of an arc in a path. path_id Foreign key the the Path table indicating which path contains these arcs. Exp The table of summary information about a computational experiment. The information in this table is mostly intended as notes for a human reader. 50 name A short name for this set of experiment. descr A summary of the experiment. goal The objective of the experiment. beta_method A summary of how the tail threshold probabilities b are to be generated. c_method A summary of how the maximum acceptable losses C are to be generated. loss_fns The loss functions involved in this experiment. computer The machine used to run the experiment. solver The program used to solve the CVARSPP models. Param The table storing solver control parameters such as time limit. filename The name of the file with the parameters for the solver to read. folder The path to the file. suffix The extension of the file indicates format. plist The Python object containing the parameters, a duplicate of the file. Run Each record represents a unique instance of the CVARSPP. name A short name for this CVARSPP instance. beta The tail threshold probability b. C The Python object representing the arc failure probabilities (a BLOB in SQL terms). loss_fn Which loss function to use; i.e., which specific model to choose (see Appendix A). solved Flag indicating whether the run has been sent to the solver or not. 51 exp_id Foreign key to Exp table. param_id Foreign key to the Param table. inst_id Foreign key to TestInstance table. OptInfo The table for information reported by the solver about an instance. If an optimal solution is found, it is stored in the Solution table. optimal Did the instance solve to optimality? grb_status The status code reported by the solver. setup The time to set up each instance. sol_time The solution time for each instance. iter The number of simplex iterations performed for all the nodes explored. nodes The number of branchandbound nodes explored. best_bound If the solver ran out of time, the best bound found during optimization. num_sol The number of integer feasible solutions during optimization. obj The Python object containing all the information about running the optimization which was returned by the solver. date Timestamp for the optimization run. run_id Foreign key for the Run table. Solution The table for information about the optimal solution to a run. obj_cost The cost of the optimal path. var The VaR for this instance. cvar The CVaR for the instance, calculated C and the slack in the CVaR constraint. 52 cvar_slack The slack in the CVaR constraint. dvs The Python object containing the values of all decision variables. run_id Foreign key for the Run table. B.2 Database Schema Figure B.1: Legend: Crow’s Foot Notation Primary Key Attribute (FK): Foreign Key (O): Optional Table Name Relationship Symbol Zero or onetozero or one Zero or onetomany Onetozero or many Onetomany 53 Figure B.2: Database Schema for Computational Experiments after Kendall and Kendall (2008) Graph id name descr n a source sink density obj Weights id name descr obj length Scenarios id name descr s scen_array ps empafp afp_id (FK) AFPs id name descr obj length TestInstance id name graph_id (FK) wts_id (FK) scen_id (FK) belongs to / is constructed from is used to generate / is generated from is used in / uses Exp id name descr goal beta_method c_method loss_fns computer solver Run id name beta C loss_fn solved exp_id (FK) param_id (FK) inst_id (FK) obj_cost var cvar cvar_slack dvs run_id (FK) OptInfo id optimal grb_status setup sol_time iter nodes best_bound num_sol obj date run_id (FK) Param id filename (O) folder (O) suffix (O) plist (O) contains / is optimal for Path id rel cost arcs type inst_id (FK) sol_id (O) (FK) PathArc id head tail path_id (FK) is comprised of / are on is solved during / solves the SSPPAF on has / belongs to of is used by / passes these to Gurobi records optmization information in / is reported by Gurobi after finds / is extracted from 54 APPENDIX C Objective Function Contour Graphs by Loss Function This appendix contains graphs of the objective function values for each test instance and loss function. 55 Figure C.1: exp0002: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C CVaR(MRP, b) CVaR(x6, b) CVaR(DSP, b) (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C 56 Figure C.2: exp0002: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 0 2 4 6 8 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 1 1.5 2 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 57 Figure C.3: branch02: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C 58 Figure C.4: branch02: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 5 102 0.1 0.15 0.2 0.25 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 240 260 280 300 320 b Average Length (c) Average CVaR for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.4 0.6 0.8 b Average CVaR (d) Infeasible instances as a function of b 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 59 Figure C.5: branch03: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 2 4 6 b Maximum Loss C 60 Figure C.6: branch03: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 10 15 20 25 30 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 1 1.5 2 2.5 3 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 61 Figure C.7: branch04: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 2 4 6 b Maximum Loss C 62 Figure C.8: branch04: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 200 250 300 350 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 1 2 3 4 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 63 Figure C.9: gnp01: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 2 4 6 b Maximum Loss C 64 Figure C.10: gnp01: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 52.5 53 53.5 54 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 0.5 1 1.5 2 2.5 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 65 Figure C.11: gnp02: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 b Maximum Loss C 66 Figure C.12: gnp02: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 50 55 60 65 70 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 0.5 1 1.5 2 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 67 Figure C.13: gnp03: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.5 1 1.5 2 b Maximum Loss C (c) Arc failure loss function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 1 2 3 b Maximum Loss C 68 Figure C.14: gnp03: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 60 80 100 120 140 160 b Average Length (c) Average CVaR for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.4 0.6 0.8 1 1.2 1.4 b Average CVaR (d) Infeasible instances as a function of b 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.4 0.6 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 69 APPENDIX D Test Network Diagrams D.1 Small Network: exp0002 The scenarios for this small network are completely enumerated. Arc failure probabilities are chosen to force the MDP to be the MRP and have the highest cost while the DSP is the least reliable path. t 1 2 s 3 4 1, .25 6, .1 2, .2 1, .25 4, .2 2, .25 3, .25 2, .2 2, .2 Figure D.1: Network exp0002 70 D.2 Branch Networks The three following networks have five disjoint paths, or branches, of six arcs connecting the source s to sink t. Figure D.2: Network branch02 The arc weights cij are drawn from DU(1, 100). One arc is allowed to fail on each branch with probability pij drawn from Unif(0, 1). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t 28 86 20 51 15 98 18 74 58 39 31 84 94 31 54 64 70 55 32 87 22 49 41 55 43, .12 8 11, .22 33, .68 8, .12 85, .20 Figure D.3: Network branch03 The arc failure probabilities and arc costs are uniform across each branch. As the cost of the path increases from the DSP to the most expensive path, so does the reliability. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t 1, .25 2, .2 3, .15 4, .1 5, .05 1, .25 2, .2 3, .15 4, .1 5, .05 1, .25 1, .25 1, .25 1, .25 2, .2 2, .2 2, .2 2, .2 3, .15 3, .15 3, .15 3, .15 4, .1 4, .1 4, .1 4, .1 5, .05 5, .05 5, .05 5, .05 71 Figure D.4: Network branch04 The arc failure probabilities pij are drawn from Unif(0, 1), and the arc weights cij are drawn from DU(1, 100). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t 31, .25 2, .47 35, .47 53, .25 91, .37 48, .8 53, .13 13, .99 84, .64 34, .67 17, .16 18, .4 67, .62 26, .96 53, .087 87, .94 74, .38 32, .072 96, .4 72, .65 67, .39 56, .34 54, .026 93, .66 15, .75 35, .29 60, .91 79, .44 88, .71 14, .032 D.3 Uniform Random Networks The final three networks are GNP(n, p) random networks (Erdös and Rényi, 1959). All three networks have n = 25 nodes. The arc failure probabilities pij are drawn from Unif(0, 1), and arc weights cij from DU(1, 100). Figure D.5: Network gnp01 — GNP(25, 0.1) t 12 11 10 9 6 7 8 5 4 3 2 0 s 13 14 15 16 17 18 19 20 21 22 23 72 Figure D.6: Network gnp02 — GNP(25, 0.25) t 12 11 10 9 6 7 8 5 4 3 2 1 s 13 14 15 16 17 18 19 20 21 22 23 Figure D.7: Network gnp03 — GNP(25, 0.4) t 12 11 10 9 6 7 8 5 4 3 2 1 s 13 14 15 16 17 18 19 20 21 22 23 73 VITA Juliana Bright Candidate for the Degree of Master of Science Thesis: ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK Major Field: Industrial Engineering and Management Biographical: Personal Data: Born in London, UK on April 22, 1980. Education: Bachelor of Arts in German from the University of Oklahoma, December 2005. Bachelor of Science in Industrial Engineering and Management from Oklahoma State University in Stillwater, Oklahoma in August 2009. Completed the Requirements for the Master of Science degree at Oklahoma State University in December 2011 with a major in Industrial Engineering and Management. Experience: Graduate research assistant, Industrial Engineering and Management, Oklahoma State University, Fall 2010Spring 2011; undergraduate and graduate teaching assistant, CEATLabs, Oklahoma State University, January 2008August 2010; staff assistant, DollarThrifty Automotive Group, Tulsa, OK, June 2004September 2006. Professional Memberships: Engineering Intern, Institute for Operations Research and Management Science (INFORMS), Institute of Industrial Engineers (IIE) Awards and Honors: Doctoral Academy Fellowship, University of Arkansas; First place in the 2010 IIE Undergraduate Student Technical Paper Competition. Name: Julie Anna Bright Date of Degree: December 2011 Institution: Oklahoma State University Location: Stillwater, Oklahoma Title of Study: ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK Pages in Study: 73 Candidate for the Degree of Master of Science Major Field: Industrial Engineering and Management Finding a shortest path in a network is a classical problem in discrete optimization. The systems underlying the network models are subjects to a variety of sources of uncertainty. The purpose of this thesis is to model the shortest path problem with probabilistic arc failures with the aim of finding short but reliable paths through the network. This thesis proposes to use Conditional ValueatRisk (CVaR) to find a path that is robust under probabilistic arc failures. CVaR, a quantitative risk measure, is roughly the mean excess loss associated with a decision. This thesis also develops three models of losses due to arc failures. Mixed integer linear programming models for the stochastic shortest path problem with probabilistic arc failures are formulated and implemented. The optimization models limit the CVaR of the loss due to arc failures asmeasured by the models of loss developed in the thesis. ADVISOR’S APPROVAL: Dr. Balabhaskar Balasundaram
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Robust Shortest Paths under Uncertainty Using Conditional ValueatRisk 
Date  20111201 
Author  Bright, Julie Anna 
Keywords  cvar, modeling, networks, optimization, shortest path, stochastic optimization 
Department  Industrial Engineering & Management 
Document Type  
Full Text Type  Open Access 
Abstract  Finding a shortest path in a network is a classical problem in discrete optimization. The systems underlying the network models are subjects to a variety of sources of uncertainty. The purpose of this thesis is to model the shortest path problem with probabilistic arc failures with the aim of finding short but reliable paths through the network. This thesis proposes to use Conditional ValueatRisk (CVaR) to find a path that is robust under probabilistic arc failures. CVaR, a quantitative risk measure, is roughly the mean excess loss associated with a decision. This thesis also develops three models of losses due to arc failures. Mixed integer linear programming models for the stochastic shortest path problem with probabilistic arc failures are formulated and implemented. The optimization models limit the CVaR of the loss due to arc failures as measured by the models of loss developed in the thesis.</ 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK By JULIE ANNA BRIGHT Bachelor of Science in Industrial Engineering and Management Oklahoma State University Stillwater, OK 74075 2009 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE December 2011 COPYRIGHT c By JULIE ANNA BRIGHT December 2011 ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK Thesis Approved: Balabhaskar Balasundaram Thesis Advisor Manjunath Kamath Ricki Ingalls Sheryl Tucker Dean of the Graduate College iii TABLE OF CONTENTS Chapter Page 1 Introduction 1 1.1 Organization of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 The Deterministic Shortest Path Problem . . . . . . . . . . . . . . . . 2 1.2.1 Formal Definition . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Classical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.3 Linear Programming Formulation . . . . . . . . . . . . . . . . 4 1.3 Stochastic Shortest Path Problems . . . . . . . . . . . . . . . . . . . . . 5 1.4 Research Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 Research Tasks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Review of Modeling Approaches 8 2.1 Stochastic Programming Models . . . . . . . . . . . . . . . . . . . . . 8 2.2 Robust Optimization Models . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Other Modeling Approaches . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Models Based on Conditional ValueatRisk . . . . . . . . . . . . . . . 11 2.5 Formal Definition of CVaR . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.6 Using CVaR in Optimization . . . . . . . . . . . . . . . . . . . . . . . . 15 3 Model Development and Analysis 17 3.1 Modeling Arc Failures . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Formulation of the CVARSPP . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.1 Additional Parameters . . . . . . . . . . . . . . . . . . . . . . . 18 iv 3.2.2 CVARSPP Optimization Model . . . . . . . . . . . . . . . . . . 18 3.3 Loss Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.1 Path Reliability . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.2 Number of Arc Failures . . . . . . . . . . . . . . . . . . . . . . 22 3.3.3 Number of Detours . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 Comparison of Loss Functions . . . . . . . . . . . . . . . . . . . . . . 29 3.4.1 Magnitude of Losses . . . . . . . . . . . . . . . . . . . . . . . . 29 3.4.2 Dependence of CVaR on b . . . . . . . . . . . . . . . . . . . . . 30 4 Computational Experiments 31 4.1 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2.1 Group 1: Complete Enumeration . . . . . . . . . . . . . . . . . 33 4.2.2 Group 2: Branch Networks . . . . . . . . . . . . . . . . . . . . 34 4.2.3 Group 3: Uniform Random Networks . . . . . . . . . . . . . . 35 4.3 Comparison of Loss Functions . . . . . . . . . . . . . . . . . . . . . . 36 5 Conclusion and Future Work 38 Bibliography 39 A Optimization Models 44 A.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 A.2 Decision Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 A.3 Path Reliability Loss Function . . . . . . . . . . . . . . . . . . . . . . . 45 A.4 Arc Failure Loss Function . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.5 Detours Loss Function . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 B Database Schema 48 B.1 Database Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 v B.2 Database Schema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 C Objective Function Graphs 55 D Test Network Diagrams 70 D.1 Small Network: exp0002 . . . . . . . . . . . . . . . . . . . . . . . . . . 70 D.2 Branch Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 D.3 Uniform Random Networks . . . . . . . . . . . . . . . . . . . . . . . . 72 vi LIST OF TABLES Table Page 3.1 Summary of Loss Functions (LF) . . . . . . . . . . . . . . . . . . . . . 29 4.1 Summary of Test Bed . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 vii LIST OF FIGURES Figure Page 2.1 An example of VaR and CVaR . . . . . . . . . . . . . . . . . . . . . . . 12 3.1 Path Reliability Loss Function: The loss for this path is 1. . . . . . . . 21 3.2 Arc Failures Loss Function: The loss for this path is 3. . . . . . . . . . 22 3.3 Detours Loss Function: The loss for this path is 2. . . . . . . . . . . . 24 4.1 Network exp0002 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Branch Network Structure . . . . . . . . . . . . . . . . . . . . . . . . . 35 B.1 Legend: Crow’s Foot Notation . . . . . . . . . . . . . . . . . . . . . . 53 B.2 Database Schema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 C.1 exp0002: Level sets of the optimal path length . . . . . . . . . . . . . 56 C.2 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 57 C.3 branch02: Level sets of the optimal path length . . . . . . . . . . . . . 58 C.4 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 59 C.5 branch03: Level sets of the optimal path length . . . . . . . . . . . . . 60 C.6 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 61 C.7 branch04: Level sets of the optimal path length . . . . . . . . . . . . . 62 C.8 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 63 C.9 gnp01: Level sets of the optimal path length . . . . . . . . . . . . . . . 64 C.10 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 65 C.11 gnp02: Level sets of the optimal path length . . . . . . . . . . . . . . . 66 C.12 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 67 viii C.13 gnp03: Level sets of the optimal path length . . . . . . . . . . . . . . . 68 C.14 Change in optimal path properties . . . . . . . . . . . . . . . . . . . . 69 D.1 Network exp0002 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 D.2 Network branch02 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 D.3 Network branch03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 D.4 Network branch04 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 D.5 Network gnp01 — GNP(25, 0.1) . . . . . . . . . . . . . . . . . . . . . . 72 D.6 Network gnp02 — GNP(25, 0.25) . . . . . . . . . . . . . . . . . . . . . 73 D.7 Network gnp03 — GNP(25, 0.4) . . . . . . . . . . . . . . . . . . . . . . 73 ix CHAPTER 1 Introduction to Shortest Path Problems Under Uncertainty The shortest path problem (SPP) is to find the path through a network with the lowest “cost” — whether that cost is measured in miles, money, minutes, or any other metric. This problem has been extensively studied over the past 50 years (Pollack and Wiebenson, 1960; Dreyfus, 1969; Deo and Pang, 1984). When the structure and costs of the network are known, then very fast and powerful tools exist to find a shortest path, such as Dijkstra’s algorithm (Dijkstra, 1959). In the real world, though, networks are afflicted with failures, and costs can fluctuate randomly. This limits the applicability of classical, deterministic SPP algorithms. Often the definition of a “shortest” path needs to be adjusted. Different mathematical models and algorithms are needed to handle uncertainty. The family of shortest path problems under uncertainty (modeled by probability distributions) are collectively called stochastic shortest path problems (SSPPs). This thesis investigates one type of SSPP. In this SSPP, the arcs of a network can fail randomly and independently with some probability. This model is useful for roads that close due to accidents or severe weather or for PERT charts where activities may not get done. The effect of arc failures on the overall network and the shortest path will be investigated in order to create measures to quantify the loss of an arc. Secondly, a quantitative measure of downside risk called Conditional ValueatRisk (CVaR) will be used to find short and robust paths in the network (Rockafellar and Uryasev, 2002). 1 1.1 Organization of this Thesis Chapter 1 discusses the classical or deterministic shortest path problem (SPP), including a formal definition and a review of the most common algorithms, before providing some motivating examples to show how uncertainty in graphs affects the SPP. This chapter also introduces the research problem. Chapter 2 reviews the modeling approaches in the literature to show what kinds of uncertainty have been captured and what kinds of uncertainty remain relatively lessstudied. This chapter also introduces the concept of CVaR. Chapter 3 formulates the CVaRConstrained Stochastic Shortest Path Problem (CVARSPP) and develops three loss functions to model loss due to probabilistic arc failures to use in the CVARSPP. This chapter also analyzes the loss functions, deriving their distributions for uniform arc failures and some implications for their behavior in CVARSPP models. Then, Chapter 4 presents some exploratory computational experiments that illustrate the behavior of the CVARSPP models with the loss functions, leading to some preliminary guidance on setting parameters. Finally, Chapter 5 presents the conclusions of this thesis and some directions for further research. 1.2 The Deterministic Shortest Path Problem 1.2.1 Formal Definition Consider the network G = (N, A), consisting of N, a set of nodes, and A N N, a set of directed arcs connecting pairs of nodes. For each arc (ij) 2 A, node i is called the tail and node j is called the head. Each arc also has an associated cost cij 2 R+. For each node i 2 N, its inneighborhood is defined as N(i) = fjj(ji) 2 Ag and its outneighborhood is defined as N+(i) = fjj(ij) 2 Ag. Adirected path P is a sequence of k distinct nodes (i1, i2, . . . , ik), where (ij, ij+1) 2 2 A for all j 2 f1, . . . , k1g. Throughout the remainder of this thesis, directed paths will be referred to simply as “paths”. The edgeinduced subgraph of the path P is denoted by G(P) = (N(P), A(P)). The path P can also be represented by an incidence vector x 2 f0, 1gjAj of arcs in P, which induces the subgraph G(x) = (N(x), A(x)). The cost, also called weight or length, of a path P is cP = å (ij)2A(P) cij. Definition 1.2.1. (Deterministic Shortest Path Problem) Given network G = (N, A) and nodes s, t 2 N, the PointtoPoint Shortest Path Problem is to find a minimum cost path from node s to t in G. The SingleSource Shortest Path Problem is to find a shortest path from a source node s 2 N to each node i 2 Nnfsg. The AllPairs Shortest Path Problem is to find a shortest path from every node i 2 N to every other node j 2 N. 1.2.2 Classical Results for Shortest Path Problem In the absence of negative cost cycles in the network, the SPP is a wellsolved problem, for which several stronglypolynomial algorithms are known (Deo and Pang, 1984). These algorithms assume that the weight of each arc in the network is constant, and network structure (nodes and arcs) are known with certainty. Algorithms for the SPP fall into two main categories: labelsetting and labelcorrecting (Deo and Pang, 1984). Dijkstra’s algorithm (Dijkstra, 1959) is an example of a labelsetting algorithm. The BellmanFord algorithm (Bellman, 1958) is a labelcorrecting algorithm derived from dynamic programming. The labelsetting algorithms require all arc costs to be nonnegative, while the labelcorrecting algorithms only require that the network have no negative cost cycles. The Floyd Warshall algorithm for the AllPairs Shortest Path Problem (Floyd, 1962) is also based on dynamic programming. The SingleSource and PointtoPoint Shortest Path Problems can also be formulated as a linear program and solved with the net 3 work simplex method when the network contains no negative cost cycles (Dantzig, 1957). 1.2.3 Linear Programming Formulation of the Shortest Path Problem The linear programming (LP) formulation of PointtoPoint SPP optimization model is worth presenting on its own since this thesis extends this formulation. The constraint matrix for the LP formulation of the SPP is totally unimodular, so the extreme point optimal solutions are integral (Ahuja et al., 1993). Inputs Network G = (N, A) is a directed network. s 2 N is the source node. t 2 N is the sink node. cij is the cost of arc (ij) 2 A. Decision Variables xij for all (ij) 2 A. At any x 2 f0, 1gjAj, xij = 1 implies that arc (ij) is in the solution and xij = 0 implies that arc (ij) is not. Formulation 1.2.2 (Deterministic Formulation). Minimize å (ij)2A cijxij (1.1) 4 subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (1.2) xij 0 8 (ij) 2 A (1.3) 1.3 Stochastic Shortest Path Problems Uncertainty in the system underlying a network model may be modeled by random variation in network parameters such as arc costs and in the structure of the network, particularly node and arc failures. Once uncertainty is allowed into the network model, there is no longer one welldefined shortest path problem. A family of Stochastic Shortest Path Problems (SSPPs) arise instead, depending on which elements in the network are uncertain and how the uncertainty is modeled. Furthermore, SSPPs generally require the redefinition of the concept of feasible and optimal solutions. One SSPP models the arc cost cij as a random variable. Applications of this model include modeling the travel time for an arc in a congested road network or the task duration in a PERT network (Ball et al., 1995a). For this problem, two ways to redefine the shortest path are in terms of expected length, which gives rise to many stochastic programming approaches, and worstcase length, in which arc costs are set at or close to their maximum possible values, which lies behind many robust optimization approaches (Ball et al., 1995b). Uncertainty about which node is the source or sink can be modeled by randomly assigning a supply of +1 to one node or 1 to another. This would be useful in modeling ambulance routing during an emergency when the driver dis 5 covers which ER has room for the patient enroute. Allowing uncertainty in the structure of the network (e.g. uncertainty in the set of nodes or arcs, or only partial knowledge) gives rise to another set of SSPPs. One source of uncertainty is limited knowledge of the network, such as not knowing the cost or head node of an arc until the tail is reached. This is frequently used to model hideandseek or adversarial search on networks (Papadimitriou, 1991) or for navigating real road networks without a map, such as after a major disaster (Papadimitriou, 1991). Secondly, the nodes may be probabilistic—that is, they can fail. This is used in modeling telecommunications networks or computer networks (Jaillet, 1992). Finally, the arcs may fail with some probability, giving rise to the subject of this thesis. Allowing arcs to fail allows modeling a wide variety of disturbances in networks. For instance, road closure in transportation networks, due to accidents, weather, or other events, can be modeled as an arc failure. Cutting through the cables in fiber optic networks, which happens during construction or fires (Ball et al., 1995a), can also be modeled as arc failure. In a social network, an onagain offagain friendship between two people can be modeled as an arc that fails after an argument and is restored after a reconciliation. It is therefore desirable to find a short, reliable path through a network. The redefinition of the shortest path depends on the trade off between reliability and shortness. 1.4 Research Problem This thesis investigates models of the Stochastic Shortest Path Problem with Probabilistic Arc Failures (SSPPPAF) using Conditional ValueatRisk (CVaR) constraints. CVaR is a risk management tool that can trade off risk aversion and performance, so it is able to handle varying levels of risk tolerance. Moreover, CVaR is based on an explicit loss function, which allows more precision in modeling the effects 6 of arc failures. The concept of CVaR is discussed in more detail in Chapter 2. This thesis will limit itself to arcs that fail randomly and independently of each other with some probability. All other graph parameters are deterministic. 1.5 Research Tasks The objective of this thesis is to develop CVaRbased optimization models for the SSPPPAF. The solutions to these models should be paths that are short and robust under arc failures; specifically, the loss due to arc failure as measured by using CVaR should be low. To achieve this objective, the tasks listed below will be pursued. 1. Model the effects of arc failures using different loss functions. 2. Develop and implement optimization models to find the shortest path subject to CVaR constraints. 3. Investigate the effect of different loss functions and CVaR bound on the feasibility of the optimization model. 4. Gain insight into the difference between models using expected value and CVaR via their behavior with different loss functions on carefullychosen test instances. 7 CHAPTER 2 A Review of Modeling Approaches for the Shortest Path Problem under Uncertainty 2.1 Stochastic Programming Models Stochastic programming is one of the approaches available to model optimization problems under uncertainty. Stochastic programming models typically use expected value to model the effect of uncertainty in the objective function and to allow for recourse actions (Sen and Higle, 1999). As a measure of risk, expected value is risk neutral. It weights gains as high as losses whereas human decisionmakers are often much more riskaverse (Canada et al., 2004), preferring to minimize losses rather than maximize gains. There have been a few stochastic programming models of the SSPPPAF. Croucher (1978) incorporated recourse decisions into a labelcorrecting algorithm for directed acyclic networks using dynamic programming. The paper attempts to model a traveler who has imperfect knowledge of a road network. In particular, the traveler does not know whether an arc is blocked until arriving at its head. The objective of the algorithm is to find a path with the shortest expected length from each node to the sink node given that the desired arc is unreliable. In the algorithm, each node is labeled with the expected shortest path length to the sink. Each arc in the network has a probability pij of being usable. At each node an outgoing arc is selected. If the desired arc is broken, one of the other outgoing arcs of the node is selected at random. When the node has only one outgoing arc, it is not allowed to fail. Each node is labeled with the expected shortest path length to the sink, which 8 is the average of the labels of the outneighbors weighted by their probability of transversal. This algorithm is quick to compute, but it relies on the assumption that at least one outgoing arc at each node will be usable. Andreatta and Romeo (1988) built on Croucher (1978)’s paper to explicitly find optimal recourse for a blocked arc, and they find a path that minimizes the expected length of the path including recourse, which they define as detours around a failed arc. The model uses stochastic dynamic programming to solve the problem. This model suffers from the curse of dimensionality as the number of scenarios and the number of failed arcs per scenario grow. Like the previous paper, this one seeks only to minimize the expected path length and does not give special consideration to the extremes of the distribution. 2.2 Robust Optimization Models Robust optimization is another approach to optimization under uncertainty. Robust optimization models seek to find solutions that are feasible no matter how the underlying uncertainty resolves (Bertsimas et al., 2008). The uncertainty is captured in an uncertainty set, which is a collection of possible future states of the system being modeled. The optimal solution must be feasible for every element in the uncertainty set. The uncertainty set need not contain every possible resolution of the uncertainty, only the subset for which the the optimal solution must remain feasible. Robust optimization tends to favor minimax objectives (Bertsimas et al., 2008). Recently Yu and Yang (1998) showed a version of the robust shortest path problem to be NPcomplete. The paper considers a network with uncertain arc costs, which are captured in a set of scenarios, and develops two robust optimization models for the problem. The first model uses a minimax objective over the scenarios. The second model seeks a path with the smallest possible range of lengths. 9 Although the authors develop a dynamic programming algorithm for networks with bounded numbers of scenarios, the algorithm is pseudopolynomial in the number of scenarios. As the number of scenarios needed to capture the behavior of a large graph grows very quickly, this algorithm is not suitable for large scale uses (Yu and Yang, 1998). 2.3 Other Modeling Approaches A conceptually similar approach to the SSSPPAF is the Minimum Cost Reliability Ratio Path Problem (MCRRPP), a bicriteria optimization model proposed by Ahuja (1988). This problem considers paths from a single source node to a single sink in directed network with arc failures. The objective is to find a path with the optimal ratio of cost to reliability, defined as the probability of no arc failures on the path. The author first develops optimality criteria for a bicriteria optimization problem with the criteria of cost and reliability. The paper develops a pseudopolynomial dynamic programming algorithm with a complexity of O(mnDlogm) to solve the problem, although the computational experiments found that averagecase complexity was much lower. The MCRRPP trades off cost and reliability in a riskneutral manner. By using reliability in the formulation of the MCRRPP, the model implicitly a loss of 100% when any arc on the path fails, and 0% otherwise. This approach is similar to the approach used in this thesis in that it uses arc failures and explicitly trades off cost and reliability. In a related approach, Jaillet (1992) incorporates node failures into a SSPP model. The goal of this model is to find an easily repairable path of minimum expected length given node failures. The paper shows that the SSSP with probabilistic node failures is in general NPcomplete although some restricted types of the problem are solvable in polynomial time. 10 2.4 Models Based on Conditional ValueatRisk CVaR was developed in the late 1990s as an improvement on an earlier financial risk measure called ValueatRisk or VaR (Acerbi and Tasche, 2002). VaR has been very popular, even being written into regulations in several countries (Rockafellar and Uryasev, 2000), but it has a number of drawbacks (Rockafellar and Uryasev, 2000). For most probability distributions, VaR is difficult to compute and is not a coherent risk measure (Artzner et al., 1999). CVaR, on the other hand, is a coherent measure of risk and is easier to calculate (Acerbi and Tasche, 2002), where it is called expected shortfall. More pertinently, Rockafellar and Uryasev presented a way to use CVaR in optimization problems via a portfolio optimization example (Rockafellar and Uryasev, 2002). Intuitively, CVaR is the average of the upper tail of the distribution of losses due to uncertainty, which are quantified in a loss function. The loss function and the threshold for the upper tail are provided by the modeler. Figure 2.1 illustrates VaR and CVaR for a negative binomial distribution with the tail threshold set at the 90th percentile. Currently CVaR is applied in several areas of finance such as insurance and credit risk evaluation (Rockafellar and Uryasev, 2000). In these applications, risk is usually specified in terms of monetary losses (Rockafellar and Uryasev, 2000) although utility can be used as well. In more general contexts, though, CVaR is defined in terms of a realvalued loss function L(x,Y), where Y is a vector of random variables that model uncertainties in the underlying system, and x is a vector of decisions. In portfolio optimization, Y would be future share prices for all companies traded in the stock market and x the portfolio allocation. For each decision x, the loss function L(x,Y) is itself a random variable. Generally gains are modeled as negative losses. In portfolio optimization, the loss is often measured in terms of 11 Figure 2.1: An example of VaR and CVaR 0 5 10 15 20 25 30 35 40 45 50 55 60 0 1 102 2 102 3 102 4 102 5 102 6 102 CVaR pmf 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 VaR cd f return on investment (ROI), where the loss is negative if the ROI is positive. A recent paper (Madadi et al., 2010) formulates a CVaRbased supply chain network model in order to minimize the risk of defective pharmaceuticals in a medical supply chain. Supplier failures and resulting defective pharmaceuticals are the sources of uncertainty in this network. The demand at the customer nodes is fixed. In the transportation network the arcs and their costs are deterministic. The objective function minimizes the CVaR of the loss from defective pharmaceuticals. The constraints trade off the cost of inspection versus the cost of defective products. They enumerate all possible scenarios and assign them a probability. This was computationally feasible because the number of suppliers is likely to be small compared to the number of customers. One strength of their model is that it allows for a state of partial failure. This paper uses CVaR for managing productionside risks. In contrast to this thesis, modeling uncertainty in the arcs of the underlying 12 transportation network is not within the scope of Madadi et al. (2010)’s paper. 2.5 Formal Definition of CVaR Definition 2.5.1 (Rockafellar and Uryasev (2002)). Let x 2 X Rn be a vector of decision variables and Y : W ! Rm be a vector of random variables defined on the sample space W. Then the loss function L(x,Y) is a function of the decisions x and the random variables Y. The random variable L(x,Y) has a family of distributions, one for each x 2 X. The sample space of L(x,Y) is denoted by WL(x). Given a decision x, denote the probability density function or the probability mass function of L(x,Y) by fL(`; x,Y) where ` 2 WL(x) and the cumulative distribution function by FL(`; x,Y). Note 2.5.2 (Finite Loss Functions). For each loss function discussed in this thesis, EfL(x,Y)g < ¥ for all x 2 X. The assumption that the uncertainties modeled by Y are not affected by the choice of x is important (Rockafellar and Uryasev, 2002). For instance, in a financial setting, this assumption implies that stock prices fluctuate randomly regardless of who owns how much of each stock. Throughout this thesis the distribution of Y, the random vector of arc failures, is assumed to be unaffected by the choice of path x. In addition, it is assumed that arc failures are independent of each other. In particular, for each arc (ij) 2 A, Yij is a Bernoulli random variable, and Y : W ! f0, 1gjAj (see Section 3.1). Modeling of loss functions in the context of the SSPPPAF is discussed extensively in Section 3.3 of this thesis. The VaR(x, b) is the quantile associated with the threshold probability b of the distribution of L(x,Y) and is defined as follows: Definition 2.5.3 (Rockafellar and Uryasev (2002)). VaR(x, b) = min f`jFL(`; x,Y) bg (2.1) 13 Definition 2.5.4 (Rockafellar and Uryasev (2002)). Given a threshold percentile b, the cumulative upper btail distribution, denoted by Fb,L(`; x,Y), is given by: Fb,L(`; x,Y) = 8>>< >>: 0 if ` < VaR(x, b) FL(`;x,Y)b 1b if ` VaR(x, b) (2.2) Definition 2.5.5 (Rockafellar and Uryasev (2002)). For a given threshold percentile b, the CVaR(x, b) of the loss associated with decision x is given by: CVaR(x, b) = E fL(x,Y)jL(x,Y) VaR(x, b)g (2.3) In this thesis the random vector Y represents arc failures, so the distribution of Y is discrete. Hence, the loss functions induced are also discrete distributions. For loss functions with finite discrete distributions, the CVaR of the loss is defined below. Proposition 2.5.6. (Proposition 8, Rockafellar and Uryasev (2002)) Assume that the sample space of Y is finite, so that for each x 2 X the distribution of L(x,Y) is likewise discrete and finite and FL(`; x,Y) is a step function with jumps at ` 2 WL(x). For a fixed x, order ` 2 WL such that `0 < `1 < < `max. Let `k be the unique element of WL such that FL(`k1; x,Y) < b FL(`k; x,Y). (2.4) The VaR(x, b) of the loss is given by VaR(x, b) = `k. (2.5) In addition, the CVaR(x, b) is given by CVaR(x, b) = 1 1 b " (FL(`; x,Y) b) `k + `max å `=`k+1 ` fL(`; x,Y) # . (2.6) 14 2.6 Using CVaR in Optimization The minimization formulas below were developed by Rockafellar and Uryasev (2000, 2002) and rely on the convexity of CVaR. Theorem 2.6.1 (Rockafellar and Uryasev (2002)). Let Fb(x, z) = z + 1 1 b E n [L(x,Y) z]+ o , (2.7) where [z]+ = maxf0, zg. Fb(x, z) is finite and convex as a function of z 2 R and CVaR(x, b) = min z Fb(x, z). (2.8) Denote the set of minimizers of Fb(x, z) by Z. This set is a nonempty, closed interval, reducing to a single point if the minimum z of Fb(x, z) is unique, and VaR(x, b) = minfz 2 Zg. (2.9) In addition, CVaR(x, b) = Fb(VaR(x, b)). (2.10) Since Fb(x, z) is convex, it may be used in optimization problems. Specifically, the next theorem and corollary show that Fb(x, z) can be used to minimize CVaR(x, b) in the objective function or part of the constraints, which is how CVaR is used in the CVARSPP. Theorem 2.6.2 (Rockafellar and Uryasev (2002)). Minimizing CVaR(x, b) over x 2 X is equivalent to minimizing Fb(x, z) over X R since min x2X CVaR(x, b) = min (x,z)2X R Fb(x, z). (2.11) In addition, (x , z ) 2 X R are minimizers of Fb(x, z) (not necessarily unique) if and only if x 2 X is a minimizer of CVaR(x, b) and z is a minimizer of Fb(x , z). 15 Since Fb(x, z) is a convex piecewise linear function of z, it can be used as a constraint or objective function in an optimization problem. This is helpful since the direct formula in Proposition 2.5.6 requires computing VaR(x, b) in advance, which is usually difficult (Rockafellar and Uryasev, 2002). Using Fb(x, z) instead allows the use of CVaR in optimization without having to compute the VaR(x, b) in advance, as the corollary below mentions. The optimal solution found by minimizing CVaR(x, b) via the minimization formula Fb(x, z) will contain information about VaR(x, b) as a byproduct. Corollary 2.6.3 (Rockafellar and Uryasev (2002)). If (x , z ) minimizes Fb(x, z) over X R, then x minimizes CVaR(x, b) over X and CVaR(x , b) = Fb(x , z ) (2.12) and VaR(x , b) z , (2.13) where equality holds if the set of minimizers Z of Fb(x, z) (see Theorem 2.6.1) reduces to a single point. 16 CHAPTER 3 Model Development and Analysis 3.1 Modeling Arc Failures Throughout this thesis the simplifying assumption is made that failures are independent of each other. The probability that arc (ij) 2 A fails is pij where 0 pij < 1. Arcs for which pij = 0 never fail. The random vector of arc failures, Y, is of length jAj. Each component of Y is a random variable Yij. If arc (ij) fails for some ys 2 W, then ys ij = 0. In other words, Yij is a Bernoulli random variable. Yij Bernoulli(1 pij) 8 (ij) 2 A (3.1) 3.2 Formulation of the CVARSPP This model is an extension of the deterministic LP model of the shortest path problem in Section 1.2.3. It simply adds a constraint on CVaR to the LP model. When solved to optimality, a shortest path that satisfies the CVaR constraint will be chosen. Unfortunately, adding a CVaR constraint generally destroys the integrality of the shortest path polytope, so it is necessary to solve the problem as a MILP, greatly increasing solution time. 17 3.2.1 Additional Parameters In addition to the notation used in Formulation 1.2.2, the following notation is used throughout the rest of this thesis. Note that the set of arcs A is random. In the formulation, we denote by A the set of arcs such that pij < 1. C is the maximum acceptable CVaR of loss. b is the percentile for the upper tail (e.g. 0.95, 0.99) S W is the set of scenarios in the model, indexed by s. S may be a random sample from W or the whole sample space if W is small. ps is the probability of scenario s, normalized so that å s2S ps = 1. For equiprobable scenarios, set ps = 1 jSj 8 s 2 S. ys is the realization of the random vector Y corresponding to scenario s 2 S. 3.2.2 CVARSPP Optimization Model Formulation 3.2.1. Minimize å (ij)2A cijxij (3.2) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.3) CVaR(x, b) C (3.4) x 2 f0, 1gjAj (3.5) Constraint 3.4 can be rewritten in terms of the minimization formula Fb(x, z) from Theorem 2.6.1 as in the formulation below. 18 Additional Decision Variables z 2 R is VaR(x, b) at optimality. zs linearizes the excess loss [L(x, ys) z]+ in scenario s 2 S. Formulation 3.2.2. Minimize å (ij)2A cijxij (3.6) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.7) Fb(x, z) C (3.8) x 2 f0, 1gjAj (3.9) z 2 R (3.10) In Constraint 3.8, Fb(x, z) = z + 1 1 b å s2S ps [L(x, ys) z]+ . (3.11) Constraint 3.8 can be linearized using standard techniques and the distribution of L(x,Y) approximated by sampling to produce Formulation 3.2.3. As the loss functions are developed in the next section, the CVARSPP models for each loss function will be based on this formulation. Formulation 3.2.3. Minimize å (ij)2A cijxij (3.12) 19 subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.13) z + 1 1 b å s2S pszs C (3.14) zs L(x, ys) z 8 s 2 S (3.15) zs 0 8 s 2 S (3.16) x 2 f0, 1gjAj (3.17) z 2 R (3.18) 3.3 Loss Functions In this section three different loss function are introduced. Each loss function corresponds to a different effect of arc failures on the path. While the concepts behind these loss functions have appeared in the literature before, explicitly formulating and analyzing them as a part of a CVaR optimization model is the main contribution of this thesis. All the loss functions measure properties of an s t path x. The first loss function is called the Path Reliability Loss Function and measures the reliability of the path. The second function is the Arc Failure Loss Function, which measures the number arc failures. The final loss function is the Detours Loss Function, which counts the number of detours around broken segments of a path. 3.3.1 Path Reliability The first loss function to be considered derives from the basic physical interpretation of the problem. It asks “Does the path break in the current scenario?” In 20 other words, the loss is 0 if no arc fails on the current path and 1 otherwise. This loss function is worth looking at because it is simple and comparable to other SSPP approaches such as Ahuja (1988). This is the loss function implicit in the Most Reliable Path Problem, which finds the path with the lowest probability of failure (Ball et al., 1995b). In Figure 3.1 the path reliability loss function measures a loss of one. Figure 3.1: Path Reliability Loss Function: The loss for this path is 1. Definition 3.3.1 (Path Reliability Loss Function). Lr(x,Y) = max (ij)2A xij(1 Yij) (3.19) Remark 3.3.2. Lr(x,Y) = 0 if no arcs fail and 1 otherwise. Lr(x,Y) Bernoulli(q0) (3.20) where q0 = Õ (ij)2A(x) (1 pij) for all (ij) 2 A. Note that the path reliability loss function is piecewise linear and convex in x, so Formulation 3.2.3 can be further linearized to give the MILP below. Formulation 3.3.3 (Path Reliability Loss Function). Minimize å (ij)2A cijxij (3.21) subject to: 21 å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.22) z + 1 1 b å s2S pszs C (3.23) zs xij(1 ys ij) z 8 s 2 S (3.24) zs 0 8 s 2 S (3.25) x 2 f0, 1gjAj; z 2 R (3.26) 3.3.2 Number of Arc Failures The second loss function counts the number of arcs that fail on the path. For a decisionmaker who is responsible for repairing the path, this loss function gives an estimate of the amount of work to be done. In Figure 3.2, which is the same path as Figure 3.1, the arc failure loss function measures a loss of three. The drawback to this loss function is that not all arc failures, even on the path, may have an equal impact on the path length. For example, in reality it may be very easy to avoid one broken link. Figure 3.2: Arc Failures Loss Function: The loss for this path is 3. Definition 3.3.4 (Arc Failure Loss Function). La(x,Y) = å (ij)2A xij(1 Yij) (3.27) Remark 3.3.5. If pij = p for all (ij) 2 A, then La(x,Y) is the sum of Bernoulli random variables and La(x,Y) Bin(jA(x)j, 1 p) (3.28) 22 Since the arc failure loss function is linear, Formulation 3.2.3 reduces to the MILP below. Formulation 3.3.6 (Arc Failure Loss Function). Minimize å (ij)2A cijxij (3.29) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.30) z + 1 1 b å s2S pszs C (3.31) zs å (ij)2A xij(1 ys ij) z 8 s 2 S (3.32) x 2 f0, 1gjAj; z 2 R jSj + ; z 2 R (3.33) 3.3.3 Number of Detours The final loss function considered in this thesis counts the number of gaps in the current path. A gap is any set of consecutive failed arcs, such as arcs (2, 3) and (3, 4) in Figure 3.3. This loss function is related to ideas from Lagrangian relaxation and derives from the formulation of the problem. Mathematically the loss function measures total violation of flow conservation due to chance in the form of arc failures. This is a useful loss function to examine because it is common to measure the violation of the constraints in stochastic programming, particularly in terms of basic recourse. Figure 3.3 shows that on the same path as Figure 3.1 and Figure 3.2 the detours loss function measures a loss of 2. 23 Figure 3.3: Detours Loss Function: The loss for this path is 2. Definition 3.3.7 (Detours Loss Function). The loss function is Ld(x,Y) = 1 2 å i2N å j2N+(i) xijYij å j2N(i) xjiYji bi (3.34) where bs = 1, bt = 1 and bi = 0 for all i 2 Nnfs, tg. This loss function measures violation of the flow balance constraints or number of detours needed. The loss at some node i 2 N is Ld i (x,Y) = 1 2 å j2N+(i) xijYij å j2N(i) xjiYji bi , (3.35) and Ld(x,Y) = å i2N Ld i (x,Y). (3.36) Since flow conservation will only be violated at the ends of each gap in the path, this loss function ends up counting the number of gaps. The loss function can be interpreted physically as the number of detours necessary to avoid all the gaps. Since one failed arc renders the following arc inaccessible, the failure of that second arc does not cause additional loss. Until all the consecutive failed arcs have all been repaired, that gap in the path must still be detoured around. Lemma 3.3.8. For any incidence vector x of an s t path, Ld(x,Y) = å i2N(x) Ld i (x,Y) Proof. A feasible s t path x corresponds to a graph G(x) = (N(x), A(x)). Nodes i /2 N(x) have no arcs incident at them, and s, t 2 N(x). Hence, 24 Ld i (x,Y) = 1 2 å j2N+(i) xijYij å j2N(i) xjiYji bi 8 i /2 N(x) (3.37) = 1 2 å j2N+(i) 0 Yji å j2N(i) 0 Yij 0 (3.38) = 0 (3.39) So Ld(x,Y) = å i2N(x) Ld i (x,Y) (3.40) Lemma 3.3.9. If x is an incidence vector of an s t path, then Lds (x,Y) = 1 2 jYsi 1j where i 2 N+(s) \ N(x), (3.41) Ldj (x,Y) = 1 2 Yij Yjk 8 j 2 N(x)nfs, tg, (3.42) where i 2 N(j) \ N(x) and k 2 N+(j) \ N(x), and Ldt (x,Y) = 1 2 1 Yjt where j 2 N(t) \ N(x). (3.43) Proof. A feasible s t path x corresponds to a graph G(x) = (N(x), A(x)) such that one arc incident must be incident to s and t and two arcs to all other nodes i 2 N(x)nfs, tg. At s, N(s) = Æ and bs = 1. Also, xsi = 0 if i /2 N(x). Lds (x,Y) = 1 2 å j2N+(s) xsjYsj 0 bs (3.44) = 1 2 j1 Ysk + 0 1j for some k 2 N+(s) (3.45) = 1 2 jYsk 1j (3.46) 25 Similarly, at t, there exists some k 2 N(t) \ N(x) such that, Ld i (x,Y) = 1 2 j0 1 Ykt (1)j (3.47) = 1 2 j1 Yktj (3.48) For j 2 N(x)nfs, tg, let N(j) \ N(x) = fig, and N+(j) \ N(x) = fkg. Ldj (x,Y) = 1 2 å u2N(j) xujYuj å v2N+(j) xjvYjv bj (3.49) = 1 2 1 Yij 1 Yjk 0 (3.50) = 1 2 Yij Yjk (3.51) Corollary 3.3.10. For all j 2 N(x)nfs, tg, Ldj (x,Y) = 1 2 if Yij 6= Yjk; that is, if one arc fails and the other does not. If Yij = Yjk, then Ldj (x,Y) = 0. Lds (x,Y) = 12 if Ysi = 0, and similarly for t. Each gap consisting of consecutive nodes v1 . . . vk 2 N(x) results in a loss of kå i=1 Ldv i (x,Y) = 1. (3.52) Corollary 3.3.11. The sample space for Ld(x,Y) is Wd = f0, 1, . . . , j jN(x)j 2 k g Distribution of Detours Loss Function Theorem 3.3.12. Suppose that all arcs have an equal probability of failure; that is, pij = p for all (ij) 2 A. In that case, the pmf of the detours loss function Ld(x,Y) is given by fL(`; x,Y) = 8>>>>>>>>>>>< >>>>>>>>>>>: qjA(x)j when ` = 0, jA(x)j`+1 å i=` jA(x)ji å j=`1 piqjA(x)ji(i1 `1)(j1 `2)(jA(x)j i j + 1) when ` = 1 . . . j jN(x)j 2 k , and 0 otherwise, (3.53) where q = 1 p. 26 Proof. From Corollary 3.3.11, we know that there are between 0 and j jN(x)j 2 k gaps in the path subgraph G(x) = (N(x), A(x)). Let us then consider the number of ways to distribute the failed arcs into gaps and the surviving arcs into “fragments”. Let Ld(x,Y) = ` 2 Wd be the loss. Suppose ` > 0. Denote by i the number of failed arcs. By Corollary 3.3.10, the failed arcs must be distributed into exactly ` gaps with at least one arc in each gap, which implies that i `. The number of ways to distribute the i failed arcs into ` gaps with at least one arc in each gap is (i1 `1). Denote by j the total number of surviving arcs between the first and last failed arcs. Between each of the ` gaps are ` 1 “fragments” of at least one surviving arc each, which implies that ` 1 j jA(x)j i. The number of ways to distribute the j surviving arcs into ` 1 fragments with at least one arc in each fragment is (j1 `2). The length of the entire path segment from the first to the last failed arc is i + j jA(x)j. This segment may start anywhere in the first jA(x)j i j +1 arcs. The total number of ways to arrange i failed arcs with j jA(x)ji surviving arcs between the failed arcs into ` gaps with ` 1 fragments between them is i 1 ` 1 j 1 ` 2 (jA(x)j i j + 1). (3.54) The probability of i arcs failing and a loss of ` is then P(L(x,Y) = ` andå (ij)2A(x) Yij = i) = piqjA(x)ji jA(x)ji å j=`1 i 1 ` 1 j 1 ` 2 (jA(x)jij+1). (3.55) Hence the total probability of a loss of ` is P(L(x,Y) = `) = jA(x)j`+1 å i=` jA(x)ji å j=`1 piqjA(x)ji i 1 ` 1 j 1 ` 2 (jA(x)j i j + 1). (3.56) 27 The case of L(x,Y) = 0 is a special case. The only possible number of arc failures is i = 0, and all arcs must fall into a single fragment. The probability of ` = 0 is therefore fL(0; x,Y) = qjA(x)j (3.57) As was the case with the path reliability loss function, the detours loss function is piecewise linear and convex. Formulation 3.3.13 (Detours Loss Function). Minimize å (ij)2A cijxij (3.58) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (3.59) z + 1 1 b å s2S pszs C (3.60) zs å i2N `s i z 8 s 2 S (3.61) `s i å j2N(i) xjiysj i å j2N+(i) xijys ij bi (3.62) `s i bi + å j2N+(i) xijys ij å j2N(i) xjiysj i (3.63) zs 0 8 s 2 S (3.64) x 2 f0, 1gjAj; ` 2 RjSj RjNj; z 2 R (3.65) Table 3.1 summarizes the different loss functions. 28 Table 3.1: Summary of Loss Functions (LF) Name Formulation Sample Space Path Reliability LF Lr(x,Y) = max (ij)2A xij(1 Yij) f0, 1g Arc Failure LF La(x,Y) = å (ij)2A xij(1 Yij) f0, . . . jA(x)jg Detours LF Ld(x,Y) = 12 å i2N j å j2N(i) xjiYji å j2N+(i) xijYij bij f0, . . . j jN(x)j 2 k g 3.4 Comparison of Loss Functions 3.4.1 Magnitude of Losses As shown in Table 3.1, the sample space for the path reliability loss function is a subset of the sample space for the detours loss function, which is in turn a subset of the sample space for the arc failure loss function. Furthermore, the loss as measured by the arc failure loss function will always be greater than or equal to the loss measured by the detours loss function, which will in turn be greater than or equal to the loss measured by the path reliability loss function for the same realization ys 2 W. Theorem 3.4.1. For all feasible paths x and realizations of the random vector of arc failures ys 2 W, Lr(x, ys) Ld(x, ys) La(x, ys). (3.66) Proof. Consider an arbitrary feasible s t path, corresponding to the incidence vector x, and ys, a realization of Y. Let `r = Lr(x, ys), `d = Ld(x, ys), and `a = La(x, ys). If no two arcs fail consecutively, then each arc failure is its own gap in the path, and `d = `a. Otherwise at least two arcs will fail in a row, and `d < `a. Therefore Ld(x, ys) La(x, ys). 29 If `d = 0, no arcs have failed, so `r = 0 and `d = `r. If `d 1, then at least one arc has failed, and `r = 1. When `d 2, therefore, `d > `m. That is, Lr(x, ys) Ld(x, ys). 3.4.2 CVaR as a function of b Before beginning the next chapter, it is useful to examine how the choice of b affects CVaR(x, b). This section gives some guidance that will be used to check the results in the next chapter. Theorem 3.4.2 shows that CVaR(x, b) is piecewise hyperbolic in b but not necessarily convex in b. Theorem 3.4.2. Let L(x,Y) be a loss function with a finite discrete distribution. Given a fixed feasible path x, the corresponding loss distribution is fL(`; x,Y) with sample space WL(x) = f0, 1 . . . `maxg. For each ` 2 WL(x), CVaR(x, b) satisfies: CVaR(x, b) = ` + E n [L(x,Y) `]+ o 1 b (3.67) when P(L(x,Y) < `) < b P(L(x,Y) b). Proof. For each ` 2 WL(x), let b` = P(L(x,Y) `), let b1 = 0, and consider CVaR(x, b) for each b 2 (b`1, b`]. Applying Proposition 2.5.6 gives: CVaR(x, b) = 1 1 b " `(FL(`; x,Y) b) + `max å i=`+1 i fL(i; x,Y) # (3.68) = 1 1 b " `FL(`; x,Y) `b + `max å i=`+1 (i ` + `)fL(i; x,Y) # (3.69) = 1 1 b " `FL(`; x,Y) `b + ` `max å i=`+1 fL(i; x,Y) + `max å i=`+1 (i `)fL(i; x,Y) # (3.70) = 1 1 b " `(1 b) + `max å i=`+1 (i `)fL(i; x,Y) # (3.71) = ` + 1 1 b " `max å i=`+1 (i `)fL(i; x,Y) # (3.72) 30 CHAPTER 4 Computational Experiments The experiments in this chapter are exploratory, meant to illustrate the use and effects of CVaR. The purpose of these computational experiments are to explore the behavior of the three loss functions as C and b vary, to compare the loss functions to each other, and to illustrate the difference between using CVaR and expected value to chose a path though an uncertain network. 4.1 Implementation Details The experiments were performed on two HP xw4600 workstations with a quadcore 2.50 GHz Intel Core 2 CPU and 8 GB of RAM. The operating system on the first computer wasWindows XP and on the otherWindows Vista. The experiments were programmed using Python 2.7 (Van Rossum, 2003) and Gurobi 4.0.2 (Gu et al., 2010). The networks were implemented using the Networkx graph package (Hagberg et al., 2008). The Numpy matrix package (Oliphant, 2006) and the Scipy scientific computing package (Jones et al., 2001) were used for some calculations and for generating scenarios. Data was collected and stored in a Sqlite3 database using the SQLAlchemy objectrelationship manager package (ORM) (Bayer et al., 2011). Finally, results were analyzed and visualized with a mix of Microsoft Access, Microsoft Excel, and Matplotlib (Hunter, 2007), a Python plotting package. 31 4.2 Experiments Test instances were implemented as a directed network with source and sink attributes, a list of integer arc weights, a list of arc failure probabilities, and a set of scenarios. The network edges, arc weights, arc failure probabilities, and scenarios were all stored in the database using Python’s builtin pickle module. Since these experiments are exploratory, smaller instances were chosen. Preliminary experiments had found that instances of several hundred nodes took a long time to solve. Some of the larger and denser networks ran out of memory as well. The final test bed consisted of seven instances divided into three groups. The instances are summarized in Table 4.1. Table 4.1: Summary of Test Bed Name Nodes Arcs Density Scenarios Notes exp0002 6 9 30.0% 512 Complete enumeration branch02 27 30 4.27% 100 Only one arc per branch fails. branch03 27 30 4.27% 100 Less expensive paths are less reliable. branch04 27 30 4.27% 100 pij Uni f (0, 1), cij DU(1, 100) gnp01 25 59 9.83% 100 pij Uni f (0, 1), cij DU(1, 100) gnp02 25 138 23.0% 100 pij Uni f (0, 1), cij DU(1, 100) gnp03 25 243 40.5% 100 pij Uni f (0, 1), cij DU(1, 100) For each instance, the deterministic shortest path (DSP), the most reliable path (MRP, the path with the lowest probability of failure), and the most direct path (MDP, the path with the fewest arcs from source to sink) were calculated and recorded. All the sets of scenarios were generated as random samples of size 100 and treated as equiprobable except exp0002. All scenarios for that instance were 32 enumerated and their exact probability used in optimization. Descriptions of each component of a test instance were also stored in the database (see Appendix B for details). The main experiment mapped the level sets of the optimal cost as a function of b and C for each of the three loss functions. The full CVARSPP formulation for each loss function are in Chapter 3 and Appendix A. The parameter b was varied from 0 to 0.99 in increments of 0.01. The parameter C was varied from 0 to `max in increments of 0.01`max. The optimization status (infeasible, optimal, or unfinished), setup time, solve time, number of nodes, and number of integer feasible solutions found during optimization were recorded for each optimization run. If the run solved to optimality, the objective value, the values of each decision variable, and the arcs, cost, and reliability of the optimal path were recorded as well. For full details of what was recorded, see Appendix B. 4.2.1 Group 1: Complete Enumeration The purpose of the first test network was to enable direct comparison with the analytical results in Section 3.4. The CVARSPP models were still small enough to solve in <1s on average, despite containing all 512 scenarios. Since the network only contained 6 nodes and 9 arcs, it was possible to completely enumerate all st paths. The arc costs and failure probabilities for each path were chosen to force the MDP to be the MRP and have the highest cost and the DSP to be the least reliable. Figure 4.1 contains a diagram of the resulting test network. The lower edges of the level sets of the optimal path length are CVaR(x, b) for each path x that was optimal for the region. The values predicted by Theorem 3.4.2 are plotted over the level sets in Appendix C and agree with the optimal solutions found during the experiment. The level sets are piecewise hyperbolic but not convex in general. The plot for the detours loss function in Figure C.2(b) shows that 33 Figure 4.1: Network exp0002 t 1 2 s 3 4 1, .25 6, .1 2, .2 1, .25 4, .2 2, .25 3, .25 2, .2 2, .2 the lower edge of the MRP contour is flat at `max = 1 for b 0.675. The other two paths represented on the plot have three arcs, so the lower edge of their level sets keep increasing toward `max = 2. 4.2.2 Group 2: Branch Networks The effect of varying the arc failure probabilities on equallength paths was examined in a branch network consisting of a source and sink node connected by five internally nodedisjoint paths of five nodes each, as shown in Figure 4.2. In the first test network, branch02, only one arc was allowed to fail in each branch. As indicated by Theorem 3.4.1, all three loss functions gave the same results for each combination of b and C (Figures C.4(a), (b), and (c)). In the second test network, branch03, the arc costs and arc failure probabilities were set so that more reliable paths also had higher cost (Figure D.2). For instance, each arc on the MRP had cij = 5 and pij = 0.05, each arc on the next most reliable path had cij = 4 and pij = .1, and so on down to the DSP, where each arc had cij = 1 and pij = 0.25. The goal of this setup was to see if each path would be optimal in some region of the (b, C) plane. This was indeed the case under all 34 Figure 4.2: Branch Network Structure 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t three loss functions. In the final branch network, branch04, arc costs and arc failure probabilities were generated randomly for all arcs in order to observe how the CVARSPP models would respond to a wider range of arc reliabilities and path costs. The arc failure probabilities pij were drawn from a Unif(0, 1) distribution and the costs cij from a DU(1, 100) distribution. Since each branch is six arcs long, the expected reliability of each path is low. Even the MRP of branch04 failed in 98 out of 100 scenarios, giving an expected loss of 0.98 under the path failure loss function. When b > 0.98, the minimum choice of C that resulted in a feasible solution to the CVARSPP model under the path reliability loss function is C = 1 as shown in Figure C.8(a). The other two loss functions discriminated among paths, as E Ld(xMRP, ys) = 1.55 and E fLa(xMRP, ys)g = 2.05. This result suggests that the path reliability loss function is better suited for networks with highly reliable arcs or shorter st paths and should potentially be used with smaller b values. 4.2.3 Group 3: Uniform Random Networks The final group of test networks are uniform random networks (Erdös and Rényi, 1959) of 25 nodes with different densities. The arc costs and failure probabilities were randomly generated from pij Unif(0, 1) and cij DU(1, 100). The purpose 35 of this group of test networks was to observe the effect of increasing densities and the behavior of the CVARSPP models on less controlled test networks. As the density increased from 0.0983 for gnp01 to 0.230 for gnp02 to 0.405 for gnp03, the number of paths that were feasible in some region of the (b, C) plane increased as well, from 2 to 4 to 6 (Figures C.10(c), C.12(c), and C.14(c)). Network gnp02 also provided an illustration of the effects of sampling. The path corresponding to the region labeled 53 had two arcs, one of which failed in all 100 scenarios. That arc had a probability of failure of 0.96, so its failure in all scenarios is an artifact of sampling. This resulted in the path being infeasible for C < 1 under the path reliability loss function but being chosen under the detours and arc failure loss functions, as in Figure C.11(d). Under the detours loss function, though, CVaR(x53, b) = 1 for all b 1. If the scenarios had been fully enumerated, the path would be feasible in some narrow region under the path reliability loss function. 4.3 Comparison of Loss Functions To expand on Section 4.2.2 above, another feature of the the path reliability loss function is that, unless the network contains a path with no arcs that can fail, above b = 1 P(YMRP = 0) the instance will be infeasible for all C < 1. The CVARSPP model under the detours loss function has a feasible solution with C < `max for all b 1 in five out of seven networks. Under the arc failure loss function the CVARSPP model has a feasible solution with C < `max for all b 1 for all seven test networks. This suggests that b should be chosen somewhat lower and C somewhat higher within the range (0, `max), when using the path reliability loss function instead of the detours or arc failure loss functions. The detours loss function can use the same values of b as the arc failure loss function for most networks. Only in networks with higher arc failure probabilities does b need to be adjusted 36 slightly downward to preserve feasibility for C < `max. The CVaR limit C, on the other hand, should be chosen lower in the range (0, `max) than under the arc failure loss function, although precise guidance must wait on further computational experiments. Figure C.2 shows how the average optimal path length, average optimal path reliability, average CVaR loss varied as b increased for the test instances using network exp0002. For each loss function, the instances which solved to optimality were grouped by b. Figure C.3(d) records the percentage of instances in the experiment that were feasible for each b and loss function. Figures C.3(a)(c) record the average optimal path reliability, optimal path length, and optimal path CVaR for the instances that solved to optimality. The results all test networks are recorded in Appendix C. 37 CHAPTER 5 Conclusion and Future Work This thesis formulates CVaRbased models for the SSPPPAF. The literature contains models for the SSPPPAF using robust optimization and expected value objectives in stochastic programming (Chapter 2). There has also been some work attempting to trade off robustness and path length (see Ahuja (1988)). Although CVaR has recently been used for supplychain optimization (Madadi et al., 2010), to our knowledge this is the first attempt to use CVaR for a SSPP. The major contributions of this thesis are formulating the CVARSPP and developing three ways of modeling loss. The results found in Chapter 3 develop the loss functions. The three loss functions measure whether the path breaks, the number of gaps, and the number of arc failures with corresponding physical interpretations. Distributions for the each loss function were found for networks with uniform arc failure probabilities. The distribution of the path failure loss function was identified for any path with probabilistic arc failures. Moreover, an analytical characterization of the shape of CVaR(x, b) in terms of b was derived. Finally, it was shown that, for each possible scenario the arc failure loss function is always at least as large as the detours loss function, which is in turn always at least as large as the path reliability loss function. The lower edge of the level sets of the the fullyenumerated test network exp0002 agreed with the values calculated analytically. The computational experiments suggest that the path reliability loss function be used for networks with more reli 38 able paths and the tail threshold b should be set lower than when using the other loss functions. In general, the CVaRbased formulations find more reliable paths than when using expected value. This thesis has focused on modeling. Future research is needed to investigate computational issues. Preliminary experiments have indicated that scaling is likely to prove challenging. A decomposition approach such as the one proposed by KünziBay and Mayer (2006) may be needed to handle larger instances. Another scaling issue is the number of scenarios. Even the small networks in this study showed the effects of sampling. As the network grows, the number of scenarios needed to approximate the actual loss distribution increases rapidly. Some investigation into a better sampling scheme might also help with scaling. Another area for further investigation lies in extending singlestage CVaRbased formulations to twostage recourse models of the SSPPPAF. Finally, more investigation into the behavior of the CVaRbased formulation over a wider and more representative collection of networks is needed. It would be particularly interesting to test the models against some realworld networks, which include their own patterns of arc failures. Relaxing the assumption of independent arc failures and extending the analysis of the loss functions is another avenue to explore. 39 BIBLIOGRAPHY Acerbi, C. and Tasche, D. (2002). On the coherence of expected shortfall. Journal of Banking & Finance, 26(7):1487–1503. Ahuja, R. (1988). Minimum costreliability ratio path problem. Computers & Operations Research, 15(1):83–89. Ahuja, R. K., Magnanti, T. L., and Orlin, J. B. (1993). Network Flows: Theory, Algorithms, and Applications. Prentice Hall, 1 edition. Andreatta, G. and Romeo, L. (1988). Stochastic shortest paths with recourse. Networks, 18(3):193–204. Artzner, P., Delbaen, F., Eber, J.M., and Heath, D. (1999). Coherent Measures of Risk. Mathematical Finance, 9(3):203–228. Ball, M. O., Magnanti, T. L., Monma, C. L., and Nemhauser, G. L., editors (1995a). Network Models, volume 7 of Handbooks in Operations Research and Management Science. Elsevier Science. Ball, M. O., Monma, C. L., and Magnanti, T. L., editors (1995b). Network Routing, volume 8 of Handbooks in Operations Research and Management Science. Elsevier Science Pub Co. Bayer, M., Kirtland, J., de Menten, G., Trier, M., Jenvey, P., Aasma, A., Johnston, P., Ellis, J., and Others (2011). SQLAlchemy 0.7.1. software. Bellman, R. E. (1958). On a Routing Problem. Quarterly of Applied Mathematics, 16:87–90. 40 Bertsimas, D., Brown, D. B., and Caramanis, C. (2008). Theory and applications of robust optimization. Canada, J. R., Sullivan,W. G., Kulonda, D. J., and White, J. A. (2004). Capital Investment Analysis for Engineering and Management. Prentice Hall, 3 edition. Croucher, J. S. (1978). A note on the stochastic shortestroute problem. Naval Research Logistics, 25(4):729–732. Dantzig, G. B. (1957). DiscreteVariable Extremum Problems. Operations Research, 5(2). Deo, N. and Pang, C.Y. (1984). Shortestpath algorithms: Taxonomy and annotation. Networks, 14(2):275–323. Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1(1):269–271. Dreyfus, S. E. (1969). An Appraisal of Some ShortestPath Algorithms. Operations Research, 17(3). Erdös, P. and Rényi, A. (1959). On random graphs, I. Publicationes Mathematicae (Debrecen), 6:290–297. Floyd, R. W. (1962). Algorithm 97: Shortest path. Commun. ACM, 5(6):345+. Gu, Z., Rothberg, E., and Bixby, R. (2010). Gurobi 4.0.2. software. Hagberg, A. A., Schult, D. A., and Swart, P. J. (2008). Exploring network structure, dynamics, and function using NetworkX. In Proceedings of the 7th Python in Science Conference (SciPy2008), pages 11–15, Pasadena, CA USA. Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science and Engineering, 9(3):90–95. 41 Jaillet, P. (1992). Shortest path problems with node failures. Networks, 22(6):589– 605. Jones, E., Oliphant, T., Peterson, P., and Others (2001). SciPy: Open source scientific tools for Python. software. Kendall, K. E. and Kendall, J. E. (2008). Systems analysis and design. Pearson/ Prentice Hall. KünziBay, A. and Mayer, J. (2006). Computational aspects of minimizing conditional valueatrisk. Computational Management Science, 3(1):3–27. Madadi, A. R., Kurz, M. B., Taafe, K., Mason, S., Pohl, E., Root, S., and Sir, M. (2010). Managing Disruptions in Pharmaceutical Supply Chain Networks. In Johnson, A. and Miller, J., editors, Industrial Engineering Research Conference. Oliphant, T. E. (2006). Guide to Numpy. Tregol Publishing, Provo, UT. Papadimitriou, C. (1991). Shortest paths without a map. Theoretical Computer Science, 84(1):127–150. Pollack, M. and Wiebenson, W. (1960). Solutions of the ShortestRoute ProblemA Review. Operations Research, 8(2). Rockafellar, R. and Uryasev, S. (2002). Conditional valueatrisk for general loss distributions. Journal of Banking & Finance, 26(7):1443–1471. Rockafellar, R. T. and Uryasev, S. (2000). Optimization of conditional valueatrisk. Journal of Risk, 2:21–41. Sen, S. and Higle, J. L. (1999). An Introductory Tutorial on Stochastic Linear Programming Models. Interfaces, 29(2). 42 Van Rossum, G. (2003). The Python Language Reference Manual. Network Theory Ltd. Yu, G. and Yang, J. (1998). On the Robust Shortest Path Problem. Computers & Operations Research, 25(6):457–468. 43 APPENDIX A Optimization Models A.1 Parameters These parameters are used in all three optimization models. Graph G = (N, A) is a directed network. S W is the set of scenarios in the model, indexed by s. S may be a random sample from W or the whole sample space if W is small. cij is the cost of arc (ij) 2 A. ys is the realization of the random vector Y corresponding to scenario s 2 S. ps is the probability of scenario s, normalized so that Så s=1 ps = 1. C is the maximum acceptable CVaR loss. b is the percentile for the upper tail (e. g. 0.50, 0.95, 0.99). A.2 Decision Variables xij indicates which arcs are chosen for the shortest path. xij = 1 if arc (ij) is chosen, and 0 otherwise. z 2 R is VaR(x, b) at optimality. zs linearizes the excess loss [L(x, ys) z]+ in scenario s 2 S. `s i is the loss at node i 2 N in scenario s 2 S and linearizes the loss function. 44 A.3 Path Reliability Loss Function Lr(x,Y) = max (ij)2A xij(1 Yij) (A.1) Formulation A.3.1 (Path Reliability Loss Function). Minimize å (ij)2A cijxij (A.2) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (A.3) z + 1 1 b å s2S pszs C (A.4) zs xij(1 ys ij) z 8 s 2 S, (ij) 2 A (A.5) x 2 f0, 1gjAj; z 2 R jSj + ; z 2 R (A.6) 45 A.4 Arc Failure Loss Function La(x,Y) = å (ij)2A xij(1 Yij) (A.7) Formulation A.4.1 (Arc Failure Loss Function). Minimize å (ij)2A cijxij (A.8) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (A.9) z + 1 1 b å s2S pszs C (A.10) zs å (ij)2A xij(1 ys ij) z 8 s 2 S (A.11) x 2 f0, 1gjAj; z 2 R jSj + ; z 2 R (A.12) 46 A.5 Detours Loss Function Ld(x,Y) = å i2N å j2N(i) xjiYji å j2N+(i) xijYij bi (A.13) Formulation A.5.1 (Detours Loss Function). Minimize å (ij)2A cijxij (A.14) subject to: å j2N+(i) xij å j2N(i) xji = 8>>>>>>< >>>>>>: 1 if i = s 1 if i = t 0 if i 2 Nnfs, tg (A.15) z + 1 1 b å s2S pszs C (A.16) zs å i2N `s i z 8 s 2 S (A.17) `s i å j2N(i) xjiysj i å j2N+(i) xijys ij bi 8 s 2 S, i 2 N (A.18) `s i bi + å j2N+(i) xijys ij å j2N(i) xjiysj i 8 s 2 S, i 2 N (A.19) x 2 f0, 1gjAj; z 2 R jSj + ; ` 2 RjSj RjNj; z 2 R (A.20) 47 APPENDIX B Schema of Experiment Database This appendix describes the database used in the computational experiments. The database contains all data generated during the experiments. The database schema is diagrammed in Figure B.2 using crow’s foot notation. Each table in the database and its fields are described below. The database was normalized to the Backus Naur Third Normal Form (Kendall and Kendall, 2008) before being selectively denormalized for ease of use. B.1 Database Tables Graph The table for the directed networks. This does not contain arc weights, arc failure probabilities, or scenarios. name The name of the directed network. descr A sentence or two describing the network to a human reader. n The number of nodes. a The number of arcs. source The source node of the network. sink The sink node of the network. density The density of the network obj The Python object containing the arc list. Weights The table for arc weights (cij). 48 name A short name for the arc weights. descr Asentence or two describing the arc weights or their generation method to a human reader. obj The Python object containing the list of arc weights. length The number of weights in the list of arc weights. graph_id Foreign key to Graphs table. AFPs The table for arc failure probabilities (pij). name A short name for the arc failure probabilities. descr A sentence or two describing the arc failure probabilities and their generation method to a human reader. obj The Python object containing the arc failure probabilities. length The number of arc failure probabilities in the list. graph_id Foreign key to Graphs table. Scenarios The table for scenarios (ys and ps for all s 2 S) name A short name for this set of scenarios. descr Asentence or two describing the scenarios and their generation method to a human reader. s The number of scenarios. scen_array The Python object containing an jAj s Numpy matrix of the scenarios. Each row is ys, the incidence vector of arc failures. ps The Python object containing a Numpy vector of scenario probabilities. Each element is corresponds to ps. empafp The Python object containing a Numpy vector of jAj entries. Each entry is the percent of scenarios in which that arc fails. 49 afp_id Foreign key to AFPs table. TestInstance The table for instances. This is a full network scenarios and with weights and failure probabilities attached to each arc. name A short name for this instance. graph_id Foreign key to Graphs table. wts_id Foreign key to Weights table. scen_id Foreign key to Scenarios table. Path The table for information about paths of interest through a network. rel The apriori reliability (i.e. Õ (ij)2A(x) (1 pij)). cost The cost of the path (i.e. å (ij)2A(x) cij). arcs Number of arcs in this path. type Classification—DSP, MRP,MDP, or an optimal solution to some CVARSPP instance. inst_id Foreign key to the TestInstance table. sol_id Optional: Foreign key to the Solution table. PathArcs The table that stores the arcs in a path. head The head of an arc in a path. tail The tail of an arc in a path. path_id Foreign key the the Path table indicating which path contains these arcs. Exp The table of summary information about a computational experiment. The information in this table is mostly intended as notes for a human reader. 50 name A short name for this set of experiment. descr A summary of the experiment. goal The objective of the experiment. beta_method A summary of how the tail threshold probabilities b are to be generated. c_method A summary of how the maximum acceptable losses C are to be generated. loss_fns The loss functions involved in this experiment. computer The machine used to run the experiment. solver The program used to solve the CVARSPP models. Param The table storing solver control parameters such as time limit. filename The name of the file with the parameters for the solver to read. folder The path to the file. suffix The extension of the file indicates format. plist The Python object containing the parameters, a duplicate of the file. Run Each record represents a unique instance of the CVARSPP. name A short name for this CVARSPP instance. beta The tail threshold probability b. C The Python object representing the arc failure probabilities (a BLOB in SQL terms). loss_fn Which loss function to use; i.e., which specific model to choose (see Appendix A). solved Flag indicating whether the run has been sent to the solver or not. 51 exp_id Foreign key to Exp table. param_id Foreign key to the Param table. inst_id Foreign key to TestInstance table. OptInfo The table for information reported by the solver about an instance. If an optimal solution is found, it is stored in the Solution table. optimal Did the instance solve to optimality? grb_status The status code reported by the solver. setup The time to set up each instance. sol_time The solution time for each instance. iter The number of simplex iterations performed for all the nodes explored. nodes The number of branchandbound nodes explored. best_bound If the solver ran out of time, the best bound found during optimization. num_sol The number of integer feasible solutions during optimization. obj The Python object containing all the information about running the optimization which was returned by the solver. date Timestamp for the optimization run. run_id Foreign key for the Run table. Solution The table for information about the optimal solution to a run. obj_cost The cost of the optimal path. var The VaR for this instance. cvar The CVaR for the instance, calculated C and the slack in the CVaR constraint. 52 cvar_slack The slack in the CVaR constraint. dvs The Python object containing the values of all decision variables. run_id Foreign key for the Run table. B.2 Database Schema Figure B.1: Legend: Crow’s Foot Notation Primary Key Attribute (FK): Foreign Key (O): Optional Table Name Relationship Symbol Zero or onetozero or one Zero or onetomany Onetozero or many Onetomany 53 Figure B.2: Database Schema for Computational Experiments after Kendall and Kendall (2008) Graph id name descr n a source sink density obj Weights id name descr obj length Scenarios id name descr s scen_array ps empafp afp_id (FK) AFPs id name descr obj length TestInstance id name graph_id (FK) wts_id (FK) scen_id (FK) belongs to / is constructed from is used to generate / is generated from is used in / uses Exp id name descr goal beta_method c_method loss_fns computer solver Run id name beta C loss_fn solved exp_id (FK) param_id (FK) inst_id (FK) obj_cost var cvar cvar_slack dvs run_id (FK) OptInfo id optimal grb_status setup sol_time iter nodes best_bound num_sol obj date run_id (FK) Param id filename (O) folder (O) suffix (O) plist (O) contains / is optimal for Path id rel cost arcs type inst_id (FK) sol_id (O) (FK) PathArc id head tail path_id (FK) is comprised of / are on is solved during / solves the SSPPAF on has / belongs to of is used by / passes these to Gurobi records optmization information in / is reported by Gurobi after finds / is extracted from 54 APPENDIX C Objective Function Contour Graphs by Loss Function This appendix contains graphs of the objective function values for each test instance and loss function. 55 Figure C.1: exp0002: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C CVaR(MRP, b) CVaR(x6, b) CVaR(DSP, b) (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C 56 Figure C.2: exp0002: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 0 2 4 6 8 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 1 1.5 2 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 57 Figure C.3: branch02: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C 58 Figure C.4: branch02: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 5 102 0.1 0.15 0.2 0.25 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 240 260 280 300 320 b Average Length (c) Average CVaR for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.4 0.6 0.8 b Average CVaR (d) Infeasible instances as a function of b 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 59 Figure C.5: branch03: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 2 4 6 b Maximum Loss C 60 Figure C.6: branch03: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 10 15 20 25 30 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 1 1.5 2 2.5 3 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 61 Figure C.7: branch04: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 2 4 6 b Maximum Loss C 62 Figure C.8: branch04: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 200 250 300 350 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 1 2 3 4 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 63 Figure C.9: gnp01: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 2 4 6 b Maximum Loss C 64 Figure C.10: gnp01: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 52.5 53 53.5 54 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 0.5 1 1.5 2 2.5 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 65 Figure C.11: gnp02: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 b Maximum Loss C (c) Arc failure loss function 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 b Maximum Loss C 66 Figure C.12: gnp02: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.2 0.4 0.6 0.8 50 55 60 65 70 b Average Length (c) Average CVaR for feasible instances 0 0.2 0.4 0.6 0.8 0.5 1 1.5 2 b Average CVaR (d) Infeasible instances as a function of b 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 67 Figure C.13: gnp03: Level sets of the optimal path length as a function of b and C (a) Path reliability loss function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8 1 b Maximum Loss C (b) Detours loss function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.5 1 1.5 2 b Maximum Loss C (c) Arc failure loss function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 1 2 3 b Maximum Loss C 68 Figure C.14: gnp03: Change in properties of the average optimal path vs. b (a) Average reliability for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8 b Average Reliability (b) Average optimal path length for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 60 80 100 120 140 160 b Average Length (c) Average CVaR for feasible instances 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.4 0.6 0.8 1 1.2 1.4 b Average CVaR (d) Infeasible instances as a function of b 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.4 0.6 b Percent Infeasible Lr(x,Y) Ld(x,Y) La(x,Y) 69 APPENDIX D Test Network Diagrams D.1 Small Network: exp0002 The scenarios for this small network are completely enumerated. Arc failure probabilities are chosen to force the MDP to be the MRP and have the highest cost while the DSP is the least reliable path. t 1 2 s 3 4 1, .25 6, .1 2, .2 1, .25 4, .2 2, .25 3, .25 2, .2 2, .2 Figure D.1: Network exp0002 70 D.2 Branch Networks The three following networks have five disjoint paths, or branches, of six arcs connecting the source s to sink t. Figure D.2: Network branch02 The arc weights cij are drawn from DU(1, 100). One arc is allowed to fail on each branch with probability pij drawn from Unif(0, 1). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t 28 86 20 51 15 98 18 74 58 39 31 84 94 31 54 64 70 55 32 87 22 49 41 55 43, .12 8 11, .22 33, .68 8, .12 85, .20 Figure D.3: Network branch03 The arc failure probabilities and arc costs are uniform across each branch. As the cost of the path increases from the DSP to the most expensive path, so does the reliability. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t 1, .25 2, .2 3, .15 4, .1 5, .05 1, .25 2, .2 3, .15 4, .1 5, .05 1, .25 1, .25 1, .25 1, .25 2, .2 2, .2 2, .2 2, .2 3, .15 3, .15 3, .15 3, .15 4, .1 4, .1 4, .1 4, .1 5, .05 5, .05 5, .05 5, .05 71 Figure D.4: Network branch04 The arc failure probabilities pij are drawn from Unif(0, 1), and the arc weights cij are drawn from DU(1, 100). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 s t 31, .25 2, .47 35, .47 53, .25 91, .37 48, .8 53, .13 13, .99 84, .64 34, .67 17, .16 18, .4 67, .62 26, .96 53, .087 87, .94 74, .38 32, .072 96, .4 72, .65 67, .39 56, .34 54, .026 93, .66 15, .75 35, .29 60, .91 79, .44 88, .71 14, .032 D.3 Uniform Random Networks The final three networks are GNP(n, p) random networks (Erdös and Rényi, 1959). All three networks have n = 25 nodes. The arc failure probabilities pij are drawn from Unif(0, 1), and arc weights cij from DU(1, 100). Figure D.5: Network gnp01 — GNP(25, 0.1) t 12 11 10 9 6 7 8 5 4 3 2 0 s 13 14 15 16 17 18 19 20 21 22 23 72 Figure D.6: Network gnp02 — GNP(25, 0.25) t 12 11 10 9 6 7 8 5 4 3 2 1 s 13 14 15 16 17 18 19 20 21 22 23 Figure D.7: Network gnp03 — GNP(25, 0.4) t 12 11 10 9 6 7 8 5 4 3 2 1 s 13 14 15 16 17 18 19 20 21 22 23 73 VITA Juliana Bright Candidate for the Degree of Master of Science Thesis: ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK Major Field: Industrial Engineering and Management Biographical: Personal Data: Born in London, UK on April 22, 1980. Education: Bachelor of Arts in German from the University of Oklahoma, December 2005. Bachelor of Science in Industrial Engineering and Management from Oklahoma State University in Stillwater, Oklahoma in August 2009. Completed the Requirements for the Master of Science degree at Oklahoma State University in December 2011 with a major in Industrial Engineering and Management. Experience: Graduate research assistant, Industrial Engineering and Management, Oklahoma State University, Fall 2010Spring 2011; undergraduate and graduate teaching assistant, CEATLabs, Oklahoma State University, January 2008August 2010; staff assistant, DollarThrifty Automotive Group, Tulsa, OK, June 2004September 2006. Professional Memberships: Engineering Intern, Institute for Operations Research and Management Science (INFORMS), Institute of Industrial Engineers (IIE) Awards and Honors: Doctoral Academy Fellowship, University of Arkansas; First place in the 2010 IIE Undergraduate Student Technical Paper Competition. Name: Julie Anna Bright Date of Degree: December 2011 Institution: Oklahoma State University Location: Stillwater, Oklahoma Title of Study: ROBUST SHORTEST PATHS UNDER UNCERTAINTY USING CONDITIONAL VALUEATRISK Pages in Study: 73 Candidate for the Degree of Master of Science Major Field: Industrial Engineering and Management Finding a shortest path in a network is a classical problem in discrete optimization. The systems underlying the network models are subjects to a variety of sources of uncertainty. The purpose of this thesis is to model the shortest path problem with probabilistic arc failures with the aim of finding short but reliable paths through the network. This thesis proposes to use Conditional ValueatRisk (CVaR) to find a path that is robust under probabilistic arc failures. CVaR, a quantitative risk measure, is roughly the mean excess loss associated with a decision. This thesis also develops three models of losses due to arc failures. Mixed integer linear programming models for the stochastic shortest path problem with probabilistic arc failures are formulated and implemented. The optimization models limit the CVaR of the loss due to arc failures asmeasured by the models of loss developed in the thesis. ADVISOR’S APPROVAL: Dr. Balabhaskar Balasundaram 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


