Network design in scarce data environment using moment-based distributionally robust optimization

Computers and Operations Research 88 (2017) 44–57

Contents lists available at ScienceDirect

Computers and Operations Research journal homepage: www.elsevier.com/locate/cor

Network design in scarce data environment using moment-based distributionally robust optimization Hideaki Nakao, Siqian Shen∗, Zhihao Chen Department of Industrial and Operations Engineering University of Michigan, Ann Arbor, MI, USA

a r t i c l e

i n f o

Article history: Received 20 September 2016 Revised 30 June 2017 Accepted 2 July 2017 Available online 3 July 2017 Keywords: Stochastic/robust network design Distributionally robust optimization Mixed-integer linear programming Polynomial linearization Cutting-plane algorithm

a b s t r a c t We consider a network design problem (NDP) under random demand with unknown distribution for which only a small number of observations are known. We design arc capacities in the first stage and optimize single-commodity network flows after realizing the demand in the second stage. The objective is to minimize the total cost of allocating arc capacities, flowing commodities, and penalty for unmet demand. We formulate a distributionally robust NDP (DR-NDP) by constructing an ambiguity set of the unknown demand distribution based on marginal moment information, to minimize the worst-case total cost over all possible distributions. Approximating polynomials with piecewise-linear functions, we reformulate DR-NDP as a mixed-integer linear program optimized via a cutting-plane algorithm. We test diverse network instances to compare DR-NDP with a stochastic programming approach, a deterministic benchmark model, and a robust NDP formulation. Our results demonstrate adequate robustness of optimal DR-NDP solutions and how they perform under varying demand, modeling parameter, network, and cost settings. The results highlight potential niche uses of DR-NDP in data-scarce contexts. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Network design problems (NDPs) refer to an important class of problems arising in a wide range of applications, including traffic control, energy dispatch, disaster relief, and supply chain management. We focus on single-commodity NDPs with uncertain demand and design arc capacities of a given network to satisfy the demand through ad-hoc flows. In particular, we consider a distributionally robust optimization approach, when a small number of demand observations are known, but the data may not be sufficient for deriving an exact distribution that is necessary for applying traditional stochastic programming approaches. In our problem, a network planner intends to build arcs and decide their capacities in a network before knowing actual realizations of random demand. The cost of building an arc is proportional to the capacity of the arc. Once the arcs are built, demand is realized and is met by flows generated from given supply locations to demand sites. The flow cost on each arc is proportional to the amount of flow on the arc. Additionally, any unmet demand is penalized proportional to the amount of demand shortage. The goal is to minimize the total cost of building arcs, flowing commodities, and unmet demand penalty. The latter two types of cost are ran∗

Corresponding author. E-mail addresses: [email protected] (H. Nakao), [email protected], [email protected] (S. Shen), [email protected] (Z. Chen). http://dx.doi.org/10.1016/j.cor.2017.07.002 0305-0548/© 2017 Elsevier Ltd. All rights reserved.

dom due to the uncertain demand. As the realized demand is not known a priori, the planner has to make decisions based on previously observed realizations. However, historical data may be scarce and may not describe the distribution accurately. To this end, we propose a distributionally robust NDP approach (DR-NDP), to conservatively allocate arc capacities as a precaution against observed data that may not be representative of the true demand distribution. Using empirical mean values and standard deviations, we construct an ambiguity set of possible demand distributions and seek optimal design of arc capacities to minimize the worst-case expected flow cost and demand shortage penalty for any distribution in the set. This approach provides more conservative solutions as compared to standard stochastic optimization methods that assume fully known demand distribution. On the other hand, the marginal moments used by DR-NDP can be easily obtained from the observed data. We show later via computational studies that using the moment information, DR-NDP also yields less conservative results than the ones of robust optimization approach relying on very little information about the uncertainty. Having a more complex formulation, it takes longer to solve DR-NDP than the stochastic NDP model constructed by using the Sample Average Approximation (SAA) method and the robust NDP formulation using a budget uncertainty set. But it takes very short amount of time (often within a few seconds) to implement all three approaches.

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

1.1. Applications The above description reveals the opportunities that this research targets. Via DR-NDP, we aim to find relatively conservative arc design solutions, of which few data points are available to accurately determine the underlying true distribution of demand. The presence of the unmet demand penalty also limits its applications to problems in which demand shortage significantly impacts service quality and system reliability. We present two examples applications below. •

•

Humanitarian relief supply network design: When providing humanitarian relief after disaster occurs, little prior data about uncertain demand for aid is likely to be known due to the rare occurrence of disaster. Using DR-NDP, one can determine the most cost-effective way to supply relief to the affected region, and can conservatively design the supply network to provide high-quality relief operations. Supply network design for new products: Releasing new products or introducing existing products to new regions requires careful planning of the supply network. Historical data for forecasting consumer demand is likely to be uncertain and scarce. Thus, using DR-NDP to design supply networks makes sure that products are launched to satisfy most customers’ needs and to keep demand shortage low.

1.2. Literature review NDPs ubiquitously arise in network planning and operational management, and there is a significant amount of literature on the related theories and applications. Magnanti and Wong (1984) describe a unifying framework for deterministic single and multicommodity NDPs, and specify the design of the framework for facility location and traffic network design problems. Specific NDPs are also studied in depth for different practical contexts, e.g., telecommunications and computer network design in Pióro and Medhi (2004) and Babonneau et al. (2013). Poss (2011) provides a detailed summary of models and algorithms for NDPs with diverse network and problem settings. Recently, parameter uncertainty becomes a major concern in NDP related studies, since the practical design of networks is usually in anticipation of future uncertain needs and circumstances. The literature of stochastic NDPs typically assumes fully known distributional information on the demand uncertainty, allowing to solve the problem by using largescale stochastic programs with a finite, but large number of realizations (i.e., “scenarios”) (see, e.g., Santoso et al., 2005; Ukkusuri and Patil, 2009). The scenarios can be generated either from certain assumed distributions (e.g., using the Sample Average Approximation (SAA) method (see Kleywegt et al., 2002)) or based on statistical information (e.g., using a moment-matching method (Høyland et al., 2003)). In cases where the demand distribution is not fully known, another avenue for optimizing NDPs is through robust optimization (Ben-Tal and Nemirovski, 1998; Bertsimas and Sim, 20 03,20 04), to optimize the worst-case objective value for any uncertain parameter realization over a pre-defined uncertainty set. They are well-suited to conservatively solve problems without exact distributional information, but only having small amounts of data to derive beliefs of the uncertainty structure. Atamtürk and Zhang (2007) propose a two-stage robust optimization approach for solving network flow and design problems with uncertain demand. Like our paper, they model the network flow decisions as recourse variables decided after realizing the demand. This allows one to control conservatism of the solutions. They also apply a budget uncertainty set (see Bertsimas and Sim (2004)) under which they provide an upper bound on the probability of a robust so-

45

lution being infeasible for any random demand vector. Altin et al. (2011), and Koster et al. (2013) describe a variety of formulations, complexities, valid inequalities, and computational results for robust NDP with uncertain (hose) demand. Álvarez-Miranda et al. (2014) discuss the complexity and propose heuristic algorithms for robust NDPs; Cacchiani et al. (2016) optimize single-commodity robust NDPs by using the branch-and-cut algorithm based on efficient separation procedures. They consider uncertainty sets of the hose demand modeled as a finite set of scenarios or as a polytope. Robust NDPs are often reformulated as two-stage mixed-integer linear programs, and can be optimized through decomposition, cutting-plane, and/or column-generation methods (see, e.g., Ayoub and Poss, 2016; Lee et al., 2013; Shen et al., 2016; Atamtürk and Zhang, 2007). Furthermore, in (Ben-Ameur and Kerivin, 2005; Mattia, 2013; Poss and Raack, 2013; Scutellà, 2009), the authors investigate robust NDPs with dynamic routing decisions made in multiple stages when sequentially realizing the demand, rather than using static network flows in the two-stage setting. Distributionally robust (DR) optimization approaches recently receive much more attention for optimization under distributional uncertainty. They plan against the worst-case distribution in an ambiguity set that is typically built based on statistical information from historical data. Bertsimas et al. (2010); Delage and Ye (2010); Jiang and Guan (2016) propose different moment-based ambiguity sets for deriving tractable reformulations of DR stochastic or chance-constrained programs. In this paper, we consider a moment-based ambiguity set that is constructed based on the empirical mean values and standard deviations of demand. To the best of our knowledge, this approach is novel for designing network arc capacities with ad-hoc flows. Later, we will also numerically demonstrate that DR-NDP yields adequately robust but not over conservative solutions that are less sensitive to changes of distribution, demand, network, and cost parameters, compared with stochastic NDP or robust NDP. 1.3. Structure of paper The remainder of the paper is organized as follows. In Section 2, we describe a generic NDP and the formulation of a marginal moment-based ambiguity set for DR-NDP. In Section 3, we reformulate DR-NDP as a single-level optimization model with mixed-integer linear programming constraints, and we approximate the quadratic objective by using a piecewise-linear function. In Section 4, we develop a cutting-plane algorithm for optimizing the reformulation. In Section 5, we demonstrate the computational results of DR-NDP and conduct sensitivity analysis by testing diverse network instances and parameters. In Section 6, we conclude the paper and discuss future research directions. 2. Problem formulation Consider a directed network G(N, A) with N being the set of nodes, and A being the set of arcs. The goal is to minimize the total cost of arc capacity allocation, network flows, and unmet demand penalty. Without loss of generality, the network has disjoint sets of supply nodes S ⊂ N and demand nodes D ⊂ N. We assume that each supply node i ∈ S has deterministic supply si > 0 and each demand node i ∈ D has uncertain demand di ≥ 0. The demands are assumed to be jointly distributed and, while the true distribution of the demand vector d = (di , i ∈ D )T is not known, we are given a set of |K| data observations of d, denoted by {dk }k ∈ K . We assume that d has a boxed-shaped support , given by  := {d ∈ |D| |D| R+ : di ≤ di ≤ di , ∀i ∈ D}, where d, d ∈ R+ represent the vectors of lower and upper bounds of d, respectively. These bounds can be obtained from historical demand data or via prior knowledge in specific problem contexts. (For example, the number of users at a

46

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

node of an internet service network is upper bounded by the local population.) The design of a box-shaped support is analogous to the use of a box-shaped uncertainty set in robust optimization (see Soyster, 1973). 2.1. Formulation of the NDP under demand uncertainty

