A two-stage robust optimization approach for the mobile facility fleet sizing and routing problem under uncertainty

Author’s Accepted Manuscript A two-stage robust optimization approach for the mobile facility fleet sizing and routing problem under uncertainty Chao Lei, Wei-Hua Lin, Lixin Miao www.elsevier.com/locate/caor

PII: DOI: Reference:

S0305-0548(15)00217-8 http://dx.doi.org/10.1016/j.cor.2015.09.007 CAOR3858

To appear in: Computers and Operation Research Received date: 24 October 2014 Revised date: 7 July 2015 Accepted date: 17 September 2015 Cite this article as: Chao Lei, Wei-Hua Lin and Lixin Miao, A two-stage robust optimization approach for the mobile facility fleet sizing and routing problem under uncertainty, Computers and Operation Research, http://dx.doi.org/10.1016/j.cor.2015.09.007 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A two-stage robust optimization approach for the mobile facility fleet sizing and routing problem under uncertainty Chao Leia,c , Wei-Hua Linb , Lixin Miaoc b

a Department of Industrial Engineering, Tsinghua University, Beijing 100084, China Department of Systems and Industrial Engineering, The University of Arizona, Tucson, AZ 85721, USA c Research Center for Modern Logistics, Graduate School at Shenzhen, Tsinghua University, Shenzhen 518055, China

Abstract We propose a two-stage robust optimization model for the mobile facility fleet sizing and routing problem with demand uncertainty. A two-level cutting plane based method is developed, which includes an algorithm to generate problem-specific lower bound inequalities in the outer level, and a hybrid algorithm in the inner level that combines heuristic and exact methods to solve the recourse problem. Numerical tests show that the design and operation from the proposed method outperforms other solution approaches. The efficiency of the proposed solution algorithm in identifying the optimal solution is quantified and the robustness of the proposed model is demonstrated for varying degrees of uncertainty in demand. Keywords: Robust optimization, Fleet management, Mobile facility 1. Introduction The Traveling Salesman Problem (TSP) and the Vehicle Routing Problem (VRP), as classical operations research problems, have been studied extensively to handle the planning and scheduling of the vehicle movement within Email addresses: [email protected] (Chao Lei), [email protected] (Wei-Hua Lin), [email protected] (Lixin Miao)

Preprint submitted to Computers & Operations Research

September 23, 2015

a network. Vehicles in those problems usually provide pickup or delivery services between facilities (e.g. depots or distribution centers) and customers (e.g. retailers or end users). However, in reality, vehicles themselves can sometimes act like a facility to provide real-time services to customers when they are stationary at one location. More importantly, they have the advantage of moving from one place to another like ordinary vehicles to cover a broader region, a feature especially important for regions that exhibits a large variability in demand over time. For instance, light trucks equipped with cellular base stations provide cellular services to areas where the established cellular network is temporarily insufficient (Halper and Raghavan, 2011). Also, vans provide fast food services to customers in different regions in different time periods of the day (Lei et al., 2014). Note that when these vehicles are in service, they behave like traditional facilities, and the duration that the vehicle stays at a customer point is directly related to the amount of service provided. In this paper, we use the term Mobile Facility (MF) to denote this “facility-like-vehicle”. The MFs considered here share the same features as those shown in Lei et al. (2014): 1) the demand served at a location is proportional to the duration of time an MF remains stationary at that location, 2) the MFs cannot provide services when they are in motion, and 3) the service begins when an MF arrives at a location and ends when it departs. C1

L1

L3

L2 C2

C3

Figure 1: Schematic representation of customers and candidate locations

We illustrate these features through a simple example as shown in Figure 1. There are three customer points (C1, C2 and C3) and three candidate locations (L1, L2 and L3) where the MF can stay. We consider a two-hour planning horizon (from 8:00AM to 10:00AM) with eight time periods (i.e., 2

Table 1: Customer demand in different time period

Customer C1 C2 C3

8:00 ∼8:15 10 0 0

8:15 ∼8:30 10 0 0

8:30 ∼8:45 10 0 0

Time 8:45 9:00 ∼9:00 ∼9:15 0 0 0 0 15 15

9:15 ∼9:30 0 0 15

9:30 ∼9:45 0 20 0

9:45 ∼10:00 0 20 0

15-minutes periods). Each customer is covered by a single location as shown in Figure 1. The travel time between any two locations is assumed to be 15 minutes. The amount of demand in each time period is given in Table 1. 8:00

8:15

8:30

8:45

9:00

9:15

9:30

9:45

10:00

L1

L1

L1

L1

L1

L1

L1

L1

L1

L2

L2

L2

L2

L2

L2

L2

L2

L2

L3

L3

L3

L3

L3

L3

L3

L3

L3

Figure 2: Route for a single MF

In Figure 2, we show the service schedule of a single MF. The horizontal arrows indicate that the MF is stationary at a certain location and is available to provide service. As shown in the figure, the MF stays at L1 for two time periods from 8:00, covering 20 units of demand from C1. It then departs from L1 at 8:30 and arrives at L3 15 minutes later. When the MF is on the road, it is unable to provide services, and the amount served from 8:30-8:45 is zero. The MF then proceeds to L3 and finally ends its services at L2 as shown in the figure. We can see that the MF will travel to another location only if the gain from covering more demands at a different location outweighs the loss resulting from the movement. Fleet sizing is one of the most important decisions as it is a major fixed investment for starting any business. Moreover, fleet sizing is a complex longterm decision which needs to be considered in conjunction with the operation level. 3

First, it is essential to come up with certain strategies to handle the uncertainty in demand in operation when making long-term investment decisions. In this paper, we focus on the robust optimization (RO) approach of modeling the proposed problem under demand uncertainty. The approach is appealing in three aspects: 1) it generates a solution immune to uncertain demand within a deterministic uncertainty set. The RO approach is quite desirable since the establishment of a fleet typically involves a significant investment in the starting stage of the business, including the cost of purchasing an MF, expenses for devices to be equipped on an MF for providing relevant service, the cost of the crew who operates the equipment, the routine maintenance cost, etc.; 2) the RO approach requires only moderate information about the uncertain demand data, rather than a detailed description of the probability distribution function or a large collection of sample scenarios; 3) the RO approach generates solutions with different level of conservativeness by manipulating the “budget of uncertainty” to guard against overly conservative results. The managers thus can be provided with more alternatives so that they can decide based on their specific needs. In this paper, the uncertain parameters are assumed to be within a polyhedral uncertainty set in line with those assumed in Bertsimas and Sim (2003, 2004). Second, the fleet size is strongly dependent on the effectiveness in fleet (such as the route and schedule planning). Moreover, owing to the “facilitylike-vehicle” feature of the MF, the allocation of demands to stationary MFs is significant for the entire system performance. Hence, it is more reasonable to integrate the strategic level fleet sizing problem, the tactical level routing and scheduling problem, and the operational level allocation problem together into the decision making process. We name this type of problem as the mobile facility fleet sizing and routing problem (MFFSRP). This paper considers the MFFSRP by applying a two-stage robust optimization modeling approach, in which the customer demand is uncertain with no information about the underlying probability distribution function. The idea of the two-stage RO approach is that the first-stage decisions (i.e., fleet sizing and routing plan) are made before the realization of uncertainty, and then the second-stage decisions (i.e., allocation of demands to the MFs) are made as recourse actions after the uncertainty has been known. The objective is to optimize the performance of the system over the entire planning horizon. The computational aspect of the proposed problem must be carefully addressed. In addition to the integration of the fleet sizing and the routing 4

problem, some other factors, such as the two-stage feature, uncertainty in demand, and multiple time periods of planning, all increase the computational complexity significantly. In this paper, we develop some new techniques to improve the performance of the conventional cutting plane method. To the best of our knowledge, this study is the first one applying the RO approach in the context of the mobile facility planning field. Specifically, the contributions of this paper are as follows: 1. A two-stage RO approach is applied to the formulation of the MFFSRP for the first time, in which the uncertainty in demand is modeled as a polyhedral set to allow flexible control of the level of conservativeness. Numerical results reveal that the proposed RO approach has some distinctive advantages in optimizing the performance of the system compared to the overly conservative solution obtained from the box uncertainty set model. 2. An enhanced version of the cutting plane method is proposed. The enhancement includes: 1) substituting the single optimality cut by multiple cuts, 2) constructing several problem specific lower bound inequalities to tighten the lower bound of the master problem, and 3) developing a hybrid algorithm solver for the recourse problem. 3. An experiment is designed to examine the effectiveness of the proposed approach. The experiment goes beyond the commonly used measures for testing solution quality and computation time. Instead, we consider a wide class of performance measures to test the model performance, including economic effectiveness, reliability, robustness against variations in demand, solution efficiency, etc. The remainder of the paper is organized as follows. In Section 2, we review some of the relevant literature with respect to the fleet sizing problem, the mobile facility planning problem, and the robust optimization approach. Section 3 first gives the mathematical formulation for the nominal model of the MFFSRP and then the two-stage robust formulation with the polyhedral uncertainty set. A two-level cutting plane type algorithm is presented in Section 4. Extensive numerical experiments are carried out in Section 5, and insights from the results are presented. This is followed by Section 6 with conclusions summarizing the main features of the proposed approach.

5

