New algorithmic framework for conditional value at risk: Application to stochastic fixed-charge transportation

New algorithmic framework for conditional value at risk: application to stochastic fixed-charge transportation Elena Fern´andeza, Yolanda Hinojosab, Justo Puertoc, Francisco Saldanha-da-Gamad aDepartamento d’Estad´ıstica i Investigaci´o Operativa, Universitat Polit`ecnica de Catalunya, Barcelona, Spain bIMUS and Departamento de Econom´ıa Aplicada I, Universidad de Sevilla, Spain cIMUS and Departamento Estad´ıstica e Investigaci´on Operativa, Facultad de Matem´aticas, Universidad de Sevilla, Spain dDepartamento de Estat´ıstica e Investiga¸c˜ao operacional and Centro de Matem´atica, Aplica¸c˜oes Fundamentais e Investiga¸c˜ao Operacional, Faculdade de Ciˆencias, Universidade de Lisboa, Portugal Abstract This paper introduces a new algorithmic scheme for two-stage stochastic mixed-integer programming assuming a risk averse decision maker. The focus is the minimization of the conditional value at risk for a hard combinatorial optimization problem. Some properties of a mixed-integer non-linear programming formulation for conditional value at risk are studied as well as their algorithmic implications. This yields to a procedure for obtaining lower and upper bounds on the optimal value of the problem that may lead to an optimal solution. The new developments are applied to a fixed-charge transportation problem with stochastic demand, and they are computationally tested. The corresponding results are thoroughly presented and discussed. Keywords: Transportation, Stochastic mixed-integer programming, CVaR 1. Introduction In this work we propose an algorithmic scheme for two-stage stochastic mixed-integer programming assuming a risk averse decision maker. As recently pointed out by Escudero et al. (2017), many risk measures can be found in the literature, among which the most popular ones are scenario immunization, semi-deviation, min-risk, stochastic dominance, value-at-risk, and the conditional value-at-risk (CVaR) (see e.g., Chelst and Canbolat, 2011, for further details on these measures). The two latter measures are possibly the most widely used in the literature. The risk measure adopted in this work is the CVaR (Uryasev and Rockafellar, 2001). Stochastic Programming has become a very important tool to hedge against uncertainty when compared to approaches derived from simplified (deterministic) counterpart models (Birge and Louveaux 2011, Wallace and Ziemba 2005). Due to applications arising in many relevant fields (e.g., logistics), an increasing amount of work is being done on two-stage, and Preprint submitted to Elsevier January 29, 2024 multi-stage, stochastic mixed-integer linear programming (e.g., Santoso et al., 2005; Hinojosa et al., 2014; Nickel et al., 2012). Initially, stochastic programming literature focused predominantly on risk-neutral decision makers, i.e., the future assets are accounted for using a mean-type objective function. More recently this trend has changed and there exists an increasing interest on models capturing a risk-averse attitude. To the best of our knowledge, the first work considering a mean-risk objective function in two-stage stochastic programming is Ahmed (2006), where the complexity of mean-risk stochastic programming considering different risk measures was investigated. F´abi´an (2008) studied two-stage CVaR minimization problems (as well as CVaR-constrained problems) and devised decomposition and solution schemes based on the linear programming representation of CVaR proposed by K¨unzi-Bay and Mayer (2006). Also in the context of two-stage stochastic linear programming problems, Miller and Ruszczy´nski (2011) consider an objective function combining several conditional risk measures. For the case of two-stage stochastic mixed-integer linear programming problems, Schultz and Tiedemann (2006) presented an explicit mixed-integer linear programming formulation for the CVaR, when the probability distribution is discrete and finite. For the same problem Carøe and Schultz (1999) devised a dual decomposition approach by relaxing the non-anticipativity constraints. The use of risk measures in two-stage stochastic mixed-integer linear programming problems has been considered in several applications. Noyan (2012) proposed a CVaR mean-risk model in the context of disaster relief, which was solved with a variation of the L-Shaped method (Van Slyke and Wets, 1969). In the context of retrofitting infrastructures in a transport network to protect against disaster events (e.g., earthquake), Liu et al. (2009) studied a mean-risk two-stage stochastic programming model in which risk is handled via a semideviation measure. Recently, Homem-de-Mello and Pagnoncelli (2016) discussed risk measures in multi-stage stochastic programming and potential applications. Also, Escudero et al. (2017) considered a mean-risk model using a so-called time stochastic dominance measure and proposed a scenario cluster Lagrangean decomposition approach. Finally, we refer to the work by Alonso-Ayuso et al. (2016) who presented a mean-risk multistage stochastic mixed-integer linear programming model for a forestry management problem. In this work we focus on the minimization of the CVaR for a hard combinatorial optimization problem. We investigate a two-stage stochastic programming formulation, derived from the direct definition of CVaR, as an alternative to previous linear representations for the case when the random vector has a discrete distribution (Rockafellar and Uryasev, 2002; Schultz and Tiedemann, 2006). In particular, we consider a mixed-integer non-linear programming formulation and study some of its properties as well as their algorithmic implications. Given 2 the combinatorial nature of the considered problem, the evaluation of each scenario requires solving a hard optimization problem, even if the number of scenarios is discrete. However, the properties of the proposed formulation can be exploited to derive a solution framework, as an alternative option to schemes based on linear representations. The new algorithmic framework has been applied to a risk-averse stochastic fixed-charge transportation problem, which generalizes the risk-neutral objective studied in Hinojosa et al. (2014). This problem incorporates three main features to a standard transportation problem. First, it is assumed that a strategic decision has to be made concerning the distribution channels that will be activated for shipping demand between origin/destination pairs. Similarly to Balinski (1961), activating a distribution channel incurs a fixed cost. Second, variable costs associated with flows, both for handling them at the origins and for sending them through the activated network, are considered as well. This extends the fixed-charge transportation problem of Balinski (1961), which does not consider handling costs at the origins. Third, we assume that demand is stochastic (reflecting volatility of current markets leading to continuously changing demands), and consider the problem as a two-stage stochastic problem with recourse. The first stage is given by the strategic decision concerning the distribution channels and the recourse action refers to the tactical decision corresponding to the transportation pattern. Examples of potential applications of such a problem are listed in Hinojosa et al. (2014) and arise in logistics systems in which previous contracts are needed in order to activate distribution channels (see, for instance, Xu and Nozick, 2009). Other examples can be found in telecommunications where traffic can be sent through some channel only if such channel has been previously made physically available (see Ahuja et al., 1993). Applications can also be found in production planning when a set of machines is available during a certain period of time for producing a set of products and one has to decide which machines will produce which products and in what quantities (see, for instance, Grieco et al., 2001, and the references therein). Finally, one potential area where the findings of this work can be applied concerns fleet management. In fact, most of the distribution companies must decide in advance about the fleet portfolio to use (type and number of vehicles) before knowing the demand. Some work in this direction includes List et al. (2003) and Meng et al. (2012). It is also possible to find some work in the literature that investigates the use of CVaR in the context of logistics/transportation systems. Hemmati et al. (2016) highlights the usage of CVaR as a risk measure in a context of energy transportation. The authors propose a two-stage stochastic mixed integer linear programming. Uncertainty is associated with wind resources, i.e., it is associated with supply. The decisions are related with scheduling storage units. Toumazis and Kwon (2016) consider CVaR for routing decisions in the context of hazmat materials transportation. In this case, uncertainty is associated with the occur- 3 ring accidents and their consequences. The concept of worst case CVaR is introduced as a way to hedge against data inaccuracy. A robust optimization modeling framework is adopted. In this paper we develop an algorithmic framework for a two-stage stochastic combinatorial optimization problem, with a finite set of scenarios, assuming a CVaR objective. The solution algorithm exploits some properties of a mixed-integer non-linear programming formulation, derived from the direct definition of CVaR. Alternative solution strategies are applied and computationally tested for a fixed-charge stochastic transportation problem that generalizes the problem studied in Hinojosa et al. (2014). The obtained numerical results are thoroughly discussed and analyzed. These results show that our proposal outperforms alternative methods for the largest tested instances. The remainder of this paper is organized as follows. Section 2 recalls the concepts of CVaR on risk-averse decision making as well as its computation via the mathematical programming representation developed in Uryasev (2000), Rockafellar and Uryasev (2002) and Schultz and Tiedemann (2006). In Section 3 we introduce the mixed integer non-linear formulation for CVaR that we use, we study some of its properties, and present the proposed solution framework. Section 4 focuses on its application to the CVaR stochastic fixed-charge transportation problem. The numerical results of the computational tests executed with our solution proposal are presented and analyzed in Section 5. The paper ends in Section 6 with some conclusions and guidelines for future research. 2. Models for risk-averse decision making Consider a two-stage stochastic problem, where stochasticity is associated with a random vector ξand is expressed via a finite set of possible scenarios Ω. Let πω,ω∈Ω, denote the probabilities of the different scenarios. In two-stage stochastic problems an a priori (firststage) solution must be defined before the actual realization of ξis known. We denote by Y the domain for a priori solutions. The second-stage recourse action associated with a given first-stage decision y∈Yrepresents the set of decisions to be made for each possible realization of the uncertain vector ξ. Let R(y,ξω) denote the cost of the recourse action under scenario ω∈Ω. Computing R(y,ξω) requires to solve some associated optimization problem on the set of second-level variables. Moreover, when an a priori solution y∈Yis considered, it is not possible to compute the actual cost of the recourse action. The reason is that when the a priori solution y∈Yis chosen, the specific scenario ω∈Ω that will actually occur in the realization of ξis unknown. In this context it is natural to consider some objective function that incorporates joint information of all possible realizations. One such objective for the case of risk-averse decision making is the value-at-risk (VaR) that we revisit below. 4 Let R(y;ξ) = f(y) + R(y,ξ) be the cost associated with an a priori solution y∈Y, where f(y) represents the here-and-now cost, and R(y,ξ) the value of the recourse function. The latter is a random variable representing the cost of the optimal second-stage decision for the a priori solution and thus, R(y;ξ) is also a random variable. For a given confidence level α, the VaR is defined as the α-quantile of the probability distribution of R(y;ξ). That is, the VaR of an a priori solution y∈Y, denoted by ηα(y), is the value of the minimum threshold ηfor which P[R(y;ξ)≤η] is at least α. More formally, ηα(y) = min {η|P[R(y;ξ)≤η]≥α}. Using as objective the minimization of the VaR for a given confidence level α, that is, miny∈Yηα(y), the decision maker aims at minimizing the threshold for the risk level 1 −α. Still, for the solution yobtained with this objective, the value of R(y;ξ) may be too high, for some of the scenarios whose cost exceeds this minimum threshold, even if their joint probability is very low (smaller than or equal to 1−α). An alternative that takes into account this consideration is to minimize the expected cost of the scenarios whose cost exceeds the threshold for a fixed confidence level α. This leads to the minimization of the α−Conditional Value at Risk (α-CVaR), which is precisely the objective that we consider in this work and we formally define below. The α-CVaR of an a priori solution y∈Y, Ψα(y), is the expected value of the total cost in the (1 −α)×100% worse scenarios for the a priori solution y, i.e., it is the expected cost conditional to the scenarios whose value exceeds ηα(y). Note that the α-CVaR generalizes the expected cost objective, which is a particular case when all the scenarios come into play. Thus, when αis suitably chosen the α-CVaR of an a priori solution yis just the expected value of R(y;ξ) (see Puerto et al., 2017). For a given solution y∈Y, it is difficult to express Ψα(y), unless an analytical representation for ηα(y) is known, because the latter is involved in the definition of the former. One possibility (see Uryasev, 2000; Rockafellar and Uryasev, 2002, for further details) is to consider the following function: Φα(y, η) = η+1 1−αE[(R(y;ξ)−η)| R(y;ξ)> η], which, when uncertainty is expressed by a finite set of scenarios reduces to Φα(y, η) = η+1 1−αX ω∈Ω (R(y;ξω)−η)+πω,(1) with a+= max{0, a}. Then, the α-CVaR of y∈Ycan be computed as: Ψα(y)=Φα(y, η(y)) = min ηΦα(y, η). 5 Thus, the problem of finding a feasible vector y∈Ywith the smallest α-CVaR reduces to: CVaRαmin y∈YΨα(y) = min y∈YΦα(y, η(y)) = min y∈Y, η,Φα(y, η),(2) which, taking into account (1), can be expressed as: CVaRαminimize η+1 1−αX ω∈Ω πωψω(3) subject to y∈Y, (4) ψω≥ R(y;ξω)−η, ω ∈Ω,(5) ψω≥0, ω ∈Ω.(6) Note that the meaning of the ψωvariables in optimal solutions to the above formulation is precisely ψω= (R(y;ξω)−η)+,ω∈Ω. When the domain Ycan be expressed by means of a set of linear constraints, the above formulation is a valid mathematical programming formulation, linear on the decision variables y. It is possible to derive an alternative formulation for CVaR by observing that when y∈Y is fixed, its α-CVaR, Ψα(y), can be expressed just in terms of y(Rockafellar and Uryasev, 2002). In order to make the paper self-contained we include the reformulation in which we base our algorithmic framework. For a given a priori solution y∈Y, let us solve the associated second-stage problem for each realization, ξωω∈Ω, of the uncertain vector ξ, then sort the scenarios by non-increasing values of R(y;ξω), and finally index them accordingly. To alleviate notation, if ωr(y) is the index of the r-th sorted scenario we will write ξr(y)and πr(y)instead of ξωr(y)and πωr(y), respectively. Hence, R(y;ξ1(y))≥. . . ≥ R(y;ξ|Ω|(y)), with probabilities π1(y), . . . , π|Ω|(y), respectively. Let `(y) be the index such that π1(y)+. . . +π`(y)−1≤1−α(7) π1(y)+. . . +π`(y)>1−α. (8) Then, according to Rockafellar and Uryasev (2002), ηα(y) = R(y;ξ`(y)). Let Ω(y) = {ωr(y)∈Ω|1≤r≤`(y)}and let β(y) be the weight vector induced by 6 Ω(y), whose r-th component is given by βr(y)=          πr(y) 1−αif ωr(y)∈Ω(y)\ {ω`(y)} 1−Pω∈Ω(y)\{ω`(y)}πω 1−αif ωr(y) = ω`(y) 0 if ωr(y)∈Ω\Ω(y). (9) Since the non-zero components of β(y) determine the conditional probabilities for the scenarios in Ω(y), the α-CVaR of ycan be expressed as: ϑ(β(y),y) := f(y) + |Ω| X r=1 βr(y)R(y,ξr(y)).(10) As pointed out by F´abi´an (2008) the above β(y) vectors correspond to the dual optimal solution of the linear programming problem emerging from (3)–(6) when y∈Yis fixed. Such weights can also be looked at as the optimal solution of the linear relaxation of a knapsack problem. This is an aspect that has been explored in the context of k-sum optimization by Bertsimas and Sim (2003, 2004); Puerto et al. (2017); and Punen (1992). 3. Algorithmic framework for minimizing CVaR In this section we exploit the results of the above section to derive our algorithmic framework for the CVaR. Note that the definition of the CVaR associated with a given a priori solution y∈Yproduces the following non-linear formulation for CVaRα: CVaRαmin ϑ(β(y),y) := f(y) + |Ω| X r=1 βr(y)R(y,ξr(y)) (11) s. t. y∈Y. (12) One difficulty of this formulation, which is independent of the representation of the domain Y, is the non-linearity of the objective function. This can be partially overcome by using a fixed weight vector β∈R|Ω|. The resulting subproblem is: H(β) = min f(y) + X ω∈Ω βωR(y,ξω) (13) s. t. y∈Y. (14) Of course, there is no guarantee that when β∈R|Ω|is fixed, H(β) produces an optimal solution to CVaRα. However, as we will see, when β∈R|Ω|is suitably chosen, from an optimal solution to H(β) we can derive both a lower and an upper bound on the optimal 7 value of CVaRα. In this context, a key issue is how to select appropriate weight vectors β∈R|Ω|that provide useful bounds. For this we will choose suitable indices `∈ {1, . . . |Ω|} and scenario sets Ω`⊆Ω with ω`∈Ω`that jointly satisfy X ω∈Ω`\{ω`} πω≤1−α(15) X ω∈Ω` πω>1−α. (16) In the following such subsets Ω`are referred to as `-tails for α. Any `-tail for αnaturally induces the following vector βthat we refer to as `-weight vector: βω=          πω 1−αω∈Ω`\ {ω`} 1−Pω∈Ω`\{ω`}πω 1−αω=ω` 0ω∈Ω\Ω`. (17) The particular case Ω`= Ω(y) is referred to as scenario subset induced by y∈Y, and the `-weight vector β=β(y) is called scenario vector induced by Ω(y). In principle, any `-weight vector induced by a subset of scenarios defining an `-tail, or any other vector in the convex hull of such vectors, can be useful choices. In our framework we will indistinctively use both weight vectors induced by generic `-tails as well as scenario vectors. Remark 1. The expected conditional cost of an a priori solution yrelative to any subset of scenarios defining an `-tail, Ω`6= Ω(y), cannot exceed the α-CVaR of y. This includes the scenario subsets Ω(y0)induced by different a priori solutions y0∈Y,y06=y. This is due to the sorting used to define the scenario subset Ω(y), which guarantees that this subset is precisely the set of worse scenarios for the α-quantile for y, and to the fact that β(y)defines the conditional probabilities for this set of scenarios. In what follows we will denote by uthe set of the second-stage variables. We further denote by Q(y) the set of feasible solutions uassociated with a given y∈Y, and by Q={(y,u)|u∈Q(y),y∈Y}the feasible domain for joint vectors. We assume that the problem we are dealing with has relatively complete recourse, i.e., Q(y)6=∅, for all y∈Y. This is motivated by the application we present in Section 4. When β=β(y) is the scenario vector induced by Ω(y) for some y∈Y,H(β) can be interpreted as the problem of finding the best solution q= (y0,u)∈Qfor the scenario set induced by y. Since y0does not necessarily coincide with y, this perspective is complementary to that of the computation of the scenario vector β(y), which, is to find the set of scenarios that define the actual CVaRα for a given solution y. 8 Next we analyze some of the properties of subproblems H(β). First we note that if q= (y,u) is an optimal solution to H(β), then q∈Q, since otherwise it would not be optimal. Let us denote by v(β,q) the objective function value (13) for solution q= (y,u)∈Qfor H(β). Note that v(β(y),q) = ϑ(β(y),y). Remark 2. (i)The domain of H(β)is exactly the same as that of CVaRα, and only their objective functions are different. Thus, any feasible solution to H(β)produces a valid upper bound to CVaRα, provided that it is evaluated with respect to the function ϑ(β(y),y). Therefore, for any q= (y,u)∈Q,v(β(y),q)is an upper bound on the optimal value to CVaRα. (ii)Let q= (y,u)∈Q, with associated scenario vector β(y). Then, vβ0,q≤v(β(y),q) for any `-weight vector β06=β(y). This property is just a restatement of Remark 1, taking into account that if β06=β(y), then the `-tail that induces β0will also be different from Ω(y). Below we prove that when βis an `-weight vector the optimal value to H(β) gives a valid lower bound for CVaRα. Proposition 1. Let v∗be the optimal value to CVaRα. Let also ¯ q= (¯ y,¯ u)∈Qbe an optimal solution to H(β)for an `-weight vector β. Then, v(β,¯ q)≤v∗. Proof. Let q∗= (y∗,u∗)∈Qdenote an optimal solution to CVaRαwith associated scenario vector β∗=β(y∗). Consider also an optimal solution to H(β), ¯ q= (¯ y,¯ u)∈Q, with value v(β,¯ q). Then v(β,¯ q)≤v(β,q∗)≤v(β∗,q∗) = v∗, where the first inequality follows from the optimality of ¯ qfor H(β) and the second one from the definition of conditional value at risk of solution q∗and Remark 2.(ii). The next result allows to determine if a given solution ¯ q= (¯y,¯u)∈Qis optimal for CVaRα. Proposition 2. Let ¯ q= (¯ y,¯ u)∈Qbe an optimal solution to H(¯ β)for an `-weight vector ¯ β. If the scenario vector induced by Ω(¯ y)coincides with ¯ β, i.e. β(¯ y) = ¯ β, then ¯ qis an optimal solution to CVaRα. Proof. Let ¯ q= (¯ y,¯ u)∈Qbe an optimal solution to H(¯ β) for the `-weight vector ¯ β. By Proposition 1, v¯ β,¯ q≤v∗. On the other hand, by Remark 2.(i), v∗≤v(β(¯ y),¯ q). Finally, since β(¯ y) = ¯ βwe have, v(β(¯ y),¯ q) = v¯ β,¯ q≤v∗≤v(β(¯ y),¯ q), and the result follows. 9 “# Iter” is the (average) number of iterations required by Algorithm 1. This information is depicted only for ALG1-2,5% since this is the only approach for which this information makes sense. Finally, three columns present different gaps. “Gap (%) UB-LB” is the (average) percentage gap considering the best upper and lower bounds obtained in the corresponding approach. For each instance, it is computed as (UB −LB)/UB ×100. “Gap (%) UB” is the (average) percentage gap computed using the objective values of the best feasible solution obtained by the corresponding approach, say ¯v, and the best feasible solution found using the standard formulation (26)–(29), say ¯vGEN. Hence, for each instance, “Gap (%) UB” is computed as [¯v−¯vGEN]/¯vGEN ×100. We do not present such gap for approach GEN since it is trivially equal to zero. “Gap (%) LB” is the (average) percentage gap computed using the objective value of the best feasible solution obtained using the standard formulation (26)– (29)—¯vGEN—and the best lower bound obtained in the corresponding approach, say v. For each instance, it is obtained as [¯vGEN −v]/¯vGEN ×100. Again, we do not present this gap for approach GEN since in this case it coincides with Gap (%) UB-LB. We note that GEN is the “standard” existing approach. This motivates the use of the corresponding upper bound for computing two of the three gaps just described. Despite the difficulty of drawing overall comparisons when observing the extensive information provided by Tables 1–3, there are a few aspects that can be highlighted. First, the time limit of two hours was often reached when using either GEN or ALG1. Nevertheless, the final gap—“Gap (%) UB-LB”—is in general below 2%. When this is not the case, we often find the final gap provided by ALG1 outperforming that of GEN (e.g., |I|= 20, |J|= 100, |Ω|= 30). The superiority of ALG1 over GEN seems stronger when we consider the largest instances tested (|I|=50; |Ω|= 20,30). Second, observing Tables 2 and 3 we see that in several cases the average “Gap (%) UB” is negative for ALG1 and ALG1-2,5% indicating a better upper bound provided by these approaches when compared to GEN. Again, this is more evident for larger instances (Table 3). Third, considering either GEN-2,5% or ALG1-2,5%, in most of the cases, a final gap— “Gap (%) UB-LB”—less than or equal to 2.5% was obtained in a fairly small computing time. Finally, regarding the quality of the best lower bound (when compared to the upper bound provided by the standard approach GEN), we observe that, in general, ALG1 performs the best. In particular, one single execution of Algorithm 1 seems to be enough for obtaining a sharp lower bound. Although true, the above observations are not enough for drawing a clear conclusion as for the superiority of our new approach when compared to the standard one. Below we present results, which provide further insight into the overall performance of the compared alternatives. 16 Table 1: Results for the instances with 10 origins. CPU time Gap (%) CPU time Gap (%) |J| |Ω|Approach (seconds) # Iter UB-LB UB LB |J| |Ω|Aproach (seconds) # Iter UB-LB UB LB 20 8 GEN 2200 — 0.04 — — 100 8 GEN 5824 — 0.06 — — ALG1 3299 — 0.03 0.00 0.03 ALG1 5827 — 0.05 0.00 0.05 GEN-2,5% <1 — 1.61 1.19 0.44 GEN-2,5% <1 — 1.80 1.71 0.12 ALG1-2,5% <1 1.00 1.11 0.63 0.48 ALG1-2,5% 1 1.00 0.94 0.82 0.12 12 GEN 1850 — 0.01 — — 12 GEN 7202 — 0.20 — — ALG1 898 — 0.29 0.26 0.03 ALG1 7201 — 0.78 0.54 0.25 GEN-2,5% 1 — 1.40 0.88 0.53 GEN-2,5% 3 — 1.86 1.60 0.30 ALG1-2,5% 1 1.00 1.43 0.96 0.48 ALG1-2,5% 2 1.20 1.90 1.59 0.34 20 GEN 6222 — 0.29 — — 20 GEN 7232 — 0.39 — — ALG1 5999 — 0.74 0.45 0.30 ALG1 7203 — 1.13 0.69 0.45 GEN-2,5% 6 — 1.68 0.97 0.73 GEN-2,5% 314 — 0.70 0.29 0.41 ALG1-2,5% 11 2.60 1.57 0.77 0.81 ALG1-2,5% 229 1.80 2.14 1.70 0.48 30 GEN 7203 — 0.43 — — 30 GEN 7203 — 0.50 — — ALG1 7203 — 0.68 0.25 0.43 ALG1 7203 — 0.63 0.11 0.52 GEN-2,5% 14 — 1.86 1.16 0.72 GEN-2,5% 882 — 0.74 0.24 0.51 ALG1-2,5% 12 1.00 1.58 0.86 0.73 ALG1-2,5% 1143 1.00 1.25 0.73 0.54 50 8 GEN 7236 — 0.11 — — ALG1 7097 — 0.11 0.00 0.11 GEN-2,5% <1 — 1.88 1.52 0.39 ALG1-2,5% <1 1.00 1.88 1.52 0.39 12 GEN 7203 — 0.42 — — ALG1 7204 — 0.69 0.24 0.45 GEN-2,5% 4 — 1.33 0.71 0.64 ALG1-2,5% 5 1.60 1.35 0.72 0.65 20 GEN 7219 — 0.60 — — ALG1 7203 — 1.37 0.78 0.61 GEN-2,5% 90 — 1.45 0.76 0.71 ALG1-2,5% 79 1.00 1.79 1.09 0.72 30 GEN 7201 — 0.73 — — ALG1 7202 — 0.95 0.19 0.76 GEN-2,5% 262 — 1.34 0.53 0.82 ALG1-2,5% 263 1.00 1.35 0.52 0.84 17 Table 2: Results for the instances with 20 origins. CPU time Gap (%) CPU time Gap (%) |J| |Ω|Approach (seconds) # Iter UB-LB UB LB |J| |Ω|Aproach (seconds) # Iter UB-LB UB LB 20 8 GEN 5329 — 0.05 — — 100 8 GEN 7230 — 0.19 — — ALG1 5025 — 0.04 -0.01 0.04 ALG1 7202 — 0.17 -0.02 0.19 GEN-2,5% <1 — 2.26 1.93 0.37 GEN-2,5% 1 — 1.82 1.52 0.32 ALG1-2,5% <1 1.00 2.21 1.75 0.49 ALG1-2,5% 3 1.00 1.36 1.06 0.32 12 GEN 5995 — 0.12 — — 12 GEN 7271 — 0.46 — — ALG1 6135 — 0.98 0.85 0.16 ALG1 7203 — 1.59 1.10 0.51 GEN-2,5% 3 — 1.01 0.61 0.40 GEN-2,5% 325 — 1.41 0.90 0.53 ALG1-2,5% 3 1.00 1.42 0.98 0.45 ALG1-2,5% 23 1.00 2.27 1.63 0.68 20 GEN 7234 — 0.27 — — 20 GEN 7219 — 1.49 — — ALG1 7203 — 1.02 0.68 0.35 ALG1 7206 — 1.17 -0.34 1.51 GEN-2,5% 9 — 1.56 0.94 0.63 GEN-2,5% 5150 — 1.86 0.37 1.50 ALG1-2,5% 14 1.00 1.36 0.85 0.53 ALG1-2,5% 4267 1.00 2.09 0.60 1.51 30 GEN 7203 — 0.42 — — 30 GEN 7200 — 5.35 — — ALG1 7203 — 0.57 0.12 0.45 ALG1 7212 — 3.20 -2.24 5.37 GEN-2,5% 34 — 1.41 0.79 0.63 GEN-2,5% 7200 — 5.35 0.00 5.35 ALG1-2,5% 55 1.20 1.14 0.60 0.55 ALG1-2,5% 7200 1.00 3.23 -2.21 5.37 50 8 GEN 7229 — 0.24 — — ALG1 7203 — 0.24 0.00 0.24 GEN-2,5% <1 — 2.28 1.71 0.60 ALG1-2,5% <1 1.00 2.15 1.54 0.64 12 GEN 7202 — 0.48 — — ALG1 7202 — 1.23 0.76 0.50 GEN-2,5% 50 — 1.10 0.51 0.60 ALG1-2,5% 63 3.20 2.30 1.72 0.63 20 GEN 7200 — 0.81 — — ALG1 7202 — 1.37 0.55 0.83 GEN-2,5% 1308 — 1.12 0.28 0.84 ALG1-2,5% 1239 1.60 1.79 0.94 0.87 30 GEN 7200 — 0.91 — — ALG1 7203 — 1.21 0.28 0.93 GEN-2,5% 2273 — 1.26 0.34 0.92 ALG1-2,5% 1285 1.20 1.70 0.75 0.96 18 Table 3: Results for the instances with 50 origins. CPU time Gap (%) CPU time Gap (%) |J| |Ω|Approach (seconds) # Iter UB-LB UB LB |J| |Ω|Aproach (seconds) # Iter UB-LB UB LB 50 8 GEN 7201 — 0.14 — — 100 8 GEN 7200 — 0.20 — — ALG1 7202 — 0.13 -0.01 0.14 ALG1 7207 — 0.20 0.01 0.20 GEN-2,5% 2 — 1.24 0.94 0.31 GEN-2,5% 5 — 1.15 0.71 0.44 ALG1-2,5% 3 1.00 1.45 1.15 0.31 ALG1-2,5% 11 1.00 1.34 0.90 0.44 12 GEN 7200 — 0.29 — — — 12 GEN 7200 — 1.58 — — ALG1 7204 — 1.49 1.20 0.31 ALG1 7214 — 1.94 0.32 1.62 GEN-2,5% 6 — 2.12 1.60 0.55 GEN-2,5% 16 — 2.04 0.13 1.92 ALG1-2,5% 46 1.40 2.04 1.51 0.56 ALG1-2,5% 667 4.80 2.07 0.30 1.77 20 GEN 7200 — 2.62 — — 20 GEN 7200 — 2.32 — — ALG1 7208 — 2.28 -0.36 2.63 ALG1 7239 — 1.79 -0.54 2.32 GEN-2,5% 5764 — 2.67 0.01 2.67 GEN-2,5% 744 — 2.43 0.00 2.43 ALG1-2,5% 445 1.20 2.34 -0.33 2.66 ALG1-2,5% 1168 1.00 2.19 -0.21 2.39 30 GEN 7200 — 3.22 — — 30 GEN 7200 — 2.96 — — ALG1 7219 — 2.61 -0.64 3.23 ALG1 7316 — 2.07 -0.92 2.97 GEN-2,5% 5913 — 3.22 0.00 3.22 GEN-2,5% 7200 — 2.96 0.00 2.96 ALG1-2,5% 18094 2.60 2.58 -0.66 3.23 ALG1-2,5% 4040 1.00 2.09 -0.92 2.99 19 50 60 70 80 90 100 0,20% 0,40% 0,75% 1,10% 1,60% 2,10% 2,60% 3,10% 3,60% % of instances GEN ALG1 (a) Whole set of instances. 50 60 70 80 90 100 0,20% 0,40% 0,75% 1,10% 1,60% 2,10% 2,60% 3,10% 3,60% % of instances GEN ALG1 (b) All instances with |Ω|= 30. 0 10 20 30 40 50 60 70 80 90 100 0,20% 0,40% 0,75% 1,10% 1,60% 2,10% 2,60% 3,10% 3,60% % of instances GEN ALG1 (c) Larger instances (|I|= 50, |J|= 100, |Ω|= 30). Figure 1: Performance profile with respect to the Gap. 5.2. Performance profiles In this section we analyze the performance profiles of the different approaches tested. We recall that a performance profile graph plots the percentage of instances solved for each value of the corresponding measure. We start with a gap performance profile for GEN and ALG1 with a time limit of two hours. For this, we consider the gaps between the best upper bound provided by each approach and the best upper bound among the two approaches. More specifically, let ¯vGEN be the objective function value for the best feasible solution found using GEN and ¯vALG1 the corresponding value for ALG1. Define v∗= min{¯vGEN,¯vALG1}. For GEN and ALG1 we compute the gap according to (¯vGEN −v∗)/v∗×100 and (¯vALG1 −v∗)/v∗×100, respectively. The gap profiles are depicted in Figure 1. As we can see in Figure 1(a), when we consider the whole set of instances the performance profile with respect to the gap of the standard formulation (26)–(29) is slightly better than the performance profile of our alternative solution approach. However for all instances with |Ω|= 30 and for the larger instances —Figures 1(b) and 1(c)— the performance profile of our new approach ALG1 is superior to that of the standard formulation. In this comparison, it is worth-noting that we are considering only one iteration of Algorithm 1. A similar conclusion was already sketched in the previous section from the analysis of the extensive results presented in Tables 1–3. Figure 2 depicts performance profiles with respect to computing times for GEN and ALG1. 20 0 10 20 30 40 50 60 70 80 90 100 1' 5' 10' 30' 60' 90' 120' >=120' % of instances GEN ALG1 Figure 2: Performance profile with respect to computing times (all instances). All instances were considered when drawing this profiles. As we can see, both approaches exhibit a similar behavior. In this figure, we also get a clear empirical evidence of the difficulty for solving CVaRαTP. As can be seen, only 10% of the instances were solved to optimality within the time limit of two hours with any of the two tested methods. Due to the difficulty in solving CVaRαTP to proven optimality, we focused on the behavior of the two variants of the proposed approaches, when a 2,5% MIP gap is set as a termination criterion when solving the MIP formulations. The corresponding performance profiles for the computing times are depicted in Figure 3. In this case, the superiority of the new modeling framework (ALG1-2,5%) as compared to the standard one (GEN-2,5%) becomes visible. Nonetheless, it is stronger for the largest instances tested as shown in Figure 3(c). The superiority of ALG1-2,5% becomes more noticeable if we consider performance profiles with respect to the Gap(%) UB-LB, as depicted in Figure 4. In the particular cases of Figures 4(b) and 4(c) we see that ALG1-2,5% provides in general smaller gaps. This fact, combined with a better performance in terms of computing times, makes clear the superiority of ALG1-2,5% when compared to the standard approach GEN-2,5%. Finally, Figure 5 gives the performance profile of ALG1-2,5% with respect to the number of iterations. As we can observe, in 90% of the cases, Algorithm 1 stopped after performing just one iteration since that was enough to obtain a solution with a MIP gap smaller than or equal to 2,5%. 5.3. Further testing—large-scale instances In order to give more evidence of the superiority of our proposed solution framework relative to the standard one, we have compared both of them on one additional set of 15 largescale instances, generated as proposed in Hinojosa et al. (2014). In all cases we considered |Ω|= 50 scenarios. For |I|= 50 we took |J|= 100 and 200; for |I|= 100 we considered |J|= 200. Five instances have been generate for each combination of |I|,|J|and |Ω|. The obtained results are summarized in Table 4. We do not report computing times, because the two hours time limit was reached in all cases. Moreover, neither the final gap at termination is reported for approaches GEN-2,5% and ALG1-2,5% with the largest instances (|I|= 100, 21 50 60 70 80 90 100 1' 5' 10' 30' 60' 90' 120' >=120' % of instances GEN-2.5% ALG1-2.5% (a) Whole set of instances. 20 30 40 50 60 70 80 90 100 1' 5' 10' 30' 60' 90' 120' >=120' % of instances GEN-2.5% ALG1-2.5% (b) All instances with |Ω|= 30. 0 10 20 30 40 50 60 70 80 90 100 1' 5' 10' 30' 60' 90' 120' >=120' % of instances GEN-2.5% ALG1-2.5% (c) Larger instances (|I|= 50, |J|= 100, |Ω|= 30). Figure 3: Performance profile with respect to computing time (MIP gap 2,5%). 55 60 65 70 75 80 85 90 95 100 2,00% 2,25% 2,50% 3,00% 3,50% 5,00% >=5% % of instances GEN-2.5% ALG1-2.5% (a) Whole set of instances. 55 60 65 70 75 80 85 90 95 100 2,00% 2,25% 2,50% 3,00% 3,50% 5,00% >=5% % of instances GEN-2.5% ALG1-2.5% (b) All instances with |Ω|= 30. 0 10 20 30 40 50 60 70 80 90 100 2,00% 2,25% 2,50% 3,00% 3,50% 5,00% >=5% % of instances GEN-2.5% ALG1-2.5% (c) Larger instances (|I|= 50, |J|= 100, |Ω|= 30). Figure 4: Performance profile with respect to Gap (MIP gap 2,5%). 22 50 60 70 80 90 100 1|Ω|/4 |Ω|/2 |Ω| % of instances ALG1-2.5% Figure 5: Performance profile with respect to the number of iterations. |J|= 200), as it is always greater than 2,5%. For these new instances, our new approach clearly outperforms the standard one in terms of the average gap attained. In some cases, particularly for all cases with |J|= 200 the difference is very significant. Table 4: Additional results. |I| |J| |Ω|Tested Approach Gap (%) UB-LB # Iter 50 100 50 GEN 5,18 — ALG1 3,57 — GEN-2,5% 5,18 — ALG1-2,5% 3,56 1,2 200 50 GEN 80,88 — ALG1 6,10 — GEN-2,5% 80,88 — ALG1-2,5% 5,36 1,4 100 200 50 GEN 22,23 — ALG1 3,46 — 6. Conclusions In this paper, we proposed a new algorithm for obtaining lower and upper bounds on the optimal value of a stochastic mixed-integer programming problem considering the minimization of the conditional value at risk. The methodological developments were applied to a fixed-charge stochastic transportation problem. The results show that the new approach is competitive when compared to the standard one. This is particularly true for the largest instances tested. One drawback of our approach is the need to solve one MIP problem in each iteration. Therefore, the competitiveness of this approach could possibly be higher if we considered a specially-tailored procedure for solving those MIPs. This depends of course on the specific application being considered. Nevertheless, this is a research direction to explore when making use of our new algorithmic scheme. 23 Related with the previous aspect, it is worth investigating the application of our new algorithm to other problems. This is relevant not only because the new approach is a promising alternative to the existing one but also because that would allow evaluating the robustness of our approach across different experimental settings. Acknowledgement The research of the first three authors was partially supported by the Spanish Ministry of Economy and Competitiveness through MINECO/FEDER grants MTM2015-63779-R and MTM2016-74983-C02-01. The research of the fourth author was partially supported by the Portuguese Science Foundation (FCT—Funda¸c˜ao para a Ciˆencia e Tecnologia) under the project UID/MAT/ 04561/2013 (CMAF-CIO/FCUL). References References Ahmed, S. (2006). Convexity and decomposition of mean-risk stochastic programs. Mathematical Programming, 106:433–446. Ahuja, R., Magnanti, T., and Orlin, J. (1993). Network Flows: Theory, Algorithms, and Applications. Prentice Hall. Alonso-Ayuso, A., Escudero, L., Guignard, M., and Weintraub, A. (2016). Risk management for forestry planning under uncertainty in demand and prices. Optimization Online. http://www.optimization-online.org/DB HTML/2016/12/5756.html. Balinski, M. (1961). Fixed-cost transportation problems. Naval Research Logistics Quaterly, 8:41–54. Bertsimas, D. and Sim, M. (2003). Robust discrete optimization and network flows. Mathematical Programming Series B, 98:49–71. Bertsimas, D. and Sim, M. (2004). The price of robustness. Operations Research, 52:35–53. Birge, J. and Louveaux, F. (2011). Introduction to stochastic Programming. Springer, Berlin- Heidelberg. Second edition. Carøe, C. and Schultz, R. (1999). Dual decomposition in stochastic integer programming. Operations Research Letters, 24:37–45. Chelst, K. and Canbolat, Y. (2011). Value-added decision making for managers. CRC Press, Boca Raton-London-New York. 24 Escudero, L., Gar´ın, M., and Unzueta, A. (2017). Scenario cluster Lagrangean decomposition for risk averse in multistage stochastic optimization. Computers & Operations Research, 85:1339–1351. F´abi´an, C. (2008). Handling CVaR objectives and constraints in two-stage stochastic models. European Journal of Operational Research, 191:888–911. Grieco, S., Semeraro, U., and Tolio, T. (2001). A review of different approaches to the FMS loading problem. International Journal of Flexible Manufacturing Systems, 13(4):361–384. Hemmati, R., Saboori, H., and Saboori, S. (2016). Stochastic risk-averse coordinated scheduling of grid integrated energy storage units in transmission constrained wind-thermal systems within a conditional value-at-risk framework. Energy, 113:762–775. Hinojosa, Y., Puerto, J., and Saldanha-da-Gama, F. (2014). A two-stage stochastic transportation problem with fixed handling costs and a priori selection of the distribution channels. TOP, 22:1123–1147. Homem-de-Mello, T. and Pagnoncelli, B. (2016). Risk aversion in multistage stochastic programming: A modeling and algorithmic perspective. European Journal of Operational Research, 249:188–199. K¨unzi-Bay, A. and Mayer, J. (2006). Computational aspects of minimizing conditional valueat-risk. Computational Management Science, 3:3–27. List, G. F., Wood, B., Nozick, L. K., Turnquist, M. A., Jones, D. A., Kjeldgaard, E. A., and Lawton, C. R. (2003). Robust optimization for fleet planning under uncertainty. Transportation Research Part E: Logistics and Transportation Review, 39:209–227. Liu, C., Fan, Y., and Ord´o˜nez, F. (2009). A two-stage stochastic programming model for transportation network protection. Computers & Operations Research, 36:1582–1590. Meng, Q., Wang, T., and Wang, S. (2012). Short-term liner ship fleet planning with container transshipment and uncertain container shipment demand. European Journal of Operational Research, 223:96–105. Miller, N. and Ruszczy´nski, A. (2011). Risk-averse two-stage stochastic linear programming: modeling and decomposition. Operations Research, 59:125–132. Nickel, S., da Gama, F. S., and Ziegler, H.-P. (2012). A multi-stage stochastic supply network design problem with financial decisions and risk management. Omega, 40:511–524. Noyan, N. (2012). Risk-averse two-stage stochastic programming with an application to disaster management. Computers & Operations Research, 39:541–559. 25