2.2. Formulation of DR-NDP with a moment-based ambiguity set

We consider the problem in two stages. First, before knowing the actual demand, we decide the capacity xi j ∈ R+ of each arc (i, j) ∈ A, each of which incurs a cost cij > 0. Since any arc capacity should not exceed the maximum total demand, we construct the feasible region of the decision vector x = (xi j , (i, j ) ∈ A )T as



|A|

x ∈ R+ : 0 ≤ xi j ≤

X :=



di ,

∀(i, j ) ∈ A .

min y,t

ai j yi j +



(i, j )∈A

x∈X

(i, j )∈A



∀(i, j ) ∈ A y ji ≤ si



yi j −

(1c)

y ji = 0

∀i ∈ N \ ( S ∪ D )

(1d)



(i, j )∈A

∀i ∈ D

( j,i )∈A

y, t ≥ 0,

(1f) |A|

where y ∈ R+ denotes the vector of flows yij on each arc (i, j) ∈ |D|

A and t ∈ R+ denotes the vector of unmet demand ti at each demand node i ∈ D. Each unit of flow on arc (i, j) incurs a flow cost aij > 0, and each unit of unmet demand at demand node i ∈ D incurs a penalty cost Gi > 0 in the objective function. Constraint (1b) restricts the flow amount on each arc being no more than the allocated capacity. Constraints (1c)–(1e) are flow balance constraints for each supply, transmission, and demand node, respectively. Constraint (1e) also determines the amount of unmet demand ti at each demand node i ∈ D. If the distribution F of the demand vector d is fully known, then a risk-neutral stochastic programming formulation of NDP is given by



min x∈X





ci j xi j + EF [g(x, d )] ,

(2)

(i, j )∈A

where EF is the expectation taken under the true distribution function F of d. If data is scarce, the observed realizations will not

.

(3)

1 

|K |

∀i ∈ D

dik

(4a)

k∈K

2 1  k d i − μi |K |

∀i ∈ D.

(4b)

 M = M0 (, μ, σ ) :=





(1e)

F ∈M

The marginal moment-based ambiguity set is formulated as



y ji − ti ≤ −di

 

k∈K

(1b)

( j,i )∈A

yi j −

σi2 :=



∀i ∈ S

ci j xi j + max EF g(x, d )

(i, j )∈A

μi :=

(1a)

( j,i )∈A

(i, j )∈A





yi j −

Gi ti

 



We construct the ambiguity set using the marginal mean vec|D| |D| tor μ ∈ R+ and variance vector σ 2 ∈ R+ of the random demand, respectively. They consist of entries

i∈D

s.t. yi j ≤ xi j 



min

In specific NDP contexts, X can include other constraints. For example, budget constraints may be considered for supply chain NDP to limit the total capacity allocation cost; for transportation NDPs, it may be desirable for certain nodes to be less congested (e.g., road junctions near schools and hospitals), which can be achieved by adding constraints that require the total capacity of arcs leading out of a node no less than the total capacity leading into the node. To keep our model generic, we only consider X in the aforementioned form in this paper. Given a first-stage decision x and demand d, the second-stage cost g(x, d) is a capacitated minimum cost flow problem that penalizes possible unmet demand, formulated as

g(x, d ) :=

We first establish an ambiguity set M of candidate distributions, from which a distribution F of random demand d is found to maximize the expected cost EF [g(x, d )] given first-stage decision x and demand d. We formulate the DR-NDP variant of model (2) as:



i∈D



be sufficient to determine F with reasonable accuracy. Next, we employ distributionally robust optimization to tackle NDPs with scarce data and use the marginal moments to construct the ambiguity set of possible F for DR-NDP.

F:

d F (d ) = 1

d i d F ( d ) = μi

(5a)

∀i ∈ D

di2 dF (d ) = μ2i + σi2

(5b)

∀i ∈ D}.

(5c)

Constraint (5a) ensures that M0 only contains valid distributions over the support set  of d, while constraints (5b) and (5c) restrict these distributions to have marginal means and variances equal to those of the observed demand, respectively. The solution time of DR-NDP using this set is independent in the number of observations – since the ambiguity set only uses μ and σ 2 , which are pre-calculated via equations (4a) and (4b), respectively. Remark 1. The ambiguity set can be constructed through other means, such as using nested confidence sets constructed with the decision maker’s prior information (Wiesemann et al., 2014), using φ -divergences from the observed density (Ben-Tal et al., 2013), or using an assumption of radial distributions on the data (Calafiore and El Ghaoui, 2006). Notably similar to the marginal momentbased ambiguity set is the cross moment-based ambiguity set, which takes into account the correlation between components of the demand vector d. We choose to use marginal moments as they are more suited for our problem context. When data is scarce, it is difficult to meaningfully obtain accurate correlation information between demand at different nodes. Furthermore, there is a lack of any cause-effect relation information between nodes in a generic NDP. As a precaution against “bad” data, we only consider marginal moments and do not impose any constraints on the cross moments, thereby computing solutions that are robust over a larger set of possible distributions.

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

3. Reformulation of DR-NDP

is equal to that of



We develop methodologies for reformulating and eventually solving the DR-NDP model (3) using the ambiguity set characterized in (5a)–(5c). We begin by dualizing the second-stage maximization problem

 

max EF g(x, d )



(6)

F ∈M

to obtain a minimization problem that can be minimized together with the first-stage arc-capacity design problem. This section details the dualization procedures and a subsequent approximation of the objective function of DR-NDP. 3.1. Dualizing the second-stage problem When M = M0 (, μ, σ ), we explicitly consider the secondstage problem as



max





F

g(x, d ) dF (d ) : (5a ) − (5c ) .

(7)

Both x and d are input parameters of model (7) and the model contains the continuous decision variable F and linear constraints (5a)–(5c). Therefore, we take the dual of (7) as a linear program with primal variable F. We associate dual variables z ∈ R, r ∈ R|D| and v ∈ R|D| with constraints (5a)–(5c), respectively. The dual is a semi-infinite linear program

min z,r,v

z+



μi r i +

i∈D

s.t. z +





μ2i + σi2 vi

di ri +

i∈D

∀ d ∈ .

min

x,z,r,v

z+



μi r i +

i∈D





ci j xi j : x ∈ X, (9 ) . μ2i + σi2 vi + (i, j )∈A

i∈D

(10) However, this formulation remains intractable as constraint (9) is a semi-infinite constraint (i.e., a constraint that must be fulfilled over the infinite set ). 3.2. Reformulating the semi-infinite constraint



z ≥ max g(x, d ) − d∈



di ri −

i∈D



v .

(11)

To further reformulate the right-hand side of (11), we formulate the dual problem of the minimization model g(x, d) to merge with the maximization over d ∈ . Lemma 1. For fixed r, v and x feasible to problem (10), the value of



max g(x, d ) − d∈

 i∈D

di ri −

 i∈D



di2 vi

di2 vi

i∈D

(14a) (14b)

πi ≥ 0

∀i ∈ S

(14c)

0 ≤ πi ≤ G i

∀i ∈ D}.

(14d)

Proof. We begin by taking the dual of g(x, d). We associate the |A| dual variables θ ∈ R+ , π i for i ∈ S, π i for i ∈ N\(S ∪ D), and π i for i ∈ D with constraints (1b)–(1e) in the problem g(x, d), respectively. We then consider g(x, d) in its dual form as:





g(x, d ) = max θ ,π

d i πi −



i∈D

s i πi −





xi j θi j : (14a ) − (14b ) .

(i, j )∈A

i∈S

(15) Substituting the dual model (15) into (12) yields







max

(θ ,π )∈P



di ri −

i∈D

d i πi −

i∈D





s i πi −



x i j θi j

(i, j )∈A

i∈S





di2 i

v .

(16)

i∈D

We merge the maximization over (θ , π ) ∈ P with the maximization over d ∈  to complete the reformulation and the proof.  Problem (13) still cannot be solved directly by off-the-shelf optimization solvers, as it contains bilinear terms di π i in its objective, with both di and π i being continuous for all i ∈ D. Since π i is bounded below and above, we rewrite di in terms of binary variables to reformulate the bilinear terms. Consequently, the following lemma further reformulates problem (13) by leveraging the fact that d has a box-shaped support. Lemma 2. For fixed r, v and x feasible to problem (10), the value of problem (13) is equal to that of



−



i∈D



s i πi −

 max di ≤di ≤di



d i πi − d i r i −

di2 i

v





x i j θi j .

(17)

(i, j )∈A

i∈S

Proof. Given the separable feasible region  of d, we separately maximize the objective function (θ , π ) ∈ P and then over d ∈ :

max

(θ ,π )∈P

−

 i∈S

(12)

di ri −



∀(i, j ) ∈ A



i∈D

i∈D



θi j ≥ 0

max



di2 i

x i j θi j −

(i, j )∈A

i∈S

π j − πi − θi j ≤ ai j ∀(i, j ) ∈ A

(θ ,π )∈P

Because constraint (9) must be satisfied for all d ∈ , it must be satisfied for the worst-case d. We move all the terms containing d to the right-hand side of constraint (9) and reformulate it as

s i πi −