2. Literature review In this section, we review the relevant literatures about the MFFSRP and the RO approach, especially the study of the mobile facility itself, the two-stage RO modeling, and the corresponding solution methodology. It should be noted that the concept of mobile facilities is not new. Bespamyatnikh et al. (2000) and Durocher and Kirkpatrick (2006) studied a problem related to mobile facility location, where both mobile facility and the customer could continuously move in the space. They aimed at finding certain trajectory for each of mobile facilities so that the distances from the set of continuously moving customers to the closest mobile facility are minimized in the space. The mobile facility in their study provides service while moving continuously, whereas in our study MFs have to be stationary to serve the demand. In addition, the customers in our study do not move around over time. Gendreau et al. (2001) studied an ambulance relocation problem whose objective was to maximize the amount of demand covered by at least two ambulances. Although the ambulance can be considered as a type of mobile facility, the underlying problem was fundamentally different from ours. The relocation of ambulances was driven by the need to maintain coverage standard when certain ambulances were dispatched to respond to a specific service request. The relocation of mobile facility in our study is driven by time dependent demand pattern of different customers. The most relevant research to our problem are studied by Halper and Raghavan (2011) and Lei et al. (2014). Halper and Raghavan presented a continuous time formulation to model the maximum covering mobile facility routing problem for both single and multiple facilities under deterministic conditions. Lei et al., proposed a two-stage stochastic programming model for the mobile facility routing and scheduling problem, in which the first stage decision determines the temporal and spatial plan of the MFs and the second stage handles the allocation of customer demands. None of them explicitly consider the fleet sizing of the MFs or incorporate the robust optimization approach into the MF planning problem under uncertainty. The MFFSRP is similar to a well-studied branch of the facility location problem: the Dynamic Facility Location Problem (DFLP). One can refer to Boloori Arabani and Farahani (2012) for a detailed review. While the MFFSRP and the DFLP both involve making time-dependent positioning decisions of certain facilities, the MFFSRP explicitly accounts for the relo6

cation time (i.e., the traveling time) of MFs whereas the DFLP does not. Another highly relevant problem is the team orienteering problem (TOP) or the selective vehicle routing problem (see the review by Vansteenwegen et al. (2011)), where a team of players collects rewards among different places to try to maximize the team-collected rewards within a certain time limit. But it does not consider the facility-like feature of the MFFSRP, and the amount of rewards collected at a certain place is not related to the length of stationary time. Fleet size planning models can be treated as a subset of capacity planning models, which involves determining the capacity needed to meet future demands. Literature review on capacity planning models can be found in Van Mieghem (2003). Fleet sizing decisions are often made in conjunction with other decisions, such as vehicle allocation (Beaujon and Turnquist, 1991; Sayarshad and Ghoseiri, 2009), vehicle routing (Hoff et al., 2010), or empty vehicle repositioning (Li and Tao, 2010; You and Hsieh, 2014). For fleet sizing models considering uncertainty, List et al. (2003) proposed a two-stage stochastic programming model for fleet sizing under uncertainty in future demands and operating conditions. Song and Earl (2008) considered a twostage stochastic model for fleet sizing and empty vehicle repositioning problem with uncertainties in arrival times for loaded vehicle and repositioning times for empty vehicles. In this paper, we consider the use of the two-stage robust optimization (RO) approach to model fleet sizing and routing of the mobile facility problem. The RO approach has also been broadly applied as a modeling framework for optimization under uncertainty in operations research, such as supply chain management (Bertsimas and Thiele, 2004, 2006; Carlsson et al., 2014), wireless networks planning (Claßen et al., 2013), vehicle routing (Agra et al., 2013), etc. The readers can refer to the review paper by Bertsimas et al. (2011) and Gabrel et al. (2014b) for further information of the RO theories and applications. Here we focus on the two-stage RO approach for its high relevance to our problem. Atamt¨ urk and Zhang (2007) pointed out that, by allowing a recourse action in the second stage, one is able to obtain less conservative solutions compared to the single-stage RO approach. However, the advantage usually comes with a price of greater computational challenge. Ben-Tal et al. (2004) proved that even a simple two-stage RO problem could be NP-hard. Two major solution strategies were proposed in the literature: the exact approach and the approximation approach. For the type of the approximation method, certain specific relationship between the second-stage 7

decisions and the uncertain parameters is typically assumed to simplify the problem. Ben-Tal et al. (2004) first introduced the affinely-adjustable robust counterpart (AARC) approach, which restricts the second-stage decision variables to be affine functions of the uncertain parameters such that the original problem can be transformed into a deterministic linear problem. Later, they applied this method to an emergency logistics planning problem (Ben-Tal et al., 2011). The relevant literatures in this category can also be found in Mudchanatongsuk et al. (2007), Chen and Zhang (2009), Erera et al. (2009) and Bertsimas et al. (2010). The focus of our literature review will primarily be on the exact approach since it is also the approach we have taken in this paper. Though exact approaches have many different variations, they are all based on the Benders decomposition method (Benders, 1962), which gradually converges to the optimal first-stage solution by using the dual information of the secondstage subproblem. Thiele et al. (2009) proposed a cutting plane algorithm, in which the second-stage recourse problem is solved in an exactly same way. Gabrel et al. (2014a) applied the Thiele et al.’s method to solve a two-stage RO model for the location transportation problem. Bertsimas et al. (2013) studied a two-stage RO model for the security constrained unit commitment problem, in which an outer approximation technique is applied to handle the bilinear issue in the recourse problem. Jiang et al. (2014) studied a similar problem as Bertsimas et al. (2013) by developing an alternating heuristic method to solve the recourse problem. 3. Mobile facility fleet sizing and routing problem Let G (V, E) be a connected graph with node set V and edge set E. The subset P ⊆ V is the set of customer nodes and the subset L ⊆ V is the set of nodes where MFs can be located. Let dij be the distance between any pair of nodes i and j. We assume that dij is symmetric, since this can be easily relaxed. We assume that the length of the planning horizon is T . For simplicity in modeling, we divide the planning horizon into |T | identical time periods, whose set is denoted as T = {τ |τ = 1, 2, . . . , |T |}. For any time period τ ∈ T , it represents the time interval between tτ −1 and tτ . For example, if τ = |T |, it begins at t|T |−1 and ends at t|T | , where t|T | = T . We denote the set of time points tτ as N = {t0 , t1 , . . . , t|T | }, in which t0 < t1 < . . . < t|T | .

8

W

W

W

W_7_

1

1

1

1

2

2

2

2

O

S j

j

j

j

_/_

_/_

_/_

_/_

Figure 3: Time-space network

Since both the spatial and temporal aspects have to be considered for the MFs, we propose to model the problem on a time-space network as shown in Figure 3. For such a network, each node represents a specific location j ∈ L at a particular time n ∈ N . Thus we denote the node set of location-time combinations in the time-space network as S, whose elements are (j, n) pairs. Let A be the set of arcs in the time-space network, which is indexed by a 4-tuple, (i, m, j, n), where (i, m) ∈ S, (j, n) ∈ S, and m < n. We use Yimjn as a decision variable for the flow over the arc (i, m, j, n). In the location-time network, the flow is the number of MFs transiting from node to node, which is a nonnegative integer variable. If i = j, Yimjn is the number of MFs staying at location i within the time interval [m, n); if i ̸= j, Yimjn is the number of MFs departing from location i at time m and arriving location j at time n. Since the fleet size needs to be decided at the very beginning of the planning horizon, we add a virtual source node and a virtual sink node into the location-time network to capture the total number of MFs entering the network. In Figure 3, there are only the outbound arcs for source node O and inbound arcs for sink node S. The outbound arcs of the source node, denoted as (O, 0, j, n), can only enter node (j, n) for n = t0 ; and in turn, the inbound arcs of sink node S, denoted as (j, n, S, T + 1), can only start from node (j, n) for n = T . ∑ The total number ∑ of MFs is defined as decision Y = variable Q which equals to j∈L YjnSN . j∈L O0jn n=t0

n=T

For each MF, a limit of capacity C is imposed on the amount of demand that can be covered in a single time period. Considering the expense of 9

purchasing or renting an MF, staffing cost, and equipment investment, we denote the fixed cost of establishing a MF as a nonnegative cost factor f . The travel time Hij from location i to j is predefined and measured as an integer multiplier of a single time period. Note that an MF in motion is unable to provide service, which is likely to result in a loss of benefit. Therefore, the explicitly consideration of the travel time would inherently discourage the uneconomic movement. To avoid repeatedly accounting for the cost of traveling, the travel cost factor is not considered in the model. Let Wpτ denote the demand at customer p in period τ . The demand of customer p can be covered by MFs located in j if only the distance dpj ≤ D, where D is the maximum coverage distance. We apply the 0-1 parameter apj to indicate the coverage relations between the customer nodes and candidate locations, i.e. either apj = 1 if dpj ≤ D, or apj = 0 otherwise. We define Xpjτ as the assignment variable representing the amount of demands of customer p, which is covered by MFs at location j in time period τ . Due to the coverage restriction and the limited capacity, it is possible that the demands cannot be fully covered by the entire MF fleet. We define a decision variable Spτ to represent the amount of unmet demands of customer p in time period τ , and we assume that a penalty cost β is incurred for each unit of unmet demand, where β is a nonnegative cost factor. The practical meaning of this penalty cost can either be considered as the opportunity cost for the unmet demands or regarded as an expense for outsourcing excess demands to other firms. The main purpose of the MFFSRP is to serve as many demands as possible in the given planning horizon at the lowest possible cost. We emphasize the tradeoff between the investment for establishing a suitable fleet of MFs (i.e., the fixed cost) and the loss of benefits for failing to satisfy customers (i.e., the penalty cost for unmet demands). The other cost factors are assumed to be relatively less important that will not be considered in the tradeoff. To get a better understanding of the tradeoff condition, we continue to use the small problem in Section 1 as the example. Assume the fixed cost factor f = 100 and the penalty cost factor β = 2. If there is only one MF in the fleet, the optimal MF route is exactly the same as the case shown in Figure 2. There will be a total of 25 units of demand lost in the entire planning horizon. Hence, the total cost should be 100 × 1 + 25 × 2 = 150. If we increase the fleet size to two MFs, the optimal routing plan is shown in Figure 4. In this case, all demands can be met so that the penalty cost is zero. The total cost equals to 100 × 2 + 0 × 2 = 200, which is greater than the case with one MF. Therefore, the result of the tradeoff is to have one MF. 10

8:00

8:15

8:30

8:45

9:00

9:15

9:30

9:45

10:00

L1

L1

L1

L1

L1

L1

L1

L1

L1

L2

L2

L2

L2

L2

L2

L2

L2

L2

L3

L3

L3

L3

L3

L3

L3

L3

L3

Figure 4: Route plan for two MFs

3.1. Nominal problem formulation Given the above notations for the parameters and decision variables, we present the nominal formulation (P ) of the MFFSRP by assuming that all parameters are deterministic and known. ∑∑ min f Q + βSpτ (1) s.t.

∑

p

τ

YO0jn ≤ Q

∀j∈L n=t0

YO0jn − ∑

∑

Yjnkq = 0,

Yimjn − YjnSN = 0, Yimjn −

∀i∈L ∀m
∑

(3)

∀j ∈ L, n = t|T |

(4)

( ) Yjnkq = 0,∀j ∈ L, ∀n ∈ N \ t0 ∪ t|T |

(5)

∀k∈L ∀q>n

Yimjn = 0, Xpjτ ≤ apj M Yjmjn , ∑

∀j ∈ L, n = t0

∀k∈L ∀q>n

∀i∈L ∀m
∑

(2)

∀i, j ∈ L, j ̸= i, ∀m, n ∈ T , n > m, and (n − m) ≤ Hij (6) ∀p ∈ P, ∀j ∈ L, ∀τ ∈ T ; m = tτ −1 , (7) and n = tτ

Xpjτ ≤ CYjmjn ,

∀j ∈ L, ∀τ ∈ T ; m = tτ −1 , n = tτ

(8)

Xpjτ + Spτ ≥ Wpτ ,

∀p ∈ P, τ ∈ T

(9)

p

∑ j

11

Q ∈ Z+ Yimjn ∈ Z+ , YO0jn ∈ Z+ , YjnSN ∈ Z+ , Xpjτ ≥ 0, Spτ ≥ 0,

∀i ∈ L, ∀j ∈ L, ∀m ∈ T , ∀n > m ∀j ∈ L, n = t0 ∀j ∈ L, n = t|T | ∀p ∈ P, ∀j ∈ L, ∀τ ∈ T ∀p ∈ P, ∀τ ∈ T

(10) (11) (12) (13) (14) (15)

The objective function (1) intends to minimize the total cost of the system throughout the planning horizon. Constraint (2) restricts that the total number of MFs used cannot exceed the fleet size. Constraints (3) state the conservation restriction for the nodes at time t0 whose inbound arcs can only originate from node O and constraints (4) are similar to constraints (3), except that they are imposed over the nodes at time t|T | when outbound arcs can only enter into the node N . Constraints (5) are the conservation constraints for all nodes except the ones at time t0 and t|T | . Constraints (6) indicate a MF can travel through an arc in the time-space network only if it satisfies the travel time restriction. Note that these constraints are related to the assumption that the travel time must be explicitly accounted for and the MFs are unavailable to provide services when they are in motion. Constraints (7) state the coverage restrictions. Here, M is a sufficient large value, i.e. M = max ∀p∈P {Wpτ } in this model. If either custom p is beyond ∀τ ∈T

the coverage of location j or there is no MF being stationary at location j in period τ , constraints (7) leads to Xpjτ = 0; otherwise, constraints (7) become redundant. Constraints (8) impose an upper limit for the service capacity of the stationary MFs at location j during the time interval [tτ −1 , tτ ). Constraints (9) account for the amount of demand from customer p in each time period that is covered and the amount of demand that is not covered. The remaining constraints preserve nonnegativity for all decision variables in the model. 3.2. Two-stage robust MFFSRP formulation Given that all parameters are known, the nominal model of the MFFSRP can be solved as a mixed-integer linear program (MILP). However, if customer demands are uncertain, the solution obtained by solving the nominal model may not be appropriate. To address the uncertainty factor, we seek to reformulate the MFFSRP by applying the two-stage RO approach, with the customer demand being bounded by a polyhedral uncertainty set. 12

First, we introduce the uncertainty set considered in this work. We define ˜ that the uncertain demand [ Wpτ of each customer ] p in each time period τ lies ˜ ¯ ˆ ¯ ˆ ¯ pτ ≥ 0 represents in an interval, Wpτ ∈ Wpτ − Wpτ , Wpτ + Wpτ , where W ˜ pτ and W ˆ pτ ≥ 0 the maximum deviation. Note that the nominal value of W the values of demands would have a negative impact on solutions only when they are greater than their nominal values in this study, i.e., one should focus ¯ pτ + W ˆ pτ . The polyhedral uncertainty only on the interval between 0 and W set U is shown as follows: { } ˜ ¯ ˆ U = W : Wpτ = Wpτ + Wpτ zpτ , ∀p ∈ P, ∀τ ∈ T , z ∈ Z , { } ∑ ∑∑ where Z = z : zpτ ≤ Γτ , zpτ ≤ Γ, 0 ≤ zpτ ≤ 1, ∀p ∈ P, ∀τ ∈ T , Γτ p

p

τ

denotes the budget of uncertainty for each time period τ and Γ denotes the total budget of uncertainty for the entire planning horizon. Note that if Γ ∈ Z+ , the value of Γ is the maximum number of customers that can simultaneously deviate from their nominal values during the entire planning horizon. For Γ = 0, all demands stay at their nominal value and the problem is equivalent to the nominal one in Section 3.1, whereas Γ = |P| × |T | results in the box uncertainty set (Ben-Tal and Nemirovski, 1999) where all customers have the largest possible amount of demands. For the two-stage RO approach, the decisions are made in two steps: first, the size of the fleet as well as its routing and scheduling plan needs to be fixed, i.e., the first-stage design solution. After that, as the actual demands become known, the allocations of the demands to customers are decided, i.e. the second-stage recourse solution. The purpose of the two-stage RO model is to seek the first-stage design solution which allows the system to obtain its best performance even for the worst possible realization of the uncertainty set. The mathematical formulation of the two-stage robust MFFSRP model, denoted as Prob , is as follows: [Prob ] min f Q + opt [RP(Y )] s.t. Constraints(2) − (6), (10) − (13) where opt [RP(Y )] is the optimal objective function value of the recourse problem, RP(Y ), which is defined as follows: ∑∑ [RP(Y )] max min (16) βSpτ z

X,S

p

τ

13

s.t. Xpjτ ≤ apj M Yjmjn , ∀p ∈ P, ∀j ∈ L, ∀τ ∈ T , and m = tτ −1 , n = tτ ∑ Xpjτ ≤ CYjmjn ,∀j ∈ L, ∀τ ∈ T ,

(17)

p

∑

Xpjτ + Spτ ≥

and m = tτ −1 , n = tτ (18) ¯ pτ + W ˆ pτ zpτ , ∀p ∈ P, τ ∈ T (19) W

zpτ ≤ Γτ ,

∀τ ∈ T

j

∑

(20)

p

∑∑ p

zpτ ≤ Γ

(21)

τ

0 ≤ zpτ ≤ 1, Xpjτ ≥ 0, Spτ ≥ 0,

∀p ∈ P, ∀τ ∈ T ∀p ∈ P, ∀j ∈ L, ∀τ ∈ T ∀p ∈ P, ∀τ ∈ T

(22) (23) (24)

Note that the objective function of the RP(Y ) is in min-max sense, which is not suitable to use in solution method design process. We present next how to reformulate the RP(Y ) into a maximization problem by replacing the ¯ inner minimization problem by its ( dual. ) For fixed Y and z¯, we can have the inner minimization problem RP Y¯ , z¯ . ( ) [RP Y¯ , z¯ ] min X,S

∑∑ p

βSpτ

τ

s.t. Constraints (17), (18),(23) and (24) ∑ ¯ pτ + W ˆ pτ z¯pτ , ∀p ∈ P, τ ∈ T Xpjτ + Spτ ≥ W j