P := {(θ , π ) :

−



d i πi −

i∈D



where P is a polyhedron described as

(9)

i∈D



(13)

d∈

Next, we replace (6) with constraints (8)–(9) in the DR-NDP model (3), and merge the minimization objective of the dualized second-stage problem with the minimization objective in the first stage to obtain the reformulation of DR-NDP model (3) under the moment-based ambiguity set M0 (, μ, σ ):



max

max

di2 vi ≥ g(x, d )



d∈,(θ ,π )∈P

(8)

i∈D



47



max d∈

s i πi −



d i πi −

i∈D



 i∈D



x i j θi j .

di ri −





di2 vi

i∈D

(18)

(i, j )∈A

Note that the order of maximization in problem (18) is swapped from that in problem (16), which does not affect the optimal solution given that set  is compact. Moreover, the maximization

48

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

over d ∈  is separable by the indices i ∈ D since the box-shaped support  is defined by independent lower and upper bounds in each dimension. We solve |D| maximization problems given in (17), each over a simple interval, to complete the reformulation and the proof.  Theorem 1. For fixed r, v and x feasible to problem (10), the value of problem (12) is equal to that of problem (17). The proof of Theorem 1 directly follows Lemmas 1 and 2. Through Lemma 2, we first consider each maximization problem



di πi − di ri − di2 vi

max di ≤di ≤di



(19)

for all i ∈ D separately. Although obtaining an analytic solution of (19) in terms of π i , ri , and vi is trivial, the resulting objective function remains bilinear and the cuts generated using this reformulation cannot be optimized directly by solvers. In our approach, we approximate the objective function of (19) as a piecewise-linear function with Ni intervals, and maximize over the piecewise-linear function. Since the maximum over the piecewise-linear function lies on exactly one of the Ni + 1 interval endpoints, the problem simplifies to finding the maximum of Ni + 1 different values. We define intervals of equal length such that

din = di +

n

Ni



di − di ,

∀n = 0, 1, . . . , Ni ,

din πi − din ri −

max

n∈{0,1,...,Ni }

2 din i

v

(20)

ρin = 1

(21a)

n=0

ρin ∈ {0, 1}, ∀n = 0, 1, . . . , Ni ,

(21b)

to ensure that exactly one endpoint is selected to maximize problem (20). Subsequently, problem (20) is equivalent to



max ρi

Ni 



ρin din πi − din ri − din2 vi : (21a ) − (21b ) 

(22)

n=0

for all i ∈ D, where ρi ∈ RNi denotes the vector of ρ in for n = 2v 0, 1, . . . , Ni . Then, the objective takes the value din πi − din ri − din i if ρin = 1 (and the remaining SOS1 variables in the same set are zero), or equivalently, when d maximizes problem (20). in

Since ρ in is binary and (14d) gives us the bounds 0 ≤ π i ≤ Gi , we can subsequently reformulate the bilinear terms ρ in π i using the McCormick inequalities (see McCormick, 1976). The following McCormick inequalities introduce auxiliary variables λin with λin ≡ ρ in π i for all n = 0, 1, . . . , Ni , i ∈ D:

λin − πi − Gi ρin ≥ −Gi

λin ≥ 0∀n = 0, 1, . . . , Ni , i ∈ D.

(23d)

When ρin = 1, constraints (23a) and (23b) bound λin below and above by π i , and when ρin = 0, constraints (23c) and (23d) bound λin below and above by 0. This results in the desired equivalence relation, namely, λin ≡ ρ in π i . Using the approximation of (19) with (22) and the above linear McCormick constraints (23a)–(23d), we approximate the optimal objective value of model (12), and consequently model (17), by a mixed-integer linear program:

h(x, r, v ) :=

s.t.

f (x, r, v, θ , π , ρ , λ )

max

θ ,π ,ρ ,λ

(24)

(14a ) − (14d ), (23a ) − (23d ), (21a ) − (21b )

where Ni 

f (x, r, v, θ , π , ρ , λ ) : =



2 din λin − din ri ρin + din vi ρin

i∈D n=0

−





s i πi −

x i j θi j .

∀n = 0, 1, . . . , Ni , i ∈ D

(23a)

∀n = 0, 1, . . . , Ni , i ∈ D

(23b)



(25)

(i, j )∈A

Therefore, we derive an approximation of the DR-NDP model (10) as

APPROX :

min

x,z,r,v

z+



μi r i +

i∈D





ci j xi j μ2i + σi2 vi +

i∈D

(i, j )∈A

(26a) s.t. z ≥ h(x, r, v )

(26b)

x ∈ X.

(26c)

APPROX, while having only linear objective and constraints with binary and continuous variables, is not a mixed-integer linear program due to constraint (26b) having a maximization problem on the right-hand side. Treating the embedded maximization problem h(x, r, v) as a separation problem, we develop a cuttingplane algorithm for solving DR-NDP, detailed in Section 4. Approximating (17) with h(x, r, v), we replace |D| continuous variables with i ∈ D Ni binary variables and i ∈ D Ni continuous  variables, and also introduce 4 i∈D Ni + |D| new constraints to the formulation. Although the binary variables comprise |D| sets of SOS1 variables, the introduction of a large number of binary variables could reduce computational efficiency. While larger Ni for each i ∈ D yields a better approximation, the accuracy of the approximation will come at a cost of longer computational time. Later in our computational studies, we demonstrate how different choices of Ni affect the optimal objective values and solution time. 4. A cutting-plane algorithm We describe a cutting-plane algorithm that iteratively passes solutions from a master problem to the maximization problem h(x, r, v) and generates a linear valid inequality with variables z, r, v, and x to the master problem, when the current master problem’s solution does not satisfy the relationship in (26b). We drop constraint (26b) from APPROX and formulate the relaxed master problem at the Mth iteration as

MASTER : min

λin − πi ≤ 0

(23c)



We define binary variables ρ in ∈ {0, 1}, for each n = 0, 1, . . . , Ni , such that ρin = 1 if the nth endpoint din yields the maximum for problem (20), and ρin = 0 otherwise, for each i ∈ D. Each set {ρin }n=0,1,...,Ni , for all i ∈ D, is a Specially Ordered Set of Type 1 (SOS1), containing binary variables that sum to one, such that for all i ∈ D, Ni 

∀n = 0, 1, . . . , Ni , i ∈ D

i∈S

is the nth endpoint for each i ∈ D. Given these endpoints, we approximate problem (19) as the largest value obtained when substituting di with one of the endpoints in {din }n=0,1,...,Ni , or equivalently,



λin − Gi ρin ≤ 0

x,z,r,v

z+

 i∈D

μi r i +





ci j xi j (27) μ2i + σi2 vi +

i∈D

(i, j )∈A

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

s.t. z +

 i∈D

μi r i +





ci j xi j ≥ 0 μ + σ vi + 2 i

i∈D

z ≥ f (x, r, v, θ m , π m , ρ m , λm )

2 i

(28)

(i, j )∈A

∀m = 1 , . . . , M

x ∈ X,

(29) (30)

where constraint (28) enforces a trivial lower bound and (29) are the cuts generated during iterations 1, . . . , M. The solution of the subproblem in the mth iteration is denoted by (θ m , π m , ρ m , λm ), for m = 1, . . . , M. The linear function f(x, r, v, θ m , π m , ρ m , λm ) in (29) follows the form of (25) with (θ , π , ρ , λ ) = (θ m , π m , ρ m , λm ). These cuts are iteratively generated by passing the current solution (zM , xM , rM , vM ) of MASTER to the subproblem h(x, r, v) and by checking whether the solution satisfies (26b). Should the current solution be infeasible to (26b), we add a new cut (29) in iteration M + 1 to MASTER. We present the algorithmic details in Algorithm 1. Algorithm 1 A cutting-plane algorithm for solving DR-NDP. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

12: 13: 14: 15: 16: 17:

for i ∈ D do Select Ni piecewise intervals. end for Select a non-negative tolerance value . Initialize iteration counter M = 0. Initialize continue = true. while continue do M =M+1 continue = false Solve MASTER and denote optimal solution by ( xM , z M , r M , vM ). Solve h(xM , r M , vM ) via branch-and-bound with each node solved via the Simplex algorithm, and denote optimal solution by (θ M , π M , ρ M , λM ), respectively. if zM < (1 − ) f (xM , r M , vM , θ M , π M , ρ M , λM ) then continue = true Add z ≥ f (x, r, v, θ M , π M , ρ M , λM ) to MASTER as the Mth cut. end if end while STOP; (xM , zM , r M , vM ) is optimal to APPROX.

Remark 2. Algorithm 1 converges to an optimal solution to APPROX in a finite number of iterations. First, note that h(x, r, v) has a feasible region independent of r, v and x. Furthermore, the feasible region at each node of the branch-and-bound algorithm for h(x, r, v) is polyhedral. If the problem at each node is solved via the Simplex algorithm, its optimal solution is associated with an extreme point. Consequently, each cut in (29) corresponds to some extreme point of the feasible region of some branch-and-bound node. After solving MASTER, we obtain a tentative optimal solution that satisfies all previously generated cuts, and a new cut is generated only if it is violated by this solution, so each cut (29) corresponds to a unique extreme point. Furthermore, there are a finite number of binary variables ρ in to branch on and the polyhedral feasible region of each branch-and-bound node must have a finite number of extreme points. Hence, Algorithm 1 will terminate after finite iterations.

robust NDP based on a budget uncertainty set. We test all the approaches using the same set of randomly generated network instances and demand data. We compare their optimal solutions under different cost settings, network topologies, demand-supply locations, and characteristics of observed demand data. We use the same reference set of demand samples to evaluate the solutions given by different approaches and demonstrate their out-of-sample performance. We will show that the DR-NDP model yields adequately robust solutions for varying demand and penalty cost. It yields smaller and less variable gap values between the estimated in-sample objective cost and the actual cost in the out-of-sample test for all the instances we test. In Section 5.1, we describe the detailed formulations and solution approaches for stochastic NDP and robust NDP. In Section 5.2, we describe our experimental design of network instances and test parameters. In Section 5.3, we compare the CPU time of solving DR-NDP, stochastic NDP, deterministic benchmark, and robust NDP, and specifically report the time and number of cuts for each replication of a DR-NDP instance. In Section 5.4, we demonstrate both in-sample and out-of-sample performance of optimal solution given by each approach. In Section 5.5, we conduct sensitivity analysis by varying the unit penalty cost, the number of intervals for constructing the piecewise-linear approximation of the quadratic objective function of DR-NDP, and the number of demand nodes. We implement all the optimization models and solution algorithms in Java and solve them using Gurobi 6.5.2. The results below are obtained using an Intel Core i7-3770 CPU @ 3.4GHz and 16GB RAM. 5.1. Benchmark approaches 5.1.1. SAA-based stochastic NDP model We employ the SAA method (Kleywegt et al., 2002) to formulate the stochastic NDP model, for which we generate i.i.d. demand samples from a pre-assumed distribution. In the stochastic NDP, the goal is to minimize the capacity design cost plus the expected total cost of network flows and unmet demand penalty. For all the benchmark approaches, we consider the same initial set of data samples with marginal mean values μi and variances σi2 for all i ∈ D. We use some prior distribution type to sample the demand realizations needed by the SAA model. We first consider Uniformly distributed demand at each demand node and sample each demand realization diω , ∀i ∈ D from a Uniform distribution with √ √ lower bound μi − 3σi and upper bound μi + 3σi , to match the empirical means and variances of the initial data set. Alternatively, we consider Gamma distributed demand and sample each realization diω from a Gamma distribution with shape μ2i /σi2 and scale σi2 /μi , ∀i ∈ D. Letting the set contain all the scenarios used by the SAA method, we denote the finite number of demand realizations by {dω }ω ∈ . Following the SAA approach, we assume that the realizations are sampled with equal probability 1/| | and formulate a scenariobased two-stage stochastic program as follows:



min x∈X

 (i, j )∈A

ci j xi j +

1 

| | ω ∈



g(x, dω ) .

(31)

As | | → ∞, the optimal objective value of model (31) converges to the one of the stochastic NDP:



5. Computational results

min We compare the following approaches, including (i) DR-NDP, (ii) stochastic NDP with pre-assumed demand distribution, (iii) a deterministic benchmark model using mean demand values, and (iv)

49

x∈X



ci j xi j + EF0





g(x, dω )



,

(i, j )∈A

where EF0 is the expectation under an assumed distribution F0 of d. (For example, in this paper, we use the above Uniform or

50

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

Gamma distribution as F0 .) This benchmark approach aims to provide risk-neutral network design solutions under the assumed distribution. The L-shaped method (or equivalently the Benders decomposition approach) (vanSlyke and Wets, 1969) is often employed to derive cuts and iteratively optimize (31), which could be large scale depending on the size of the sample set (see, e.g., Botton et al., 2013; Fortz and Poss, 2009). The Benders cuts underestimate the second-stage value function in the first-stage decision variables by using affine functions and thus can be considered using a simple affine rule to obtain/approximate the optimal objective value of stochastic NDP.

We let ui ≡ qi π i and use the McCormick inequalities (McCormick, 1976) to obtain an equivalent mixed-integer linear programming formulation for the inner problem:

5.1.2. Robust NDP using a budget uncertainty set We describe another benchmark approach using the two-stage robust NDP model introduced by Atamtürk and Zhang (2007), which aims to design a network to minimize the worst-case (i.e., maximum) second-stage cost of network flows and demand shortage over all possible demand. We consider a parametric budget uncertainty set of the demand vector d, first introduced by Bertsimas and Sim (2004), and specify it as

u i ≥ πi − ( 1 − q i ) G i

 ( ) : = d ∈ R|+D| : di ≤ di ≤ di , ∀i ∈ D ⊆ D, |D | ≤ ,  d i = μi , ∀ i ∈ D \ D .

 ( ) : = d ∈ R|+D| : μi − (μi − di )qi ≤ di ≤ μi + (di − μi )qi , ∀i ∈ D,



 qi ≤ .

(32)

i∈D

We formulate the first-stage problem of the robust NDP as



min x∈X





ci j xi j + max {g(x, d )} ,

(33)

d∈( )

(i, j )∈A

π ,θ ,q

d∈( )

−

 i∈S

s.t.



( μ i + ( d i − μ i ) q i ) πi

i∈D

s i πi −



x i j θi j

(34a)

(i, j )∈A

(14a ) − (14d ) 

qi ≤

(34b)

i∈D

qi ∈ {0, 1}

∀i ∈ D,

p(x, θ , π , u )

(35a)

(14a ) − (14d ), (34b ) − (34c )

s.t.

∀i ∈ D

ui ≤ Gi qi u i ≤ πi

(35b)

∀i ∈ D

(35c)

∀i ∈ D

(35d)

∀i ∈ D,

ui ≥ 0 where

p(x, θ , π , u ) =



μ i πi +

i∈D



( d i − μi ) u i −

i∈D

(35e)  i∈S

s i πi −



x i j θi j ,

(i, j )∈A

and Gi corresponds to the unit cost for penalizing unmet demand, defined in (1a), which is the upper bound of the dual variable π i , ∀i ∈ D. We optimize the robust NDP model (33) using the Benders decomposition method (see Birge and Louveaux, 2011; vanSlyke and Wets, 1969), for which at iteration M, we solve a relaxed master problem



min

x∈X, η≥0





ci j xi j + η : η ≥ p(x, θ , π , u ), ∀m = 1, . . . , M , (36) m

m

m

(i, j )∈A

with optimality cuts η ≥ p(x, θ m , π m , um ), m = 1, . . . , M, each generated by solving the subproblem (35a)–(35e) using the optimal solution xm to the relaxed master problem at iteration m. We

continue the iterative process until at some iteration M , ηM ≥



M M M M p(x , θ , π , u ) is satisfied. Again, note that the cuts are linear given the linear function p in (35f). Therefore, the above procedures also use an affine rule to iteratively estimate the optimal objective value of the robust NDP model. Since the stochastic NDP and the robust NDP have different second-stage objective functions, the forms of the corresponding Benders cuts are also different although they are both affine. We will compare the solution method for DR-NDP with the affine rule used by the robust NDP later. 5.2. Experimental setup

where ( ) is in the form of (32) and the second-stage problem g(x, d) is given by (1a)–(1f). Recalling the dual of g(x, d) derived in (15), we reformulate the inner maximization problem in (33) as:

max {g(x, d )} : = max

π ,θ ,q,u

(35f)

Parameter ∈ Z+ indicates the total number of demand nodes allowed to have a demand realization that deviates from its mean value. We use the same notation in DR-NDP for denoting parameters d i , di , and μi , which represent the upper bound, lower bound, and empirical mean value of the random demand di , respectively, for all i ∈ D. Using the techniques in Bertsimas and Sim (2004), we define binary variables qi for each i ∈ D, such that qi = 1 if the realization of demand di deviates from its mean value, and qi = 0 otherwise. The above uncertainty set ( ) can be reformulated by using linear constraints as follows:

qi ∈ {0, 1},

max {g(x, d )} := max

d∈( )

(34c)

where the first-stage network design decision x is an input parameter, and bilinear products qi π i , i ∈ D appear in the objective.

5.2.1. Network settings We mainly test the NOBEL-US network from the SNDlib (Orlowski et al., 2010), which is a library of test instances for Survivable fixed telecommunication Network Design1 . It represents a real-world network based on a fiber-optic backbone telecommunication network in the United States. For some tests, we also consider the NOBEL-EU network and the PIORO40 network from SNDlib, and a 6x6 grid network. Table 1 summarizes the size of each network given by the number of nodes and arcs. For the default case, we randomly locate three supply nodes and three demand nodes in each network and make sure that the supply nodes and demand nodes do not overlap. We do not differentiate the supply nodes and make each supply node i ∈ S have si = 20 units of capacity. We generate 10,0 0 0 i.i.d. samples of the demand vector d by sampling realizations of di at each demand 1

The detailed network structures are presented at http://sndlib.zib.de

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57 Table 1 Sizes of test networks. Network

# of nodes

# of arcs

nobel-us nobel-eu pioro40 grid 6x6

14 28 40 36

42 82 178 120

Table 2 Parameters used in each NDP approach variant. NDP Variant

DR-NDP

Stochastic

Deterministic

Robust

Parameter

d i , di , μi , σ i , Ni

| | samples from a distribution

di = μi

d i , di , μi ,

node i ∈ D following a Gamma distribution with shape 4 and scale 5. However, if any sampled demand value is less than 0 or larger than 50, we discard it and re-sample. Thus, the mean value and standard deviation of demand at each demand node are roughly 20 and 10, respectively, and the total supply capacity approximately equals to the expected value of the total demand. These 10,0 0 0 samples are considered to form a reference set to present all the demand realizations in reality. (Therefore, the true demand distribution is Gamma.) Later we select 20 demand samples from the reference set to compute their empirical demand mean values and standard deviations used by DR-NDP. Assuming certain distribution type (not necessarily Gamma), we use the same empirical means and standard deviations to generate data in the test sample set , used in the SAA model (31) for stochastic NDP. The deterministic benchmark model directly uses the empirical mean to replace the random demand in NDP model. The robust NDP uses the same demand bounds as the ones used by DR-NDP but not the empirical moment information. We design the cost parameters used by DR-NDP, stochastic NDP, deterministic benchmark, and robust NDP as follows. For each test instance, we randomly sample the cost cij of allocating one unit capacity to arc (i, j) from a Normal distribution N (40, 36 ) with 40 being the mean and 6 being the standard deviation, and sample the unit flow cost aij on each arc (i, j) from N (10, 4 ). For both costs, we discard any negative cost value and re-sample. These two cost parameters are randomly sampled to avoid bias from certain cost ratios, while we keep cij significantly larger than aij on average and with larger variance. We let the unit penalty cost Gi of unmet demand at each demand node i ∈ D be 10(|N| − 1 ), where 10 is the mean value of the unit flow cost and |N| is the number of nodes in a given network. In this way, the unit penalty cost in the default case increases as the network size increases, to compensate potential increase of total flow cost in larger networks. The cost values used in our computational studies are unit-free as the relationship between the costs is more relevant to this study than the absolute cost values. Later in our sensitivity analysis, we vary the unit penalty Gi , ∀i ∈ D via multiplying the default penalty by a certain integer and examine how Gi affects the optimal solutions and their performance of different approaches. 5.2.2. Model settings and design of experiments In addition to the shared parameters, in Table 2 we recall the parameters used by each NDP variant. For both DR-NDP and robust NDP, we use di = 0 and di = 50 for all i ∈ D. For DR-NDP, we use Ni = 25 intervals for each i ∈ D in the default setting, yielding cut points 0, 2, . . . , 50 of each demand interval [di , di ] for the piecewise-linear approximation of the quadratic objective in model (19). We also use Ni = 5 and 10 to vary the interval size and examine the tradeoff between the time for solving DR-NDP and the solution quality obtained by using the piecewise-linear approximation. For the robust NDP, we

51

use = 1, 2 in (34b) representing cases where one and two demand nodes can deviate their realized demand values from given mean values (i.e., 20 as the true demand mean), respectively. For the SAA model (31) for stochastic NDP, we use the sample size | | = 50 and = 10 0 0. We follow a distribution type (either Uniform or Gamma) and use the empirical mean values and standard deviations of 20 demand observations described below to generate i.i.d. demand samples in the set . For each reference set with 10,0 0 0 demand data, we draw three sets of 20 data points, representing three different cases of scarce demand data given to a decision maker. The first regular demand set is drawn such that each marginal mean of the observations in the set deviates no more than 7.5% from that of the reference set. The second very low demand set comprises the 20 realizations in the reference set with the lowest total demand. The third low demand set is drawn such that difference between the total marginal means of realizations in the set and that of the very low demand set is at most 70% of the difference between the total true marginal means and that of the very low demand set. In all the tables, we name the DR-NDP approach “drndp”, the stochastic NDP approach using Gamma or Uniform distribution and the test sample size being 50 or 10 0 0 scenarios as “gam50”, “gam10 0 0”, “uni50”, and “uni10 0 0”, respectively. We name the deterministic benchmark model using mean demand values “det” and the robust NDP approach with = 1, 2 “rndp1” and “rndp2”, respectively. The four network names are presented in lower case. 5.3. Results of CPU time Table 3 uses instances of each network type with the default parameter setting and a test sample of very low demand data to compare the average CPU time of running all the replications by using different approaches. We also demonstrate the time fluctuation of running all the ten replications for each instance, by presenting the maximum offset (in percentage) of either the maximum or the minimum solution time that deviates from the average time. It is evident that although the DR-NDP approach takes longer time than all the other approaches, all the approaches are very efficient to implement, including DR-NDP, which uses no more than 5 seconds for computing any instance. In general, the computational time of each approach increases as we increase the network size and thus PIORO40 instances take the longest time to solve. We note that the DR-NDP also has smaller variance in solving all the replications as compared to other approaches, reflected by the maximum CPU time offset values reported for all the networks. This is because DR-NDP does not depend on samples of a certain distribution, which could vary in different replications, and it also does not depend on only one or a few parameters like in the deterministic model and in robust NDP, respectively. The above CPU time trends are similar in all the test instances and for brevity, we omit the detailed results for all the other parameter settings and demand cases. We continue to investigate the details of implementing DR-NDP, using instances of the NOBEL-US network with default parameter settings. In Table 4, we report the results for both regular and very low demand samples and show the total computational time (in Column “time (s)”), the total time of solving the second-stage separation problem (in Column “2ndT (s)”), their ratios (in Column “ratio (%)”), and the total number of cuts (29) generated, which is the same as the number of iterations (in Column “#cuts”). Clearly, most of the computation time is devoted to solving the separation problem (24), which is a mixed-integer linear program taking significant amount of time to optimize. The number of cuts is similar among all the instances, which indicates that the deviation in the solution time for optimizing the separation problem is mainly due to the branch-and-bound procedure.

52

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57 Table 3 CPU time of different approaches for all network instances with very low demand. Approach

drndp gam50 gam10 0 0 uni50 uni10 0 0 det rndp1 rndp2

Average CPU time (s)

CPU time maximum offset (%)

nobel-us

nobel-eu

pioro40

grid

nobel-us

nobel-eu

pioro40

grid

1.242 0.005 0.005 0.006 0.005 0.005 0.025 0.036

2.225 0.013 0.011 0.013 0.012 0.015 0.111 0.122

3.328 0.037 0.035 0.037 0.035 0.036 0.487 0.412

2.585 0.020 0.019 0.021 0.019 0.029 0.229 0.304

28.30 50.94 88.68 106.90 50.94 66.67 25.98 17.32

31.22 68.99 68.14 68.00 66.10 73.68 69.52 72.95

57.09 87.17 87.50 83.29 84.44 130.56 62.46 75.80

54.77 133.83 143.39 128.16 118.75 135.29 114.41 72.07

Table 4 DR-NDP solution time and iteration details for default NOBEL-US instance. Inst.

1 2 3 4 5 6 7 8 9 10

very low demand sample

regular demand sample

time (s)

2ndT (s)

ratio (%)

#cuts

time (s)

2ndT (s)

ratio (%)

#cuts

1.286 1.098 1.216 1.332 1.279 1.188 1.593 1.356 1.056 1.012

1.248 1.062 1.186 1.302 1.257 1.153 1.566 1.327 1.034 0.992

97.05 96.72 97.53 97.75 98.28 97.05 98.31 97.86 97.92 98.02

75 73 73 69 63 76 63 80 71 60

3.729 5.314 29.018 3.417 10.132 4.535 4.743 5.594 3.924 4.992

3.663 5.242 28.958 3.385 10.093 4.504 4.699 5.544 3.882 4.948

98.23 98.65 99.79 99.06 99.62 99.32 99.07 99.11 98.93 99.12

74 73 86 69 72 73 85 90 80 81

5.4. Comparison of in-sample and out-of-sample performance We implement all the approaches using each of the 20-data sample set to obtain an optimal capacity design solution xˆ. We also record the in-sample performance of each optimal solution, namely, the optimal objective costs in the first and second stages of each model. Then, we optimize the capacitated minimum cost flow problem (1a)–(1f) with fixed capacity x = xˆ, for each of the 10,0 0 0 demand vectors in the reference set. The results of the 10,0 0 0 trials are averaged to obtain the out-of-sample expected costs (e.g., expected flow cost and expected penalty cost of unmet demand). We test 10 replications of each instance, each of which uses identical parameter settings, except that we randomly select supply and demand nodes in each network to avoid possible location-related bias. In this section, we demonstrate the in-sample and out-ofsample results of solutions for NDP obtained by using different approaches. We test the NOBEL-US network with default parameter settings and examine the two cases when the 20-observation test samples contain regular and very low demand data, respectively. Table 5 presents the first-stage cost of allocating capacities on arcs (“1st-capaC”), the total in-sample objective cost (“in-Total”), the optimal objective value of the second-stage problem having network flow cost and unmet demand penalty (“in-2nd”), given by each model listed under “Approach”. We then fix the optimal solution of each approach and evaluate the out-of-sample performance using all 10,0 0 0 demand realizations in the reference set. We report the expected cost of sending ad hoc flows in “outE[Flow]” and the expected penalty cost of demand shortage in “outE[Pena]”. We sum up the two cost values and present the out-of-sample expected cost yielded by the optimal solution of each approach in “out-2nd”. We also calculate the gap between the in-sample estimated cost for the second stage (see “in-2nd”) and the out-ofsample cost realization (see “out-2nd”), and report the result in “I/O gap”. We make the following observations from Table 5. First, the results of stochastic NDP method can vary a lot depending on the assumed distribution type and the sample size | | used by the SAA

method. It could lead to higher or lower cost of allocating arc capacities in the first stage than DR-NDP, but always yields better insample estimated cost for the second stage compared with all the other methods except for the deterministic model. Second, robust NDP indeed takes a pessimistic view and results in both higher cost of allocating capacities and higher in-sample objective value for the second-stage problem as compared to all the other approaches. DR-NDP sits in between stochastic NDP and robust DNP, by providing adequately robust solutions. This is because DR-NDP uses more distributional information (e.g., σ i ) from data than robust NDP does and thus becomes less conservative. Third, as robust NDP prepares for the worst case, it yields lower expected penalty cost for unmet demand than DR-NDP, stochastic NDP, and deterministic benchmark. It yields higher expected flow cost, but still lower expected total cost in the second stage of the out-of-sample test. Lastly, DR-NDP has the most accurate estimate of the total cost under uncertainty as the differences between its in-sample and out-of-sample second-stage costs are the smallest for both regular and very low demand cases. On the other hand, stochastic NDP and the deterministic model are over optimistic, reflected by the large, positive percentage values given in “I/O gap”; robust NDP is over conservative, reflected by the large, negative percentage values in “I/O gap”. The stochastic NDP has much larger gaps between the two cost values if given very low demand observations; in the same case, robust NDP with = 1 can perform much better (i.e., improving the cost estimation gap from −52.19% to −12.68%), while the result of robust NDP with = 2 is not sensitive to the data set change. We further take the ratios of entries of each column for each approach to the corresponding ones of DR-NDP in Table 5, and report the results in Table 6. We can view the above trends more clearly by noting that (i) the two robust NDP models yields > 1 ratios for all the costs expect the expected penalty cost given by the out-of-sample test; (ii) the results of the stochastic NDP models vary significantly and the ratios are mostly < 1 for the in-sample costs, but become larger and even > 1 in the out-of-sample results. Therefore, we conclude that DR-NDP is more suitable for balancing in-sample and out-of-sample performance, and do not over

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

53

Table 5 Average in-sample and out-of-sample results of different approaches for default NOBEL-US instance with regular and very low demand samples. Approach

regular demand sample 1st-capaC

in-Total

in-2nd

outE[Flow]

outE[Pena]

out-2nd

I/O gap

drndp gam50 gam10 0 0 uni50 uni10 0 0 det rndp1 rndp2

2486.51 2173.45 2398.52 2500.26 2136.48 3113.02 3162.53 3342.32

5393.68 4788.81 4894.92 4983.97 5289.97 4156.15 8243.30 12225.79

2907.17 2615.36 2496.40 2483.71 3153.49 1043.14 5080.77 8883.47

567.11 510.17 551.30 567.79 492.50 665.64 671.04 689.31

2400.48 2767.21 2502.08 2384.68 2864.56 1801.21 1758.11 1632.48

2967.59 3277.38 3053.38 2952.47 3357.06 2466.85 2429.15 2321.79

2.08% 25.31% 22.31% 18.87% 6.46% 136.48% −52.19% −73.86%

drndp gam50 gam10 0 0 uni50 uni10 0 0 det rndp1 rndp2

very low demand sample 980.63 1746.19 1024.83 1709.94 835.73 1408.05 1551.42 902.22 945.32 1576.91 1031.77 1296.78 1393.37 6942.18 3328.60 11938.11

5765.55 685.10 572.32 649.20 631.59 265.01 5548.81 8609.51

241.54 247.41 196.72 216.65 228.97 249.61 314.57 656.01

5103.12 5004.27 5490.46 5277.33 5155.27 5158.59 4530.55 1912.74

5344.66 5251.68 5687.18 5493.98 5384.24 5408.20 4845.12 2568.75

−7.30% 666.56% 893.71% 746.27% 752.49% 1940.75% -12.68% −70.16%

Table 6 Average performance ratios of different approaches to DR-NDP for default NOBEL-US instance with regular and very low demand samples. Approach

regular demand sample 1st-capaC

in-Total

in-2nd

outE[Flow]

outE[Pena]

drndp gam50 gam10 0 0 uni50 uni10 0 0 det rndp1 rndp2

0.87 0.96 1.01 0.86 1.25 1.27 1.34

0.89 0.91 0.92 0.98 0.77 1.53 2.27

0.90 0.86 0.85 1.08 0.36 1.75 3.06

0.90 0.97 1.00 0.87 1.17 1.18 1.22

1.15 1.04 0.99 1.19 0.75 0.73 0.68

drndp gam50 gam10 0 0 uni50 uni10 0 0 det rndp1 rndp2

very low demand sample 1.05 0.98 0.12 0.85 0.81 0.10 0.92 0.89 0.11 0.96 0.90 0.11 1.05 0.74 0.05 1.42 3.98 0.96 3.39 6.84 1.49

1.02 0.81 0.90 0.95 1.03 1.30 2.72

0.98 1.08 1.03 1.01 1.01 0.89 0.37

address the issue of distributional ambiguity as robust NDP, by utilizing the empirical moments of data to construct the ambiguity set.

5.5. Results of sensitivity analysis We conduct sensitivity analysis by varying the unit penalty cost Gi , the number Ni of intervals for discretizing the quadratic objective of DR-NDP, and the number |D| of demand nodes in Sections 5.5.1–5.5.3, respectively. We test the following cases for sensitivity analysis: •

•

We set the unit penalty cost being 10, 20, 50 times the default penalty, i.e., being 10 0(|N| − 1 ), 20 0(|N| − 1 ), 50 0(|N| − 1 ), respectively. We vary the number of demand nodes from 2 to 5 and randomly pick demand nodes for each case. The true distribution of the demand at each node for generating the reference data set is assumed to follow a Gamma distribution with shape 4 × (3/|D|) and scale 5, in order to make the distribution of the total demand be equivalent among all the cases. The number of

•

supply nodes and the design of supply node capacity are kept the same. We change Ni to 5 and to 10 for the DR-NDP and adjust the cut points of intervals [di , di ], ∀i ∈ D accordingly.

5.5.1. Sensitivity to unit penalty cost Gi It is often difficult to determine the right amount of penalty for not meeting the demand and thus we analyze how varying the value of Gi , ∀i ∈ D can affect the optimal solution of each approach and their performance. In Fig. 1, we abstract the information and plot the average out-of-sample expected total cost (named “E[total cost]”) given by each approach as the percentage of the one using a regular demand sample. We report the results for cases where the test sample for each approach contains regular, low, and very low demand data. Each subfigure corresponds to a choice of Gi . In Fig. 1, when the observed demand is lower, the out-of-sample expected cost given by each approach increases except for the robust NDP and the cost is especially high for the very low demand case. For robust NDP, when we use the default penalty Gi = 10(|N| − 1 ), the expected cost based on lower demand observations is slightly higher, but it is lower when we increase the unit penalty to 10, 20, and 50 times the default setting. The increases of the out-of-sample cost of the stochastic NDP and the deterministic model are much more significant than the cost increase by using optimal solutions given by DR-NDP for Gi settings. Therefore, the DR-NDP approach is able to provide less variable results than the two approaches when we happen to observe the very low demand data, and meanwhile provides less conservative and cheaper design of arc capacities. We use NOBEL-US network with default parameter as an example and demonstrate the optimal solutions of DR-NDP under the four different penalty settings in Fig. 2. The shaded nodes 2, 6, 13 are supply nodes and the nodes 7, 12, 14 with dark circles are the demand nodes. As the unit unmet demand penalty Gi , ∀i ∈ D increases, we intend to build more arcs and allocate more capacities to flow more units to satisfy the total demand. We also build more arcs from the supply node that is further away from the three demand nodes (i.e., node 2 in this case) when Gi increases. Next, we compare the optimal solutions of DR-NDP, stochastic NDP (using “gam10 0 0”), and robust NDP (using “rndp2”) for Gi = 500(|N| − 1 ). (The solution patterns for other penalty settings are similar and thus we only use Gi = 500(|N| − 1 ) as a representative.) We plot the three solutions in Fig. 3. The robust NDP builds the most number of arcs with the highest capacities among the

54

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

Fig. 1. Relative out-of-sample performance for different settings of unit penalty Gi , ∀i ∈ D.

Fig. 2. Optimal solutions of DR-NDP for different settings of unit penalty Gi , ∀i ∈ D.

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57

55

Table 9 Average solution time of DR-NDP for default NOBEL-US instance with |D| = 2, 3, 4, 5.

|D| = 2 |D| = 3 |D| = 4 |D| = 5

demand

time (s)

2ndT (s)

#cuts

regular low very low regular low very low regular low very low regular low very low

2.202 1.661 0.507 3.729 4.122 1.286 16.007 3.298 2.188 20.327 9.868 2.018

2.157 1.634 0.491 3.663 4.087 1.248 15.906 3.254 2.138 20.178 9.771 1.947

55 57 54 74 77 75 100 87 106 122 135 125

three, to satisfy the demand at all the three nodes 7, 12, and 14. The stochastic NDP only uses supply nodes 6 and 13, but does not use node 2 even the unit penalty increases to Gi = 500(|N| − 1 ). The DR-NDP builds fewer number of arcs with smaller capacities than the robust NDP, but more than the stochastic NDP, to balance the cost of allocating capacities and flowing commodities and the penalty cost of demand shortage.

Fig. 3. Optimal solutions of DR-NDP, stochastic NDP, robust NDP with Gi = 500(|N| − 1 ).

Table 7 Average solution time of DR-NDP using Ni = 5, 10, 25 for default NOBEL-US instance. #Interval

demand

time (s)

2ndT (s)

#cuts

Ni = 5

regular low very low regular low very low regular low very low

0.207 0.321 0.032 0.533 0.752 0.274 3.729 4.122 1.286

0.178 0.299 0.029 0.497 0.729 0.245 3.663 4.087 1.248

42 62 16 51 58 59 74 77 75

Ni = 10

Ni = 25

5.5.2. Sensitivity to the number Ni of intervals used by DR-NDP In DR-NDP, we discretize the demand bound [di , di ] into Ni intervals with equal length and reformulate the quadratic objective function as a piecewise-linear function for computation. Intuitively, increasing Ni yields more accurate approximation but more variables and thus longer solution time. We analyze the sensitivity results of changing Ni and examine the tradeoff between CPU time and solution quality by using DR-NDP. Table 7 reports the total time, time for solving the separation problem, and the number of cuts generated for solving one DRNDP default instance when using Ni = 5, 10, 25 under different demand samples. We note that the CPU time increases almost linearly in Ni , while the majority of the solution time is again devoted to solving the separation problem. The number of cuts increases slightly as we increase Ni from 5 to 25. We then examine the solution quality change for different Ni choices. We report the average results of the in-sample and out-ofsample cost performance of DR-NDP for solving replications of the default NOBEL-US instance under Ni = 5, 10, 25 in Table 8. The results are similar for all Ni values we test, and Ni = 25 yields smaller capacity design cost but higher estimated cost for the secondstage problem as compared to Ni = 5. Using larger Ni can also yield smaller gap between the estimated in-sample cost and the out-ofsample expected cost (given by the difference between values in “out-2nd” and “in-2nd”). 5.5.3. Sensitivity to the number |D| of demand nodes Recall that in the default case, we have three demand nodes randomly distributed in each network. Here we vary the number of demand nodes and test |D| = 2, 3, 4, 5. We first show the CPU time of using DR-NDP for optimizing a NOBEL-US network instance with default settings and regular, low, and very low demand cases in Table 9. We observe that the total solution time increases much

Table 8 Average cost performance of DR-NDP using Ni = 5, 10, 25 for default NOBEL-US instance with regular demand observations. #Interval

in-Total

1st-capaC

in-2nd

out-Total

out-2nd

outE[Flow]

outE[Pena]

Ni = 5 Ni = 10 Ni = 25

5616.37 5745.72 5783.81

2200.60 1835.07 1908.88

3415.78 3910.65 3874.93

6125.48 6050.52 6051.83

3924.89 4215.45 4142.95

440.95 391.75 403.91

3483.93 3823.70 3739.04

56

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57 Table 10 Average cost performance of DR-NDP for default NOBEL-US instance with regular demand observations and |D| = 2, 3, 4, 5. Approach

1st-capaC

in-Total

in-2nd

out-Total

out-2nd

out-E[Flow]

out-E[Pena]

I/O gap

drndp

|D| = 2 |D| = 3 |D| = 4 |D| = 5

3006.80 1908.88 1888.44 2182.64

6135.23 5783.81 6407.90 5686.72

3128.43 3874.93 4519.46 3504.08

5996.30 6051.83 6474.37 5315.35

2989.50 4142.95 4585.93 3132.71

667.14 403.91 383.45 381.17

2322.36 3739.04 4202.48 2751.54

−4.44% 6.92% 1.47% −10.60%

gam10 0 0

|D| = 2 |D| = 3 |D| = 4 |D| = 5

3006.80 1767.99 2041.56 2286.91

5758.20 5357.18 5542.44 4997.29

2751.40 3589.19 3500.88 2710.38

5996.30 6020.60 6428.30 5177.80

2989.50 4252.61 4386.74 2890.89

667.14 391.48 428.63 420.41

2322.36 3861.13 3958.11 2470.48

8.65% 18.48% 25.30% 6.66%

det

|D| = 2 |D| = 3 |D| = 4 |D| = 5

3709.75 2116.96 2079.97 2338.66

5383.99 4957.86 5509.12 4213.22

1674.23 2840.90 3429.15 1874.56

6162.19 6080.36 6515.61 5414.68

2452.44 3963.40 4435.64 3076.03

758.92 436.71 411.59 405.25

1693.52 3526.69 4024.05 2670.78

46.48% 39.51% 29.35% 64.09%

rndp2

|D| = 2 |D| = 3 |D| = 4 |D| = 5

4776.00 3926.77 3952.95 3169.00

11059.20 16555.77 23074.80 28621.20

6283.20 12629.00 19121.85 25452.20

6833.88 7462.18 8430.01 7334.63

2057.88 3535.41 4477.06 4165.63

826.02 506.84 395.54 194.07

1231.86 3028.56 4081.52 3971.56

−67.25% −72.01% −76.59% −83.63%

more significantly if the given test sample consists of regular demand when we increase the number of demand nodes, as compared to if we are given low or very low demand data. This is due to the increased number of cuts (and thus iterations) needed for larger |D|. In Table 10, varying |D| = 2, 3, 4, 5, we report the average insample and out-of-sample cost performance of solutions given by different approaches for solving 10 replications of the default NOBEL-US instance with regular demand observations. We arrange all the columns similar to the ones in Table 5 and compute the gap between in-sample and out-of-sample costs using (out-2nd − in-2nd )/in-2nd × 100%. We only report results of “gam10 0 0” as a representative of the stochastic NDP approach and results of “rndp2” as a representative of the robust NDP, since the result trends are similar for the others. In Table 10, |D| = 2 has significantly higher capacity design cost than other choices of |D| for all the approaches. This is because given higher demand in two demand nodes, more arcs with much larger capacities are needed to connect all the three supply nodes to the two nodes; but when |D| = 3, 4, 5, the relatively lower demand in each node can be satisfied by one or two supply nodes and thus we can construct fewer arcs with smaller capacities. As a result, we see smaller expected flow cost in the out-of-sample test for |D| = 2. The comparison results among the four different approaches are similar to the ones of Table 5. Based on the column “I/O gap”, we observe that (i) the in-sample estimated cost is much closer to the out-of-sample cost realization in the reference data set for DR-NDP. The stochastic NDP and the deterministic model both under-estimate the second-stage cost, while the robust NDP significantly over-estimate the out-of-sample cost. Using DRNDP, we over-estimate the expected cost in the out-of-sample test when |D| = 2, 5, and under-estimate when |D| = 3, 4. However, the gap magnitudes are much smaller compared with the other three approaches. 6. Conclusions and future research In this paper, we proposed a DR-NDP approach to handle NDPs with random demand under distributional ambiguity. We developed a cutting-plane algorithm for solving DR-NDP, by leveraging the structure of the approximate reformulation to generate effective cuts. Our results showed that the solutions of DR-NDP were more conservative as compared to the traditional stochastic optimization approach and a deterministic model using nominal demand data, but much less conservative as compared to robust NDP us-

ing a budget uncertainty set. Crucially, the results of DR-NDP were insensitive to the changes of observed demand, the unit penalty cost setting for unmet demand, the number of demand nodes and their distribution in a network. In general, the DR-NDP can have much better prediction about the expected cost in reality (by testing its optimal solution in the reference data set), with smaller and less variable in-sample and out-of-sample cost gaps. Via examining networks with diverse sizes and structures, we conclude that DRNDP is more suited for niche uses, in which data is scarce, but yet the network planner cannot afford to have too much unmet demand due to “bad” data. For future research, we plan to investigate more variants of DR-NDP that utilize more statistical information of the random data for constructing the ambiguity set of the unknown distribution. Moreover, it will be instructive to apply DR-NDP to real-world problems. A possible application that could fit this niche is humanitarian relief supply network design (see, e.g., Hong et al., 2015). A decision maker may have very little data to react to a disaster and many parameters are random such as travel time and demand for relief. Meanwhile, the demands are urgent and need to be satisfied with high priority. These factors make such a problem a very suitable candidate for potential application. Acknowledgement The authors are grateful for all the referees and the Associate Editor for their constructive feedback and helpful suggestions. This research was funded in part by the National Science Foundation grant CMMI-1433066 and by the Department of Defense, Department of the Army grant W911NF-17-1-0102. References Altin, A., Yaman, H., Pinar, M.c., 2011. The robust network loading problem under hose demand uncertainty: formulation, polyhedral analysis, and computations. INFORMS J. Comput. 23 (1), 75–89. Álvarez-Miranda, E., Cacchiani, V., Lodi, A., Parriani, T., Schmidt, D.R., 2014. Single– commodity robust network design problem: complexity, instances and heuristic solutions. Eur. J. Oper. Res. 238 (3), 711–723. Atamtürk, A., Zhang, M., 2007. Two-stage robust network flow and design under demand uncertainty. Oper. Res. 55 (4), 662–673. Ayoub, J., Poss, M., 2016. Decomposition for adjustable robust linear optimization subject to uncertainty polytope. Comput. Manage. Sci. 13 (2), 219–239. Babonneau, F., Vial, J.-P., Klopfenstein, O., Ouorou, A., 2013. Robust capacity assignment solutions for telecommunications networks with uncertain demands. Networks 62 (4), 255–272. Ben-Ameur, W., Kerivin, H., 2005. Routing of uncertain traffic demands. Optimiz. Eng. 6 (3), 283–313. Ben-Tal, A., denHertog, D., De Waegenaere, A., Melenberg, B., Rennen, G., 2013. Robust solutions of optimization problems affected by uncertain probabilities. Manage. Sci. 59 (2), 341–357.

H. Nakao et al. / Computers and Operations Research 88 (2017) 44–57 Ben-Tal, A., Nemirovski, A., 1998. Robust convex optimization. Math. Oper. Res. 23 (4), 769–805. Bertsimas, D., Doan, X.V., Natarajan, K., Teo, C.-P., 2010. Models for minimax stochastic linear optimization problems with risk aversion. Math. Oper. Res. 35 (3), 580–602. Bertsimas, D., Sim, M., 2003. Robust discrete optimization and network flows. Math. Program. 98 (1), 49–71. Bertsimas, D., Sim, M., 2004. The price of robustness. Oper. Res. 52 (1), 35–53. Birge, J.R., Louveaux, F.V., 2011. Introduction to stochastic programming. Springer, New York, NY. Botton, Q., Fortz, B., Gouveia, L., Poss, M., 2013. Benders decomposition for the hop– constrained survivable network design problem. INFORMS J. Comput. 25 (1), 13–26. Cacchiani, V., Jünger, M., Liers, F., Lodi, A., Schmidt, D.R., 2016. Single-commodity robust network design with finite and hose demand sets. Math. Program. 157 (1), 297–342. Calafiore, G.C., El Ghaoui, L., 2006. On distributionally robust chance-constrained linear programs. J. Optim. Theor. Appl. 130 (1), 1–22. Delage, E., Ye, Y., 2010. Distributionally robust optimization under moment uncertainty with application to data-driven problems. Oper. Res. 58 (3), 595–612. Fortz, B., Poss, M., 2009. An improved benders decomposition applied to a multi-layer network design problem. Oper. Res. Lett. 37 (5), 359–364. Hong, X., Lejeune, M.A., Noyan, N., 2015. Stochastic network design for disaster preparedness. IIE Trans. 47 (4), 329–357. Høyland, K., Kaut, M., Wallace, S.W., 2003. A heuristic for moment-matching scenario generation. Comput. Optim. Appl. 24 (2–3), 169–185. Jiang, R., Guan, Y., 2016. Data-driven chance constrained stochastic program. Math. Program. Ser. A 158 (1), 291–327. Kleywegt, A.J., Shapiro, A., Homem-deMello, T., 2002. The sample average approximation method for stochastic discrete optimization. SIAM Journal on Optimization 12 (2), 479–502. Koster, A.M., Kutschka, M., Raack, C., 2013. Robust network design: formulations, valid inequalities, and computations. Networks 61 (2), 128–149.

57

Lee, C., Lee, K., Park, S., 2013. Benders decomposition approach for the robust network design problem with flow bifurcations. Networks 62 (1), 1–16. Magnanti, T.L., Wong, R.T., 1984. Network design and transportation planning: models and algorithms. Transportation Science 18 (1), 1–55. Mattia, S., 2013. The robust network loading problem with dynamic routing. Comput Optim Appl 54 (3), 619–643. McCormick, G., 1976. Computability of global solutions to factorable nonconvex programs: Part I–convex underestimating problems. Math. Program. 10 (1), 147–175. Orlowski, S., Wessäly, R., Pióro, M., Tomaszewski, A., 2010. Sndlib 1.0 - survivable network design library. Networks 55 (3), 276–286. Pióro, M., Medhi, D., 2004. Routing, flow, and capacity design in communication and computer networks. Elsevier. Poss, M., 2011. Models and algorithms for network design problems. 4OR 10 (2), 215–216. Poss, M., Raack, C., 2013. Affine recourse for the robust network design problem: between static and dynamic routing. Networks 61 (2), 180–198. Santoso, T., Ahmed, S., Goetschalckx, M., Shapiro, A., 2005. A stochastic programming approach for supply chain network design under uncertainty. Eur. J. Oper. Res. 167 (1), 96–115. Scutellà, M.G., 2009. On improving optimal oblivious routing. Oper. Res. Lett. 37 (3), 197–200. Shen, S., You, M., Ma, Y., 2016. Single-commodity stochastic network design under demand and topological uncertainties with insufficient data. Working paper. Soyster, A.L., 1973. Convex programming with set-inclusive constraints and applications to inexact linear programming. Oper. Res. 21 (5), 1154–1157. Ukkusuri, S.V., Patil, G., 2009. Multi-period transportation network design under demand uncertainty. Transp. Res. Part B 43 (6), 625–642. vanSlyke, R., Wets, R., 1969. L-shaped linear programs with applications to optimal control and stochastic programming. SIAM J. Appl. Math. 17, 638–663. Wiesemann, W., Kuhn, D., Sim, M., 2014. Distributionally robust convex optimization. Oper. Res. 62 (6), 1358–1376.