( ) The RP Y¯ , z¯ ( can be)easily decomposed into |T | independent smaller sub-problems, RP) Y¯ , z¯, τ , with respect to each time period τ . We observe ( that RP Y¯ , z¯, τ is in fact a linear transportation problem, and it is feasible( for all theory, we write out the dual of the ) ) the demands. By ( duality RP Y¯ , z¯ , denoted as DRP Y¯ , z¯ : ∑∑ ∑∑∑( ) ( ) (C Y¯jmjn )γjτ [DRP Y¯ , z¯ ] max − apj M Y¯jmjn λpjτ − λ,γ,η

p

j

τ

j

14

τ

+

∑∑( p

s.t.

) ¯ pτ + W ˆ pτ z¯pτ ηpτ W

(25)

τ

− λpjτ − γjτ + ηpτ ≤ 0, ∀p ∈ P, ∀j ∈ L, ∀τ ∈ T ηpτ ≤ β, ∀p ∈ P, ∀τ ∈ T λpjτ ≥ 0, ∀p ∈ P, ∀j ∈ L, ∀τ ∈ T γjτ ≥ 0, ∀j ∈ L, ∀τ ∈ T

(26) (27) (28) (29)

where λ, γ, and η are the multiplies of the constraints (17),(18) and (19), ( ) ¯ respectively. We denote the feasible region of DRP Y , z¯ as Ω. Note that Ω is a bounded polyhedral and the optimal solution v ∗ := (λ∗ , γ ∗ , η ∗ ) is an extreme point of Ω. Then the second stage recourse problem RP(Y ) can be equivalently ( ) transformed into a bilinear maximization program, denoted as DRP Y¯ . ∑∑∑( ∑∑ ( ) ) [DRP Y¯ ] max − apj M Y¯jmjn λpjτ − (C Y¯jmjn )γjτ λ,γ,η,z

p

+

j

τ

∑∑( p

)

j

τ

¯ pτ + W ˆ pτ zpτ ηpτ W

τ

s.t. Constraints (20)-(22) and (26)-(29) Proposition 1. DRP (Y ) is a convex piecewise linear function in Y . Proof. We denote the extreme points set of Ω and U as I and K, respectively. Since both Ω and U are bounded polyhedral, i.e.,( |I|) < ∞ and |K| < ∞, there exists an optimal solution (v ∗ , z ∗ ) for DRP Y¯ such that v ∗ ∈ I and z ∗ ∈ K (Horst and Tuy, 1996). Therefore, we can have: [ ( ) ] DRP (Y ) = max − λTi B1 + γ Ti B2 Y + η Ti W z k v i ∀i∈I z k ,∀k∈K

which shows that DRP (Y ) is the maximum of a finite number of affine functions in Y . The proposition is proved.2 + Proposition 2. If Γ ∈ Z and Γτ ∈ Z+ , ∀τ ∈ T , there exists an optimal ) ( ∗ ∈ {0, 1}, ∀p ∈ P, ∀τ ∈ T . solution (v ∗ , z ∗ ) of DRP Y¯ such that zpτ

Proof. As mentioned (in )the proof of Proposition 1, an optimal solution (v ∗ , z ∗ ) exists for DRP Y¯ while z ∗ is an extreme point of U. When Γ and Γτ are all integers, the extreme points of U can only be either 1 or 0. The proposition is thus proved.2 15

4. Solution methodology Given the two-stage characteristics of the problem, it is natural to attempt to solve the problem with a two-stage algorithm. In this section, we develop a solution approach based on the cutting plane type method used by Thiele et al. (2009), which belongs to the category of the Benders decomposition method (Benders, 1962). An overview of the framework of this cutting plane method is given below. Step 1. Solve the master problem (MP) to obtain the solution to the firststage problem. The lower bound value of the problem is updated in this step. The formulation of the MP is shown as follows, where constraints (30) (also known as optimality cuts or Benders cuts) are deduced on the base of Proposition 1. [MP] min f Q + θ s.t. Constraints (2)-(6) and (10)-(13) ∑∑∑ ∑∑ k ¯ k )Yjmjn − θ≥− (apj M λ (C γ¯jτ )Yjmjn pjτ p

+

j

∑∑( p

τ

)

k k ¯ pτ + W ˆ pτ z¯pτ W η¯pτ ,

j

∀k ≤ s

τ

(30)

τ

Step 2. The obtained first-stage solution is then input into the recourse problem. Because of the complete recourse nature of the RP(Y¯ ) which implies that any possible realization of uncertain demand for a given solution Y¯ must be feasible, ( ) the feasibility checking process is omitted. Solve the DRP Y¯ optimally and update the upper bound value. Check the gap between the lower bound value and upper bound value. If optimality is reached, the algorithm stops; otherwise, go to Step 3. Step 3. Add the newly generated optimality cut to the master problem, and go back to Step 1. Instead of applying the Thiele et al.’s algorithm directly, we attempt to improve this cutting plane method by further analyzing the master problem and the recourse problem. For the former, we attempt to enhance the lower bounding effect of the master problem. For the latter, we aim at developing a hybrid algorithm by considering a combination of both heuristic and exact 16

algorithm. The general idea behind these improvements is to extract more useful information from the recourse problem to tighten the feasible region of the master problem so as to reduce the overall computation time. In the following, we describe the master problem solver and the recourse problem solver separately. 4.1. Master problem solver Since the master problem is obtained by relaxing the second-stage problem, the tightness of the lower bound is the key of convergence efficiency. Here we propose to improve the master problem solver in two aspects: 1) generating multiple optimality cuts rather than a single one in each iteration, and 2) constructing several lower bound inequalities for the master problem by further exploiting the structure of the recourse problem. 4.1.1. Multiple optimality cuts Since the recourse problem can be considered as a number of independent subproblems once the first-stage solutions and the uncertain demands are known, we may construct optimality cuts with respect to each subproblem, which is similar as the multicut version L-shaped method introduced in Birge and Louveaux (1988). With that, the information contained in the optimality cuts is more specific, which may result in better lower bounding effect and faster convergence of the algorithm. ∑ Let θpτ denote the penalty variable for each time period τ , where θ = τ θτ and θτ ≥ 0. Denote sτ as the optimality cut set with respect to each τ , indexed by m. The newly developed optimality cuts are shown as the following inequalities: ∑∑ ∑ k ¯ k )Yjmjn − θτ ≥ − (apj M λ (C γ¯jτ )Yjmjn pjτ p

+

j

∑(

)

k k ¯ pτ + W ˆ pτ z¯pτ W η¯pτ ,

j

∀k ≤ sτ , ∀τ

(31)

p

Note that the original optimality cut is simply the summation of the multiple cuts (31) for all τ ∈ T . Therefore, more or at least an equal amount of information is incorporated into the master problem in each iteration by using (31), compared with the original cutting plane method.

17

4.1.2. Lower bound inequalities Lei et al. (2014) proposed to tighten the lower bound of the master problem by constructing several lower bound inequalities in solving two-stage stochastic programs. The rationale behind this technique is to incorporate part of the second-stage information into the master problem through certain kind of constraints. In this section, we try to apply this technique into our algorithm design, where the lower bound inequalities are uniquely constructed by exploiting the specific characteristics of the recourse problem. Similar to generating multiple optimality cuts in a more specific manner, we attempt to construct the lower bound inequalities with respect to each time period. Proposition 3. The inequalities ( ) ∑ ∑ ¯ pτ − θτ ≥ β W CYjmjn , p

∀τ, m = tτ −1 , n = tτ

(32)

j

are lower bound inequalities for the MFFSRP. Proof. The left hand side of inequalities (32) is the penalty variable for time period τ . Based on the definition of the uncertainty set, the demands ¯ pτ . The of each customer p in period τ is no less than its mean value, i.e., W total amount of demands the MFs can ∑ provide will not exceed the maximum capacity of the entire fleet, i.e., j CYjmjn . Hence, the lowest possible penalty cost in period τ is computed as the right hand side of inequalities (32).2 Proposition 4. The inequalities ( ) ∑ ¯ pτ − θpτ ≥ β W δapj Yjmjn , θτ ≥

∑

∀p, τ, m = tτ −1 , n = tτ

(33)

j

θpτ

(34)

p

are lower bound inequalities for the MFFSRP, where { ¯ pτ ≥ C C if W δ= ¯ pτ otherwise W 18

Proof. We define a new nonnegative continuous variable, θpτ , as the penalty variable for customer p in each time period. The inequality (34) is thus true, because the total penalty cost in one time period should be no less than the summation of the penalty cost from all customers. Similar to Proposition 3, we know that the lowest possible amount of demands for customer p in period τ is its mean value. The upper limit of amount of demands the MF ∑ fleet can possibly provide to customer p is no greater than j Capj Yjmjn if ¯ pτ ≥ C, or ∑ W ¯ pτ apj Yjmjn if W ¯ pτ < C. Therefore, the penalty cost for W j customer p in time period τ is greater than or at least equal to the right hand side of inequalities (33), which justifies the correctness of the lower bound inequalities (33) and (34).2 Note that the lower bound inequalities (32)-(34) are appealing because they can improve the quality of the lower bound of the master problem with no significant increase in complexity and no influence on the optimality condition of the original problem. This feature is especially desirable when the master problem is nontrivial to solve. After making the above two enhancements, we can present the new formulation of the master problem, denoted as MP′ . ∑ [MP′ ] min f Q + θτ τ

s.t. Constraints (2)-(6), (10)-(13), (31)-(34) 4.2. Recourse problem solver ( ) The recourse problem DRP Y¯ is a bilinear maximization program with linear constraints, which belongs to a class of convex maximization problem proven to be NP-hard (Vavasis, 1991). According to Gabrel et al. (2014a), it is usually the most computationally intensive part in the entire solving process. Both exact methods (Gabrel et al., 2014a; Thiele et al., 2009) and heuristic methods (Bertsimas et al., 2013) have been developed before. The exact methods emphasize on the optimality of the solution but usually take more computational efforts. On the other hand, the heuristic methods aim at obtaining some feasible solution in a short time, whereas the solution quality cannot be guaranteed. In this study, we attempt to develop a hybrid algorithm by taking advantage of both types of algorithms. The main framework of the method is illustrated in Figure 5. At first, we try to come up with some feasible solutions of the recourse problem in a fast way by applying a certain heuristic 19

method. If these solutions indicate that the current first-stage solution is not optimal, we will then add optimality cuts to the master problem and the exact solver can be skipped; otherwise, we resort to an exact solution method (e.g., branch-and-bound algorithm) to warrant the global optimality of the first-stage solution. Note that the recourse problem solver would incur the exact solver at least once to finish the algorithm. Hence, it is essentially an exact algorithm such that the optimality can be guaranteed. Yet, through the steps in the dash line box in Figure 5, we are able to quickly identify and reject many inferior first-stage solutions. The computation time thus can be saved whenever the relatively slower exact solver is skipped. Heuristic solver

Test first-stage solution optimality

No

Generate cut, go back to master problem solver

Pass Exact solver

No Test first-stage solution optimality

Pass Optimality is reached, stop.

Figure 5: Framework of the recourse problem solver

20

4.2.1. Heuristic method For the heuristic method, we apply the outer approximation (OA) algorithm to solve the second-stage recourse problem. A detailed introduction of the OA algorithm can be found in Duran and Grossmann (1986); Fletcher and Leyffer (1994). Bertsimas et al. (2013) first apply the OA method in solving the bilinear program for the robust optimization problem. The key of the OA algorithm is to linearize the bilinear term in the objective function around intermediate solution points, and then substitute the linear approx( ) ¯ imation for the bilinear term. Since the DRP Y is non-concave and the approximation of the bilinear term is performed in a heuristic way, only a local optimum is guaranteed. ˆ z is linearized at a known point (η k , z¯k ) as follows: The bilinear term η T W ˆ zk + ηT W ˆ (z − z k ) + η T W ˆ zk . Lk (η, z) = (η − η k )T W k k ( ) ( ) Here we present an approximation to DRP Y¯ , denoted as L-DRP Y¯ , z k , in a compact form shown as follows. ( ) ( ) ¯z+∆ [L-DRP Y¯ , z k ] max − λT B1 + γ T B2 Y¯ + eT W s.t. Ll (η, z) ≥ ∆, ∀l = 1, . . . , k −λ−γ+η ≤0 λ ≥ 0, γ ≥ 0, η ≤ β, z ∈ U An initial value of z¯1 ∈ U has to be set before the start of the OA algorithm. In this paper, this initial value is computed by solving the following linear program, denoted by Z-INIT. ∑∑ ˆ pτ zpτ [Z-INIT] max W p

τ

s.t. Constraints (20)-(22) Now, we present a framework of the OA algorithm. ) ( Step 1. Solve DRP Y¯ , z k optimally, get the optimal solution v ∗k := (λ∗k , γ ∗k , η ∗k ), and update ( ) the lower bound value of the dual of the recourse problem DRP Y¯ . ) ( Step 2. Linearize the bilinear term and solve L-DRP Y¯ , z k (to )obtain a new z¯ ∈ U. Also, after the upper bound value of DRP Y¯ is updated, go to Step 3. 21

Step 3. Check the optimality criterion by comparing lower bound value and the upper bound value. If optimality is reached, the algorithm terminates; otherwise, go to Step 1. 4.2.2. Exact method In this section, we reformulate the recourse problem as a mixed-integer program to enable the use of off-shelf MILP solver, e.g., CPLEX. ( ) We replace the bilinear term ηpτ zpτ in the objective function of DRP Y¯ with a new variable ζpτ by forcing ζpτ = ηpτ if zpτ = 1, and ζpτ = 0 otherwise, as we know zpτ is binary variable. The( mixed-integer reformulation of the recourse ) ¯ problem, denoted as MIP-DRP Y , can be expressed as follows: ∑∑∑( ∑∑ ( ) ) apj M Y¯jmjn λpjτ − [MIP-DRP Y¯ ] max − (C Y¯jmjn )γjτ λ,γ,η,z

p

+

j

∑∑ p

τ

¯ pτ ηpτ + W

τ

∑∑ p

j

ˆ pτ ζpτ W

τ

(35)

τ

s.t. Constraints (20)-(22) and (26)-(29) ζpτ ≤ ηpτ , ∀p ∈ P, ∀τ ∈ T ζpτ ≤ βzpτ , ∀p ∈ P, ∀τ ∈ T zpτ ∈ {0, 1}, ∀p ∈ P, ∀τ ∈ T ζpτ ≥ 0, ∀p ∈ P, ∀τ ∈ T

(36) (37) (38) (39)

Constraints (37) are originally written as ζpτ ≤ M zpτ for each customer and each time period, where M is a sufficiently large constant representing the bound of ηpτ . Since we already have ηpτ ≤ β as in (27) and zpτ ∈ {0, 1}, it is natural for us to set the value of M to β as in (37). To further improve the computational performance, we make several additional adjustments to the solution process. First, we apply the best solution ( ) obtained from the OA algorithm as the initial solution to the MIP-DRP Y¯ . With this, some of the nodes in the branch-and-bound tree are directly pruned. Second, instead of solving the model to optimality without considering the computation time, we place a time threshold to the algorithm. If the algorithm reaches the optimality within the time threshold, the optimal solution is obtained and the algorithm stops; otherwise, the algorithm will stop once it finds a feasible solution whose objective function value is ¯ greater than the current value of θ. Another benefit of the exact method ( ) is that we can always track the upper bound value of the MIP-DRP Y¯ , which enables us to update the 22

( ) ¯ upper bound value of the MFFSRP. ( ) Denote ( )U BRP Y as the best upper bound value of the MIP-DRP Y¯ and F Y¯ as the fixed costs of the fleet for the first-stage solution Y¯ . Thus, we can have an upper bound for the MFFSRP as described in Proposition 5. ( ) ( ) Proposition 5. F Y¯ + U BRP Y¯ is an upper bound for the MFFSRP. ( ) Proof. Assume that we solve the MIP-DRP Y¯ (to optimality ) ( ) and obtain its optimal objective function value θ∗ ((Y¯ ).) Then F Y¯ +θ∗ Y¯ is an upper bound for the MFFSRP. Since U BRP Y¯ is the worst penalty cost ( )possible ∗ ¯ ¯ ¯ value for( the ) U BRP Y ≥(θ ()Y ) and obvi( ) ) first-stage ( solution ) ( Y) , we (have ously F Y¯ +U BRP Y¯ ≥ F Y¯ +θ∗ Y¯ . Therefore, F Y¯ +U BRP Y¯ is definitely an upper bound for the MFFSRP.2 4.3. The enhanced cutting plane method By integrating the above master problem solver and the recourse problem solver into the framework of the Thiele et al.’s algorithm, we have obtained the enhanced version of the cutting plane method, denoted as E-CP. It can be characterized as a two-level algorithm, in which the outer level includes the master problem solver and the inner level contains the recourse problem solver. The detailed procedures of the algorithm are presented in the Appendix A, where the outer level algorithm is shown in Algorithm 1, inner level algorithm in Algorithm 2, and the OA algorithm in Algorithm 3. 5. Numerical study The primary purposes of the numerical study are to gain some insights into: 1) the computational efficiency of the E-CP method in comparison with the traditional CP method, 2) the effectiveness of the two-stage RO approach. All of the experiments are run on a PC with 2.6GHz CPU and 16GB RAM. The CPLEX 12.5 optimization package has been used as the MILP solver for the master and recourse problems. 5.1. Parameters generation In this section, we introduce the relevant parameters and data in constructing the numerical experiments. The candidate locations and customer nodes uniformly spread over a L × W rectangular region. The distances between nodes are measured in Euclidean sense, and without loss of generality, 23

they are normalized into values between 0 and 1. The coverage distance is chosen as the minimum distance which can guarantee that each customer node is covered by at least one candidate location. The nominal value of uncertain demand for each customer in each period is randomly generated from a uniform distribution U [w, w], where 0 ≤ w ≤ ˆ pτ is given as a product of mean value W ¯ pτ with a w. The deviation W ˆ pτ = ϵτ W ¯ pτ . Since we attempt to generate a solution for coefficient ϵτ , i.e. W a whole planning horizon at a given point, we need to make the prediction over uncertain demands for the entire time periods following the decision making point. It is reasonable to assume that better prediction can be made for the demands at time periods closer to the decision making point than those in the distant future. Therefore, we let ϵτ = µϵτ −1 with ϵ0 being the minimum level of deviation and µ ∈ [1.0, 1.5], which means that ϵτ is a non-decreasing function in τ . The budget of uncertainty for each time period Γτ takes values in the range of 0 to |P|, and the budget of uncertainty for entire planning ∑ horizon Γ is given as a fraction of the summation of all Γτ , i.e., Γ = ρ τ Γτ and ρ ∈ (0, 1]. In the experiments, we set different values of the total uncertainty budget Γ = ξ|P||T |. ξ is a fraction ranging from 0 to 1, and we name it as the uncertainty level in the following experiments. Then the value of Γτ for each time period is randomly generated within the range of 0 to |P|, ∑ while satisfying the restriction Γ = ρ τ Γτ . Note that both Γ and Γτ are restricted to be integer values throughout the numerical studies. The values of the fixed cost factor f , the penalty cost factor per unit demand β, the capacity C for an MF, and other parameters are given as in Table 2. Without loss of generality, we assume that all of the cost parameters are calculated in terms of present monetary value in this paper. The values of parameters mentioned above are listed in Table 2. Table 2: Parameter values in computational experiments

Parameter Value

L

W

f

β

C

ρ

ϵ0

w

w

100

100

100

1.0

100

1.0

0.1

50

100

The test instances are generated with different number of customers (|P| ∈ [5, 10, 15, 20, 25, 30]) and time periods (|T | ∈ [10, 20]). The number of candidate locations |L| equals to the cardinality of the customer nodes set, i.e. |L| = |P|. Therefore, a total of 12 different test cases are generated. For 24

each test case, 10 instances are independently generated and solved. The average values are computed and have been employed as the performance measurement of the algorithm. We use the combination (|P|,|L|,|T |) as the code marks for different test sets. For example, (10, 5, 20) represents the instance with 10 customers, 5 candidate locations and 20 time periods. The scale of the instance is then determined as the product of the number of customers, the number of candidate locations, and the number of time periods (e.g. the scale of test case (10, 10, 20) is 2000). The convergence tolerance and the time limit for the outer level algorithm, inner level algorithm and the master problem are shown in Table 3. Table 3: Tolerance gap and computation time limit

Scale 0-500 500-100 1000-2000 2000-4000 >4000

Outer level algorithm Gap(%) Time (s) 0.01 100 0.1 300 1.0 500 2.0 1000 3.0 3000

Inner level algorithm Gap(%) Time (s) 0.001 5 0.01 5 0.1 10 0.2 20 1.0 30

Master problem Gap(%) Time (s) 0.005 5 0.05 5 0.5 10 1.5 20 2.5 30

5.2. Computational efficiency In this section, we aim at quantifying the benefits of the proposed method in three aspects: 1) the improvement in computation time of the modified multiple cuts over the single cut method, 2) the effect of the lower bound inequalities on the convergence process, and 3) the efficiency of the hybrid solver in solving the recourse problem. The value of the budget level of robustness Γ is set to be 0.5 for all of the test instances in this section. The other parameters are given as described in Section 5.1. 5.2.1. Efficiency of the lower bound inequalities To study the effect of the lower bound inequalities (32)-(34) on the performance of the algorithm, we separately applied the E-CP method to the test cases with and without the lower bound inequalities, and then compare the performance of the algorithm in terms of the computation time and gap values. The results are displayed in Table 4. The columns tagged with “After” and “Before” refer to the results obtained by using the typical E-CP method and the algorithm without the lower bound inequalities, respectively. 25

This is done by keeping all other components the same as the E-CP method. We observe that by including the lower bound inequalities, the computation time is substantially reduced and the quality of the gap values improved. In comparison, without the lower bound inequalities, the performance of the algorithm is very likely to deteriorate, especially when the scale of the problem becomes relatively large. For example, when the scale of instance is increased to 10000, the algorithm without lower bound inequalities can barely improve the feasible solution in 3000 seconds, whereas the one including lower bound inequalities is able to reach a tolerance gap of 3% in less than 1500 seconds. Table 4: Computational results before and after adding lower bound inequalities

Scale 250 500 1000 2000 2250 4000 4500 6250 8000 9000 12500 18000

Test case |P| |L| |T | 5 5 10 5 5 20 10 10 10 10 10 20 15 15 10 20 20 10 15 15 20 25 25 10 20 20 20 30 30 10 25 25 20 30 30 20

CPU time (s) Before After 0.52 0.28 0.96 0.85 17.15 6.90 51.91 41.42 88.24 48.74 208.87 81.39 460.87 177.58 974.93 276.32 2340.43 414.39 2737.45 595.83 3000.00 1223.60 3000.00 1582.85

Gap (%) Before After 0.00 0.00 0.00 0.00 0.00 0.00 0.33 0.33 0.67 0.33 0.75 0.75 18.83 1.83 80.40 2.20 61.00 2.60 83.67 3.00 100.00 3.00 100.00 3.00

To get a better understanding of how the lower bound inequalities affect the quality of the lower bound, we consider one of the test instances with 20 customer points, 20 candidate locations and 10 time periods to show the difference between the lower bound values throughout the solving process with and without the inequalities. As shown in Figure 6, it is quite obvious that the lower bounds converge faster when the lower bound inequalities (32)(34) are introduced. Moreover, because of the better lower bounding effect, the algorithm may converge to the optimal in less number of iterations which in turn reduce the computational efforts. For instance, in Figure 6, we can see that the algorithm converges to the optimum in 50 iterations if the lower bound inequalities are added, whereas it takes 72 iterations to reach the same lower bound value if we do not consider them. 26

Lower bound value

1200 1100 1000 900 800 700 600 500 400 300 200 5

10

15

20

25

30

35 40 Iteration

45

50

55

60

65

70

Without lower bound inequalities With lower bound inequalities

Figure 6: Comparisons of lower bound value before and after the inequalities added

5.2.2. Efficiency of the multiple optimality cuts To compare the performance of multiple cuts with the single cut, we test the single cut version first and then apply the multicut version on the same test cases. The computation time and gap values are presented in Table 5. We observe that the multicut version is significantly better than the single cut version with respect to either computation time or solution quality. The single cut version fails to reach the tolerance gap within the time limit for

Table 5: Computational results for single cut version and multicut version

Scale 250 500 1000 2000 2250 4000 4500 6250 8000 9000 12500 18000

Test case |P| |L| |T | 5 5 10 5 5 20 10 10 10 10 10 20 15 15 10 20 20 10 15 15 20 25 25 10 20 20 20 30 30 10 25 25 20 30 30 20

CPU time (s) Single Multiple 0.91 0.28 3.56 0.85 173.89 6.90 500.00 41.42 500.00 48.74 1000.00 81.39 1000.00 177.58 1000.00 276.32 3000.00 414.39 3000.00 595.83 3000.00 1223.60 3000.00 1582.85

27

Gap (%) Single Multiple 0.00 0.00 0.00 0.00 68.00 0.00 100.00 0.33 100.00 0.33 100.00 0.75 100.00 1.83 100.00 2.20 100.00 2.60 100.00 3.00 100.00 3.00 100.00 3.00

1200

Lower bound value

1150 1100 1050 1000 950 900 850 800 5

10

15

20

25 Iteration

Multiple cuts

30

35

40

45

50

Single cut

Figure 7: Comparisons of lower bound value for single cut version and multicut version

most cases when the scale is greater than 1000. On the contrary, the multicut version exhibits to be quite efficient to obtain a good enough solution. We believe this improvement is the result of better lower bounding effect of the multiple optimality cuts. Take the same test instance in Section 5.2.1 as an example. We show the lower bound values throughout the solving process for both versions in Figure 7. It is easy to see that the single cut version performs rather poorly in increasing the lower bound. The results clearly demonstrate the importance of containing more specific information in the optimality cuts. 5.2.3. Efficiency of the hybrid recourse problem solver In this section, we analyze the performance of the hybrid algorithm by comparing it with the method which solves the MILP of the recourse problem to its optimality. We conduct experiments for different test cases for both methods, and the computational results are shown in Table 6. The results clearly indicate a distinct advantage of the hybrid method in solving the recourse problem in terms of effectiveness and efficiency. As shown in Table 6, the computation times in solving recourse problem, as well as the total computation times, are consistently reduced by applying the proposed hybrid algorithm for all test cases. When the scales are small, the two methods are comparable. However, along with the increase in the scale of the problem, the computation time of the exact method significantly increases. In addition, the gap values by using the exact approach are quite large for some of the test cases. Moreover, the exact algorithm fails to con28

Table 6: Computational results for the hybrid method and the exact method

Scale 250 500 1000 2000 2250 4000 4500 6250 8000 9000 12500 18000

Test case |P| |L| |T | 5 5 10 5 5 20 10 10 10 10 10 20 15 15 10 20 20 10 15 15 20 25 25 10 20 20 20 30 30 10 25 25 20 30 30 20

Total time (s) Hybrid Exact 0.3 0.4 0.9 0.8 6.9 9.5 41 403 49 430 81 472 178 1141 276 1238 414 3000 596 3000 1224 3000 1583 3000

Time to solve RP(s) Hybrid Exact 0.2 0.3 0.6 0.6 4.3 7.8 38 401 42 425 61 455 153 1109 151 1232 255 2581 242 2763 326 2904 348 2970

Gap (%) Hybrid Exact 0.0 0.0 0.0 0.0 0.0 0.0 0.3 2.0 0.3 8.0 0.8 14.8 1.8 26.3 2.2 26.8 2.6 28.6 3.0 30.3 3.0 42.7 3.0 44.6

verge to the preset tolerance gap within the given time limit, which does not happen by applying the hybrid method. 5.3. Solution effectiveness In this section, we focus on the numerical studies on the solution quality from the two-stage RO model with polyhedral uncertainty set against the nominal and box uncertainty case. We compare the performance of these three approaches in four aspects: 1) economic effectiveness of the first-stage design solutions, 2) the reliability of the first-stage solutions, 3) the robustness of the first-stage solutions over different demand distributions, and 4) the robustness of the first-stage solutions for different choice of parameters. The tests are performed mainly in two steps: 1) solve the two-stage RO model for different values of the total uncertainty budget to obtain the first-stage design solutions, then 2) for each first-stage design solution, solve the recourse problem for different sets of 500 randomly generated demand scenarios to obtain the penalty costs and percentage of unmet demands. Two sets of demand scenarios are generated. In the first set, demand scenarios are generated following a uniform distribution in the interval ] [ ¯ pτ − W ˆ pτ , W ¯ pτ + W ˆ pτ for each customer p at time τ . In the second set, W ¯ pτ and standard normally distributed demands are generated with mean W ˆ pτ /1.44, which makes an 85 percentile of demands fall in the deviation W 29

[ ] ¯ pτ − W ˆ pτ , W ¯ pτ + W ˆ pτ (Bertsimas et al., 2013). Since the values range of W of demands should not be negative, sampled demands with negative values are truncated. Note that the two randomly generated data sets are not necessarily bounded by the budgeted uncertainty set, which allow us to test the “robustness” of the proposed model when the uncertainty is inaccurately defined. We have conducted experiments for multiple instances with different size to support the conclusions in the following part of Section 5.3. Considering the fact that the results obtained for different instances consistently share some common features, we decide to illustrate the conclusion of the experiments with one representative instance (test case (10, 10, 10)) to keep our presentation simple. The results for the other test instances are included in the supplementary file for reference. Moreover, since the computational results for uniformly and normally distributed demands are very similar to each other, we only present the results of the uniformly distributed demands except in Section 5.3.4, where the results from both distributions are presented for comparison. 5.3.1. Economic effectiveness To compare the economic effectiveness of the first-stage design solutions of the two-stage RO model with the nominal and box uncertainty case, we present the performance results for two demand scenario sets when the budget of uncertainty set varies in different levels. We take the test case (10, 10, 10) as an example and the results are shown in Table 7. The first column is the value of budget level ξ. In particular, when ξ = 0, the two-stage RO model is equivalent to the nominal model. Similarly, the two-stage RO model is equivalent to the box uncertainty model for ξ = 1. In Table 7, we have displayed the number of MFs used, the average penalty costs, the average percentage of unmet demands, and the total costs (from column 2 to column 5). The last two columns present the quantity of relative improvement of the two-stage RO solution with respect to the nominal and box uncertainty case, denoted as rn and rbox : rn =

On − Oro Obox − Oro and rbox = , On Obox

where Oro is the summation of the fixed cost and the average penalty cost obtained from the test of a randomly generated data set for the two-stage 30

Table 7: Test results for uniformly distributed demands (ϵ0 = 0.1 )

Budget level (ξ) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 a

MF

Penalty cost

8 8 8 9 9 9 10 11 11 12 12

482.96 473.41 457.19 317.03 305.48 303.73 205.51 114.76 118.02 80.87 80.87

Unmet demand(%) 8.09 7.90 7.57 5.55 4.15 3.80 3.52 2.77 2.59 1.10 1.10

Total cost

rn (%)

rbox (%)

1282.96 1273.41 1257.19 1217.03 1205.48 1203.73 1205.51 1214.76 1218.02 1280.87 1280.87

0.00 0.74 2.01 5.14 6.04 6.18 6.04 5.32 5.06 0.16 0.16

-a 0.58 1.85 4.98 5.89 6.02 5.88 5.16 4.91 0.00 0.00

The ratio is displayed as “-” if it is negative.

RO solution. Similarly, we denote On and Obox as the corresponding costs for nominal and box uncertainty case. As expected, the number of MFs increases along with the increase in the budget level of the uncertainty set while the penalty cost and the percentage of unmet demands decrease, since with more MFs incorporated into the system, the supply capacity becomes larger so that more demands can be covered. However, we observe that the two-stage RO approach improves the system performance against either the nominal or box uncertainty case for the budget level ranging from [0.1, 0.9]. Also, the total costs reach the lowest value when the budget level is around 0.5. Take ξ = 0.5 as an example. The total cost decreases by 6.18% (6.02%) with respect to the nominal case (the box uncertainty case). Meanwhile, the average penalty costs decrease by 37.26% for the nominal case and the unmet demands drop from 8.09% to 3.80%. Therefore, the first-stage design solutions obtained from the two-stage RO model are able to yield more favorable cost efficiency while retaining the conservativeness of the model in a relatively low level. Note that the increase of the budget level and the number of MFs are not linear. For example, when the budget level is increased from 0.0 to 0.2, MF stayed in 8 but the penalty cost and unserved unmet demand ratios become smaller. The reason for the number of MFs staying the same is that the change in demand is not large enough to justify the cost for adding a new 31

MF. However, the routing plan is improved to accommodate the increase in uncertainty. This reveals that the MFFSRP has two ways to handle the robustness: one is to increase the fleet size and the other is to adjust the routing plan. The fleet size will increase only if the adjustment of the routing plan is inadequate in dealing with the uncertainty. This feature is desirable in practice, since the expansion of the fleet should be considered only if the same goal cannot be accomplished by better managing the available resources. 5.3.2. Reliability In this section, we examine how the solution generated by the two-stage RO model may affect the volatility of the performance of the recourse operations. Here we use the same test case as in Section 5.3.1 to illustrate the point. In Table 8, we present the mean and standard deviation of the penalty cost and unmet demands ratio for uniformly distributed demands, respectively. As shown in Table 8, with the increase of the budget of uncertainty set, there are significant reductions in the means and standard deviations for both penalty cost and unmet demand ratio. For the box uncertainty case, i.e., ξ = 1, the solution appears to be the most reliable but at the expense of a reduction in cost efficiency. If we want to avoid being overly conservative, we may choose lower level of uncertainty, such as ξ = 0.5. As a result,

Table 8: Mean and standard deviation of the penalty cost and unmet demand ratio (ϵ0 = 0.1)

Budget level (ξ) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Penalty cost Mean 482.96 473.41 457.19 317.03 305.48 303.73 205.51 114.76 118.02 80.87 80.87

Std 56.45 56.17 56.29 52.81 50.86 47.38 50.30 47.18 44.75 35.97 35.97

32

Unmet demand(%) Mean Std 8.09 1.10 7.90 1.07 7.57 1.07 5.55 1.02 4.15 0.98 3.80 0.94 3.52 0.99 2.77 0.94 2.59 0.89 1.10 0.78 1.10 0.78

the standard deviation of the unmet demand ratio by 15% compared to the nominal case. That means the reliability of the system performance increases by around 15% if we adopt the two-stage RO approach. Moreover, for ξ = 0.5, we observe that the unmet demand ratios mostly lie within a range between 2.86% to 4.74%, which is less than 5% and reduced by almost half compared to the nominal case (i.e., 6.99% to 9.19%). Given the fact that the entire system has reached its lowest total cost when ξ = 0.5, we believe the unmet demand ratio (less than 5%) is acceptable for most reallife commercial businesses. 5.3.3. Robustness against demand variations We also test the performance of the two-stage RO model for different levels of demand variations using the same test case in the last two sections. In Figure 8, we show the results for the total cost, penalty cost and the unmet demand ratio for three different minimum demand variation values, ϵ0 ∈ {0.1, 0.3, 0.5}. From Figure 8a, we can see that the two-stage RO approach retains its advantage of cost efficiency against both nominal and box uncertainty case. In particular, in comparison with the nominal case, the performance of the two-stage RO model shows to be more stable. For example, the total costs for the three variations range from 1200 to 1600 when the budget level is around 0.5, whereas they vary between 1200 and 2100 for the nominal case. From Figure 8b and 8c, we can see that the average penalty costs and unmet demand ratios are quite close to each other when the budget level is greater than 0.2, which is not the case for the nominal model. These corresponding results reveal the advantage of the two-stage RO approach in retaining the robustness for different levels of demand variation. The total cost reduction ratios for the three different levels of demand variations are given in Table 9. First, the best cost efficiency is attained when the budget level is somewhere between 0.4 and 0.5 for three different values of ϵ0 , which is consistent with the results shown in Section 5.3.1. Second, the improvement in cost efficiency of the two-stage RO model becomes more significant for higher levels of demand variations with respect to both the nominal and box uncertainty case. Here we take ξ = 0.5 as an example. The total cost decreases by 6.18% compared with the result from the nominal model when ϵ0 = 0.1. A reduction of 20% and 27.66% is expected when ϵ0 = 0.3 and ϵ0 = 0.5, respectively, indicating that the two-stage RO model is even more attractive at a higher level of demand variations. Another issue worth noticing is that rbox is negative for ξ = 0.1 when ϵ0 = 0.3 and for 33

2200

ε0 = 0.1 ε0 = 0.3 ε0 = 0.5

2000

Total cost

1800 1600 1400 1200 1000 0

0.1

0.2

0.3

0.4

0.5 0.6 Budget level

0.7

0.8

0.9

1

(a) Total cost 1400

ε0 = 0.1 ε0 = 0.3 ε0 = 0.5

1200

Total cost

1000 800 600 400 200 0 0

0.1

0.2

0.3

0.4

0.5 0.6 Budget level

0.7

0.8

0.9

1

(b) Average penalty costs ε0 = 0.1 ε0 = 0.3 ε0 = 0.5

Unmet demand ratio (%)

20

15

10

5

0 0

0.1

0.2

0.3

0.4

0.5 0.6 Budget level

0.7

0.8

0.9

1

(c) Average unmet demand ratios Figure 8: Comparison of the results for different level of demand variations

34

Table 9: Comparison of the two-stage RO solution with nominal and box uncertainty case for different demand variations

Budget level (ξ) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

ϵ0 = 0.1 0.00 0.74 2.01 5.14 6.04 6.18 6.04 5.32 5.06 0.16 0.16

rn (%) ϵ0 = 0.3 0.00 8.35 17.53 18.88 18.78 19.57 18.03 16.29 15.98 12.24 11.98

ϵ0 = 0.5 0.00 10.89 23.57 27.62 28.10 27.66 27.52 27.08 25.05 21.38 21.29

ϵ0 = 0.1 0.58 1.85 4.98 5.89 6.02 5.88 5.16 4.91 0.00 0.00

rbox (%) ϵ0 = 0.3 6.30 7.83 7.72 8.62 6.87 5.58 4.54 0.10 0.00

ϵ0 = 0.5 8.04 8.65 8.09 7.91 7.36 4.78 0.12 0.00

ξ = 0.1 and 0.2 when ϵ0 = 0.5. It indicates that when the demand variation level becomes higher, one should be more conservative in choosing the budget level to provide reliable services. 5.3.4. Robustness against demand distributions We illustrate the average penalty costs and unmet demand ratios for uniformly and normally distributed demands in Figure 9a and 9b, respectively. It shows that when budget of uncertainty is greater than 0, they are not sensitive to the underlying distribution. The relative difference of penalty cost (unmet demand ratio) for the two distributions is within the range of 0.35% to 4% (0.39% to 3%) for the budget level from 0.1 to 1. However, for the nominal case, the relative difference is much higher (9.8% for penalty costs and 11.5% for unmet demand ratio). The results indicate that the two-stage RO modeling approach has more stable performance over different demand distributions compared to the nominal one. This property is highly desirable, since it is not easy to accurately characterize the probability distribution of the uncertainty in demand in the real world. 6. Conclusions The primary goal of this paper is to construct a two-stage RO model of the MFFSRP and solve it in an effective and efficient way. Compared 35

Uniform distribution Normal distribution

500

Penalty cost

400

300

200

100 0

0.1

0.2

0.3

0.4

0.5 0.6 Budget level

0.7

0.8

0.9

1

0.9

1

(a) Average penalty costs 10 Uniform distribution Normal distribution

Unmet demand ratio (%)

9 8 7 6 5 4 3 2 1 0

0.1

0.2

0.3

0.4

0.5 0.6 Budget level

0.7

0.8

(b) Average unmet demand ratio Figure 9: Comparison of the results for uniformly and normally distributed demands

to the classic box uncertainty approach, the two-stage RO approach offers more flexibility to the decision maker in dealing with uncertainty for two reasons: 1) the polyhedral uncertainty set is less conservative than the box uncertainty approach since it only considers a limited number of demands deviating from the nominal values, and 2) the second-stage decision variables can be determined adaptively once uncertainty is revealed. However, the two-stage RO model turns out to be computationally very demanding, especially for the recourse problem. To address this drawback, we have made a few ad-hoc enhancements with a proposed two-level cutting plane type method, including 1) adding multiple optimality cuts instead of a single one in each iteration, 2) designing several new lower bound inequalities by utilizing the special structure of the recourse problem and 3) developing 36

a hybrid algorithm by combining an outer approximation algorithm with an exact method to solve the recourse problem. Tests have been conducted to examine the performance of the proposed model. Here are the key findings from the result. First, the proposed E-CP method, with the enhancements made in this work, is shown to be a more efficient solution approach considering a wide range of measures ranging from reliability to economic effectiveness. In particular, the hybrid algorithm for solving the recourse problem can be utilized as a general solution approach to solving other two-stage RO problems. Second, the two-stage RO model of the MFFSRP is shown to be more cost-effective than the nominal and the box uncertainty approach yet flexible enough to maintain a desired level of robustness. Third, in comparison with the nominal case, the two-stage RO modeling approach exhibits a higher level of stability with respect to the magnitude of demand variations or different distributions. Even though the solutions obtained from the box uncertainty approach tends to be the most stable one as expected, they are likely to suffer from poor system-wide performance. The flexibility and reliability of the two-stage RO approach are two very desirable properties in dealing with uncertainty in real-life operations, since, realistically speaking, it is in general very difficult to accurately characterize the uncertain parameters in the first place.

37

Appendix A. Pseudo-code of the E-CP method Algorithm 1 Outer level algorithm Input: Necessary parameters Output: Y¯opt 1: initialization: ν ← 0, LB 0 ← −∞, U B 0 ← +∞ 2: while U B ν − LB ν > ϵ do 3: ν ←ν+1 ( ) 4: solve the MPν to ϵM P -optimality, get the best feasible solution Y¯ ν , θ¯ν , ν ν get the best lower bound value LBM P and the fixed cost value F ν ν−1 5: if LBM P ≥ LB then ν ν 6: LB ← LBM P ◃ Update the lower bound. 7: else 8: LB ν ← LB ν−1 9: end if 10: solve the RP(Y¯ ν ) by applying Algorithm 2, get the set of Benders cuts B and the upper bound value of the recourse problem U BRP 11: if B = ∅ then 12: if (F ν + U BRP ) ≤ U B ν−1 then 13: U B ν ← (F ν + U BRP ) ◃ Update the upper bound. 14: Y¯opt ← Y¯ ν ◃ Update the best solution. 15: else 16: U B ν ← U B ν−1 17: end if 18: else 19: add the newly generated cuts to the MP 20: end if 21: end while

38

Algorithm 2 Inner level algorithm

( ) Input: Y¯ ν , θ¯ν : the solution of the MP Output: set of Benders cuts B and the upper bound value of the recourse problem U BRP ¯ 1: solve DRP(Y ( ) ) using Algorithm 3, get the best solution of the OA algorithm ˆ λ, ˆ γˆ , ηˆ, zˆ θ, 2: if θˆ ≥ θ¯(then ) ˆ γˆ , ηˆ, zˆ to generate the Benders cuts, U BRP ← +∞, return to 3: use λ, Algorithm 1 4: else ( ) ˆ 5: invoke the exact solver, use λ, γˆ , ηˆ, zˆ as the start solution 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:

while ϵRP -optimality is not reached do keep searching the branch-and-bound ( )tree ˜ ˜ if a new feasible solution θ, λ, γ˜ , η˜, z˜ is found then ( ) ( ) ˆ λ, ˆ γˆ , ηˆ, zˆ ← θ, ˜ λ, ˜ γ˜ , η˜, z˜ θ, update U B as the best upper bound value of the branch-and-bound tree if θˆ > θ¯ and CPU time > TRP then break end if end if end while if θˆ ≥ θ¯(then ) ˆ γˆ , ηˆ, zˆ to generate the Benders cuts, U BRP ← U B, return to use λ, Algorithm 1 end if end if

39

Algorithm 3 OA algorithm

( ) Input: Y¯ ν , θ¯ν ( ) Output: solution λ, γ, η, z 1: initialization:κ ← 1, LBbest ← −∞, U Bbest ← +∞ 2: solve the problem Z-INIT, get z¯κ 3: while U Bbest − LBbest > ϵsub do ) ( κ 4: solve the DRP(Y¯ , z¯κ ), get optimal solution λ , γ κ , η κ and objective function value S¯κ 5: if LBbest ≤ S¯κ then 6: LBbest ← S¯κ ◃ Update the lower bound. 7: end if 8: solve the L-DRP(Y¯ ), get z˜ and the objective function value S˜κ 9: if U Bbest ≥ S˜κ then 10: U Bbest ← S˜κ ◃ Update the upper bound. 11: end if 12: κ←κ+1 13: z¯κ ← z˜ 14: end while

Acknowledgements This research is supported by National Natural Science Foundation of China (Grant No. 71272030). The support is gratefully acknowledged. The comments from two anonymous referees are greatly appreciated. References Agra, A., Christiansen, M., Figueiredo, R., Hvattum, L. M., Poss, M., Requejo, C., 2013. The robust vehicle routing problem with time windows. Computers & Operations Research 40 (3), 856–866. Atamt¨ urk, A., Zhang, M., 2007. Two-stage robust network flow and design under demand uncertainty. Operations Research 55 (4), 662–673. Beaujon, G. J., Turnquist, M. A., 1991. A model for fleet sizing and vehicle allocation. Transportation Science 25 (1), 19–45. Ben-Tal, A., Chung, B. D., Mandala, S. R., Yao, T., 2011. Robust optimization for emergency logistics planning: Risk mitigation in humanitarian re40

lief supply chains. Transportation Research Part B: Methodological 45 (8), 1177–1189. Ben-Tal, A., Goryashko, A., Guslitzer, E., Nemirovski, A., 2004. Adjustable robust solutions of uncertain linear programs. Mathematical Programming 99 (2), 351–376. Ben-Tal, A., Nemirovski, A., 1999. Robust solutions of uncertain linear programs. Operations Research Letters 25 (1), 1–13. Benders, J. F., 1962. Partitioning procedures for solving mixed-variables programming problems. Numerische Mathematik 4, 238–252. Bertsimas, D., Brown, D., Caramanis, C., 2011. Theory and applications of robust optimization. SIAM Review 53 (3), 464–501. Bertsimas, D., Iancu, D. A., Parrilo, P. A., 2010. Optimality of affine policies in multistage robust optimization. Mathematics of Operations Research 35 (2), 363–394. Bertsimas, D., Litvinov, E., Sun, X. A., Jinye, Z., Tongxin, Z., 2013. Adaptive robust optimization for the security constrained unit commitment problem. Power Systems, IEEE Transactions on 28 (1), 52–63. Bertsimas, D., Sim, M., 2003. Robust discrete optimization and network flows. Mathematical Programming 98 (1-3), 49–71. Bertsimas, D., Sim, M., 2004. The price of robustness. Operations Research 52 (1), 35–53. Bertsimas, D., Thiele, A., 2004. A robust optimization approach to supply chain management. In: Bienstock, D., Nemhauser, G. (Eds.), Integer Programming and Combinatorial Optimization. Vol. 3064 of Integer Programming and Combinatorial Optimization. Springer-Verlag, Berlin, pp. 86–100. Bertsimas, D., Thiele, A., 2006. A robust optimization approach to inventory theory. Operations Research 54 (1), 150–168. Bespamyatnikh, S., Bhattacharya, B., Kirkpatrick, D., Segal, M., 2000. Mobile facility location. In: Proceedings of the 4th international workshop on 41

Discrete algorithms and methods for mobile computing and communications. Proceedings of the 4th international workshop on Discrete algorithms and methods for mobile computing and communications. ACM, Boston, Massachusetts, United States, pp. 46–53. Birge, J. R., Louveaux, F. O. V., 1988. A multicut algorithm for two-stage stochastic linear programs. European Journal of Operational Research 34 (3), 384–392. Boloori Arabani, A., Farahani, R. Z., 2012. Facility location dynamics: An overview of classifications and applications. Computers & Industrial Engineering 62 (1), 408–420. Carlsson, D., Flisberg, P., Rnnqvist, M., 2014. Using robust optimization for distribution and inventory planning for a large pulp producer. Computers & Operations Research 44 (0), 214–225. Chen, X., Zhang, Y., 2009. Uncertain linear programs: Extended affinely adjustable robust counterparts. Operations Research 57 (6), 1469–1482. Claßen, G., Koster, A. M. C. A., Schmeink, A., 2013. A robust optimisation model and cutting planes for the planning of energy-efficient wireless networks. Computers & Operations Research 40 (1), 80–90. Duran, M., Grossmann, I., 1986. An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Mathematical Programming 36 (3), 307–339. Durocher, S., Kirkpatrick, D., 2006. The steiner centre of a set of points: stability, eccentricity, and applications to mobile facility location. International Journal of Computational Geometry & Applications 16 (04), 345– 371. Erera, A. L., Morales, J. C., Savelsbergh, M., 2009. Robust optimization for empty repositioning problems. Operations Research 57 (2), 468–483. Fletcher, R., Leyffer, S., 1994. Solving mixed integer nonlinear programs by outer approximation. Mathematical Programming 66 (1-3), 327–349. Gabrel, V., Lacroix, M., Murat, C., Remli, N., 2014a. Robust location transportation problems under uncertain demands. Discrete Applied Mathematics 164 (1), 100–111. 42

Gabrel, V., Murat, C., Thiele, A., 2014b. Recent advances in robust optimization: An overview. European Journal of Operational Research 235 (3), 471–483. Gendreau, M., Laporte, G., Semet, F., 2001. A dynamic model and parallel tabu search heuristic for real-time ambulance relocation. Parallel Computing 27 (12), 1641–1653. Halper, R., Raghavan, S., 2011. The mobile facility routing problem. Transportation Science 45 (3), 413 –434. Hoff, A., Andersson, H., Christiansen, M., Hasle, G., Lkketangen, A., 2010. Industrial aspects and literature survey: Fleet composition and routing. Computers & Operations Research 37 (12), 2041–2061. Horst, R., Tuy, H., 1996. Global optimization: Deterministic approaches. Springer-Verlag, Berlin, Heidelberg, NewYork. Jiang, R., Zhang, M., Li, G., Guan, Y., 2014. Two-stage network constrained robust unit commitment problem. European Journal of Operational Research 234 (3), 751–762. Lei, C., Lin, W.-H., Miao, L., 2014. A multicut l-shaped based algorithm to solve a stochastic programming model for the mobile facility routing and scheduling problem. European Journal of Operational Research 238 (3), 699–710. Li, Z., Tao, F., 2010. On determining optimal fleet size and vehicle transfer policy for a car rental company. Computers & Operations Research 37 (2), 341–350. List, G. F., Wood, B., Nozick, L. K., Turnquist, M. A., Jones, D. A., Kjeldgaard, E. A., Lawton, C. R., 2003. Robust optimization for fleet planning under uncertainty. Transportation Research Part E: Logistics and Transportation Review 39 (3), 209–227. Mudchanatongsuk, S., Ordonez, F., Liu, J., 2007. Robust solutions for network design under transportation cost and demand uncertainty. Journal of the Operational Research Society 59 (5), 652–662.

43

Sayarshad, H. R., Ghoseiri, K., 2009. A simulated annealing approach for the multi-periodic rail-car fleet sizing problem. Computers & Operations Research 36 (6), 1789–1799. Song, D.-P., Earl, C. F., 2008. Optimal empty vehicle repositioning and fleet-sizing for two-depot service systems. European Journal of Operational Research 185 (2), 760–777. Thiele, A., Terry, T., Epelman, M., 2009. Robust linear optimization with recourse. Technical Report. Van Mieghem, A. J., 2003. Capacity management, investment, and hedging: Review and recent developments. Manufacturing & Service Operations Management 5 (4), 269–302. Vansteenwegen, P., Souffriau, W., Oudheusden, D. V., 2011. The orienteering problem: A survey. European Journal of Operational Research 209 (1), 1– 10. Vavasis, S. A., 1991. Nonlinear Optimization: Complexity Issues. Oxford University Press, Inc., New York, USA. You, P.-S., Hsieh, Y.-C., 2014. A study on the vehicle size and transfer policy for car rental problems. Transportation Research Part E: Logistics and Transportation Review 64 (0), 110–121.

44