Implementation of a three-stage approach for the dynamic resource-constrained shortest-path sub-problem in branch-and-price

Computers & Operations Research 40 (2013) 385–394

Contents lists available at SciVerse ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

Implementation of a three-stage approach for the dynamic resource-constrained shortest-path sub-problem in branch-and-price Xiaoyan Zhu a,n, Wilbert E. Wilhelm b a b

Industrial and Information Engineering Department, University of Tennessee, 412 East Stadium Hall, Knoxville, TN 37996-0700, USA Industrial and System Engineering Department, Texas A&M University, TAMUS 3131, College Station, TX 77843-3131, USA

a r t i c l e i n f o

abstract

Available online 13 July 2012

The resource-constrained shortest-path problem (RCSP) is often used as a sub-problem in branch-andprice because it can model the complex logic by which many actual systems operate. This paper addresses two special issues that arise in such an application. First, RCSP in this context is dynamic in the sense that arc costs are updated at each column-generation iteration, but constraints are not changed. Often, only a few arc costs are updated at an iteration. Second, RCSP must be solved subject to arcs that are forbidden or prescribed as corresponding binary variables are fixed to 0 or 1 by the branching rule. A reoptimizing algorithm for dealing with a few arc-cost changes and a method for dealing with fixed arcs are proposed and incorporated into a three-stage approach, specializing it for repeatedly solving the dynamic RCSP as a sub-problem in branch-and-price. Computational tests evaluate the effectiveness of the proposed algorithms. & 2012 Published by Elsevier Ltd.

Keywords: Resource-constrained shortest-path problem Three-stage approach Re-optimization Fixed arcs Column generation Branch-and-price

1. Introduction The resource-constrained shortest-path problem (RCSP) is a variant of the shortest-path problem (SPP). Let R be a set of nonnegative, discrete-valued resources and T¼(T1,y,T9R9) be an 9R9dimensional vector specifying the limited supply of each resource. Considering an acyclic, directed graph G, let n ¼9V(G)9 and m¼9A(G)9, and order nodes topologically so that every arc (i,j)AA(G) satisfies ioj. Note that VðdÞ and AðdÞ denote the set of nodes and arcs of graph d, respectively. Each arc (i,j) has an associated cost cij and a vector uij ¼ ðuij1 ,    ,uij9R9 Þ, indicating the amount of each resource required to traverse the arc. Let path v1  vj denote a series of consecutive arcs from start node v1 to node vj (such a path may not be unique). The requirement of resource r (cost) on a path is the sum of the requirements of resource r (costs) associated with all arcs on that path. RCSP is to find the shortest path (i.e., the path with the least total arc cost) from start node v1 to end node vn with a total requirement of each type of resource rAR that observes a given limited supply Tr. RCSP is often used as a sub-problem in branch-and-price (B&P) [1] and has been successful in solving well-known problems like crew scheduling [2,3], prescribing the content and timing of products upgrades [4], optimizing picking and placing operations on dual-head placement machines [5–7], and multi-commodity

n

Corresponding author. Fax: þ1 865 9740588. E-mail addresses: [email protected] (X. Zhu), [email protected] (W.E. Wilhelm).

0305-0548/$ - see front matter & 2012 Published by Elsevier Ltd. http://dx.doi.org/10.1016/j.cor.2012.07.007

flow problems [8]. Many applications are modeled using an acyclic graph and/or multiple resource restrictions (e.g., [2–7] mentioned above, nurse scheduling [9], crew pairing [10], simultaneous scheduling of pilots and aircraft flights [11], and transfer line balancing [12]). A resource, for example, can be time, distance, capacity, money, workload, or reliability requirement. Motivated by these applications, Zhu and Wilhelm [13] studied RCSP with one or more resource constraints on an acyclic graph in the context of column generation in which arc costs change at each iteration but resource constraints do not. Zhu and Wilhelm [13] proposed a three-stage approach (TSA), analyzed its worst-case complexity, and evaluated it computationally, comparing its performance with that of CPLEX as well as the state-ofthe-art label-setting algorithm of Dumitrescu and Boland [14], showing that TSA is effective in repeatedly solving the dynamic RCSP (e.g., in column generation). When RCSP is used as a sub-problem in column generation, arc costs (i.e., objective function coefficients) are updated using new (optimal) values of dual variables from the master problem but resource constraints remain the same at each column-generation iteration. Consequently, RCSP must be reoptimized with respect to these new arc costs, giving rise to the dynamic property of RCSP. We refer to the dynamic RCSP as the underlying graph and resource constraints are fixed, but arc costs are changed dynamically at each column-generation iteration. Even though label setting algorithms can apply effective dominance criteria to restrict the number of alternative paths they analyze, they suffer from the need to solve the entire problem from scratch each time

386

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

Table 1 Acronyms. B&B MDFA RCSP SPP

Branch-and-bound Method for dealing with fixed arcs Resource-constrained SPP Shortest-path problem

arc costs change. In contrast, TSA deals only with resources in the first two stages, so they must be implemented only once in solving the overall problem, and the third stage, which is solved repetitively, is a shortest path problem in which arcs costs are changed at each iteration. The current paper extends this line of research, investigating two issues with the goal of effectively solving the dynamic RCSP in a B&P context. First, when RCSP is used as a sub-problem in column generation, it is typical that costs are updated for a different subset of arcs on each column-generation iteration. We propose a reoptimizing algorithm (ROA) to generate the new solution by updating the last solution using revised (i.e., different) arc costs and incorporate ROA in TSA. Second, when RCSP is used as a sub-problem in B&P, some arcs in the graph may be fixed: forbidden or prescribed (i.e., associated decision variables fixed to 0 or 1, respectively, by the branching rule). We propose a method for dealing with fixed arcs (MDFA) and incorporate it in TSA. The research objectives of this paper are (i) specialized algorithms (i.e., ROA and MDFA) for RCSP on an acyclic graph that can handle issues related to reoptimization as well as arcs that are fixed by branching within B&P; (ii) complexity analysis of these algorithms; (iii) computational evaluation; and (iv) implementation issues for solving the dynamic RCSP in B&P. These specialized algorithms are incorporated in TSA to make it suitable for solving the dynamic RCSP in B&P. Consequently, they can speed up B&P to solve problems (e.g., the scheduling problems in [2,3]) that involve RCSP as a sub-problem. It is worth mentioning that the goal of our computational tests is to evaluate the performances of ROA and MDFA imbed within TSA to solve the dynamic RCSP, rather than column generation or B&P itself. For reader convenience, Table 1 summarizes the acronyms used in this paper. The remainder of this paper is structured as follows. Section 2 reviews related literature and Section 3 gives an overview of TSA and related propositions. Section 4 proposes ROA and presents a version of TSA that incorporates ROA. Section 5 proposes MDFA and presents a version of TSA that incorporates MDFA. Section 6 presents a set of tests designed to evaluate the effectiveness of ROA and MDFA. Finally, Section 7 gives conclusions and discusses implementation issues of these algorithms in B&P.

2. Literature review RCSP is NP-hard [15], even if the graph is acyclic, only one resource constraint is involved, and all resource requirements and costs are positive [14]. RCSP on an acyclic graph is NP-hard in the ordinary sense and can be solved in pseudo-polynomial time [17]. The solution to RCSP is guaranteed to be elementary (no node is visited more than once) if arc costs are nonnegative and resource constraints have only upper bounds, or if the graph is acyclic as in this paper [16]. Because a number of important practical problems imbed RCSP, it has been studied rather extensively. Handler and Zang [15], Beasley and Christofides [16], and Mehlhorn and Ziegelmann [18] used methods that involve solving a relaxed problem using Lagrangian or linear relaxation. Others (e.g., [14,19–26]) used dynamic programming. These methods did not explicitly consider repeated solution of the dynamic RCSP as a sub-problem in

B&P OA ROA TSA

Branch-and-price Optimizing algorithm used in stage 3 of TSA Reoptimizing algorithm used in stage 3 of TSA Three-stage approach

column generation. In this context, the dynamic RCSP must be solved repeatedly, each time with an updated set of arc costs but the same set of resource constraints. To reoptimize at each column-generation iteration, prior algorithms must be employed starting from scratch. In contrast, several approaches exploit the column-generation context. Zhu and Wilhelm [13] proposed TSA for use in repeatedly solving the dynamic RCSP on an acyclic graph. Wilhelm et al. [4] reported a similar, but much less efficient, two-phase approach for a layered acyclic graph in which each arc is incident from a node in one level to another node in the next level. However, neither study specially investigated the issues of reoptimization and fixed arcs in the context of B&P. This paper proposes methods to deal with these two issues and specializes TSA, incorporating the proposed methods for repeatedly solving the dynamic RCSP. Section 3 gives a brief review of TSA so that this paper is self-contained. In early work related to reoptimizing SPP, Goto and Sangiovanni-Vincentelli [27] investigated the problem of updating shortest paths from all nodes to a selected set of nodes for the case in which the cost on each of a subset of arcs is decreased. Their method, which is based on LU factorization of arc cost matrix, {cij}, requires a considerable amount of memory and is effective only if matrix {cij} is sparse. Gallo [28] adapted Dijkstra’s shortest-path algorithm [29] for reoptimizing a SPP in two cases: (i) a different node is selected to be the start node of SPP, and (ii) exactly one arc is assigned a new cost that is less than the previous value. Fujishige [30] proposed another Dijkstra-like method for the case in which each of a set of arcs incident to one node is assigned a new cost that is less than the previous value. Recently, Buriol et al. [31] proposed a technique that can reduce the sizes of heaps used by several reoptimization algorithms [32–34] for the case in which a single arc is assigned a new cost that is either smaller or larger than the previous value. They presented a comprehensive computational evaluation of their technique and provided a survey of research that dealt with the case in which the cost of a single arc is changed. These methods are very restrictive and apply only to special cases (e.g., exactly one arc is assigned a new cost, or each of a set of arcs incident to one node is assigned a new cost). Although they may be applied iteratively if several arc costs change, the resulting algorithms are computationally impractical when a number of arc costs are assigned new values. Pallottino and Scutella [35] proposed a methodology to reoptimize SPP on an arbitrary graph for the case in which the cost on each of a specified subset of arcs is either increased or decreased. They generalized previous work [28,30], devising a two-phase method. The dual phase sequentially reoptimizes arcs with increased costs. The primal phase dynamically decomposes the set of arcs with decreased costs into disjoint subsets and reoptimizes subsets sequentially. As a sub-problem in column-generation, RCSP must be reoptimized at each column-generation iteration with respect to new arc costs but subject to the same set of resource constraints. The optimal solution (i.e., a shortest path tree rooted at start node) from the last iteration is available. This paper presents ROA to reoptimize the dynamic RCSP on an acyclic graph if the costs of any subset of arcs are changed to values that are either smaller or larger than the previous ones.

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

The issue related to fixed arcs arises when a B&P approach employs RCSP sub-problem(s). The traditional method induces the RCSP solution to exclude (include) the decision variable fixed to 0 (1) by assigning a large (small) cost to the corresponding forbidden (prescribed) arc (e.g., [36]). However, such a method has the same complexity whether arcs are forbidden and prescribed or not, even though these constraints can be used to reduce the size of RCSP to improve computational effectiveness. This paper presents a method (i.e., MDFA) that can detect infeasibility with respect to resource limitations before solving RCSP and reduce computational effort by exploiting the set of fixed arcs. MDFA is suitable for repeatedly solving the dynamic RCSP, subject to a set of fixed arcs at a node in the branch-and-bound (B&B) tree.

GE can be identified easily by tracing back using pðs1n Þ. If arc ðshi ,skj Þ is on the shortest path for SPP on GE, then arc (i,j) is on the optimal path for RCSP on G. In this paper, we propose an alternative method, ROA, which can be used in stage 3 (in place of OA) in column generation applications. For solving RCSP as a sub-problem in column generation, TSA is described as follows, and its effectiveness has been tested and verified in [13]. Initially, TSA runs stage 1 and 2 algorithms (only once) to generate GE; for a solution at each column-generation iteration with respect to a new set of arc costs, only stage 3 (i.e., OA) is used to solve SPP on GE and the solution for SPP is the solution for RCSP. Proposition 1 below [13] gives an upper bound on the size of GE and Proposition 2 [13] gives the complexity of each stage and of the entire TSA. Proposition 1. The number of nodes (arcs) in GE is bounded by nPrAR(Tr þ1) (mPrAR(Tr þ1)).

3. Description of the three-stage approach TSA comprises a preprocessing stage (stage 1), a setup stage (stage 2), and a solution stage (stage 3). Stages 1 and 2 deal with the resource constraint(s) and are invoked only once for a dynamic RCSP instance; only stage 3 deals with arc costs and must be applied for a solution with respect to each different set of arc costs. Thus, TSA is suitable for solving RCSP as a sub-problem in column generation. Stage 1 algorithm (preprocessing stage) reduces the size of an instance in two ways. First, it reduces input graph G to reduced graph (i.e., sub-graph) GR by deleting bottleneck nodes and arcs. An arc or a node is bottleneck if it cannot be on any v1 vn path that is feasible with respect to the resource limitation(s). A bottleneck arc (i,j) is one that renders the resource available at vj insufficient to continue a path to end node vn. A node is bottleneck if all incoming and outgoing arcs incident to it are bottlenecks. Second, the stage 1 algorithm formulates a tight window (i.e., a lower and upper bounds) for each resource at each node vj in reduced graph GR, making stage 2 of TSA work much better than if it used the trivial (vector) window [0,T] at each node. The computational evaluation in [13] showed that these preprocessing techniques can reduce computational effort, perhaps significantly, thus facilitating solution. Following stage 1, the stage 2 algorithm (setup stage) releases the resource constraints by expanding GR to form expanded graph GE. Let skj A VðGE Þ denote the kth node in GE that is associated with specific node vjAV(GR) in GR and Sj be the set of nodes skj that are associated with vj. Let ykj ¼ ðykj1 ,    ykj9R9 Þ be a cumulative resource requirement vector at skj . Basically, the stage 2 algorithm processes nodes vjAV(GR) in the order of increasing index, calculating ykj , associating a node skj ASj with each unique vector ykj , and connecting shi to skj if ykj is calculated from yhi . The specially-designed technique for calculating ykj makes the stage 2 algorithm computationally effective [13]. After stages 1 and 2, RCSP on G is reduced to SPP (from start node s11 to end node s1n ) on GE. The cost on arc ðshi ,skj ÞAA(GE) for shi A Si and skj A Sj is the cost on arc (i,j)AA(G), cij. The task of stage 3 (solution stage) is to optimize (or reoptimize) SPP on acyclic graph GE. Zhu and Wilhelm [13] proposed an optimizing algorithm (OA) for stage 3, as described below. Let pðskj Þ be the label on node skj , that is, the minimum cost among s11 skj paths in GE. Let pðskj Þ denote the predecessor of skj on the shortest path from s11 to skj (break ties arbitrarily). The values of pðskj Þ and pðskj Þ can be determined in increasing order of skj indices j and k as follows:

pðs11 Þ ¼ 0 and pðs11 Þ ¼ 0 ðstart node s11 has no predecessorÞ; n

n

pðskj Þ ¼ minpðshi Þ þ cij : ðshi ,skj Þ A AðGE Þ ¼ pðshi Þ þ ci j and pðskj Þ ¼ shi : n

387

n

n

ð1Þ ðs1n Þ,

which is also the The cost of an optimal path on GE is z ¼ p optimal solution value for RCSP on G. Given z, a shortest path on

Proposition 2. Initialization of TSA (preprocessing and setup stages) requires O(9R9mþ 9R9mPrAR(Tr þ1) þmSrAR(Tr þ1)) and O(9R9m2 þ9R9mPrAR(Tr þ1)þmSrAR(Tr þ1)) time for RCSP with 9R9¼1 and 9R94 1, respectively. Each solution (OA) requires O(mPrAR(Tr þ1)) time. 4. Reoptimizing algorithm In this section, we explore the opportunity to improve TSA for solving the dynamic RCSP, by using ROA in stage 3 of TSA in place of OA. Subsection 4.1 proposes ROA, which works to reduce computational effort when only a few arcs are assigned new costs, and analyzes its complexity. Based on that, Subsection 4.2 presents a version of TSA that incorporates ROA for solving RCSP as a sub-problem in column generation. 4.1. Description of ROA Rather than solving from scratch like OA to prescribe each solution, ROA only updates labels pðskj Þ that are affected by new arc costs. The shortest path tree found at the last iteration (i.e., using the previous arc costs) is available; it comprises the shortest path from start node s11 to each of the other nodes in GE. Let B DA(G) be a set of arcs that have updated (i.e., new) costs. These new costs may be smaller or larger than the old ones. Let cij denote the old cost on arc (i,j)AA(G); and c0ij , the new cost on the arc. Correspondingly, let pðskj Þ and p0 ðskj Þ denote the old and new minimum costs among s11 skj paths, respectively. Recall that pðskj Þ denotes the predecessor of skj in the shortest path tree. Let BSE ðskj Þ and FSE ðskj Þ denote the sets of the predecessors (backward stars) and successors (forward stars) of skj in GE, respectively. BSE ðskj Þ ¼ shi : ðshi ,skj Þ A AðGE Þ and FSE ðskj Þ ¼ shi : ðskj ,shi Þ A AðGE Þ. Additionally, let H be a heap that stores a set of arcs in G that have updated costs; and HE, a heap that stores a set of arcs in GE, in which the tail of each arc has the updated pðskj Þ value. Using this notation, Fig. 1 details ROA given set B. ROA initializes heap H with the arcs in B (step 1) and continues until H¼| and HE ¼| (step 2). On each iteration of step 2, the arc with the smallest index j from either H or HE is selected for processing (step 2(i)). If the selected arc is from H, then steps 2(ii– x) apply; otherwise, steps 2(xi–xviii) apply. Because arc (i,j) in G might correspond to several arcs (shi ,skj ) in GE, each with tail in Si and head in Sj, respectively, the processing of arc (i,j) in H involves dealing with all arcs (shi ,skj ) in GE (see step 2(iv)). In steps 2(v–viii), each arc (shi ,skj ) that corresponds to arc (i,j) in H is processed as follows (j here is designated as j1 in ROA of Fig. 1). If c0ij 4cij and pðskj Þ ashi (step 2(v)), arc (i,j) is not on the shortest path prescribed at the last column-generation iteration and will not affect the

388

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

Fig. 1. ROA.

shortest path since c0ij 4cij. If c0ij 4cij and pðskj Þ ¼ shi (step 2(vi)), then p0 ðskj Þ and pðskj Þ are updated to be minp0 ðshi Þ þc0ij : shi A BSE ðskj Þ and argminsh A BSE ðsk Þ p0 ðshi Þ þ c0ij , respectively. If c0ij ocij and pðskj Þ a shi i j (step 2(vii)), and, in addition, if p0 ðshi Þ þ c0ij o p0 ðskj Þ, then update 0 k 0 h 0 p ðsj Þ ¼ p ðsi Þ þ cij and pðskj Þ ¼shi ; otherwise, the current shortest path is still optimal. If c0ij ocij and pðskj Þ ¼ shi (step 2(viii)), then the current shortest path is still optimal but the optimal cost on the path is decreased by amount of cij c0ij . If pðskj Þ is changed during the processing of (shi ,skj ), all arcs outgoing from skj are inserted into HE (step 2(ix)). After it is processed, arc (i,j) is removed from H (step 2(x)). A similar analysis applies to steps 2(xiii–xvii) where pðski Þ differs from p0 ðski Þ, but we do not detail this part of ROA to conserve space. ROA is tailored for reoptimizing RCSP on an acyclic graph. It processes each node in Si before any node in Sj for j4i. Further, because no pair of nodes in set Si is connected, heap HE can be stored as an n-dimensional array in which each element is a linked-list of unordered arcs. All arcs with head in Si are inserted into linked-list HE[i] if labels (i.e., the minimum cost pðskj Þ) on their tails are changed by steps 2(ix) and (xvii). With heap structure HE, heap operations (extraction in steps 2(i) and (x) and insertion in steps 2(ix) and (xvii)) can be done in constant time. ROA processes arcs in HE[i] before processing arcs in HE[j] for j4i. The following proposition establishes the complexity of ROA. Proposition 3. ROA runs in O(9A(GE)9)¼ mPrAR(Tr þ1) time in the worst-case.

Proof. Since each operation on heap HE takes a constant amount of time, each arc in heap HE is processed in constant time. In the worst-case, each arc in GE must be processed once; thus, the total run time of ROA is O(9A(GE)9) ¼mPrAR(Tr þ1) according to Proposition 1. Note that ROA and OA are two methods for solving SPP in stage 3 of TSA. Even though they have the same worst-case complexity of mPrAR (Tr þ1), ROA may require longer run time since ROA must maintain a heap. On the other hand, the actual run time of ROA relates to set B, and increases with 9B9 on average. For small 9B9, ROA is more likely to be faster than OA. Section 6.1 presents computational tests that compare average-case performances of ROA and OA. 4.2. ROA imbedded in TSA Fig. 2 details a version of TSA that incorporates ROA in stage 3 to reoptimize RCSP as a sub-problem in column generation. In stage 3, l is a parameter specified by the analyst. After the first column-generation iteration, if the number of arc costs that are assigned new values on an iteration is greater than l, SPP is solved from scratch using OA; otherwise, SPP is reoptimized using ROA. l ¼0 implies OA is implemented at every iteration. Section 6 designs experiments to estimate an appropriate value for l for the instances we test. Fig. 2 shows that stages 1 and 2 are implemented only once in solving an instance of the dynamic RCSP. At each iteration, only stage 3 is needed and it runs in mPrAR(Tr þ1)

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

389

Fig. 2. TSA incorporating ROA.

Fig. 3. Dashed arcs ((a)–(c)) are forbidden due to prescribed arc (i,j).

time. Thus, TSA is suitable for repeatedly solving RCSP as a subproblem in column generation.

Proof. If arc (i,j) is prescribed, then any arc (i0 ,j0 ) identified by cases (a), (b), and (c) cannot be on a v1  vn path together with arc (i,j) (see Fig. 3).

5. Method for dealing with fixed arcs

Corollary 1. Suppose F1 a|. If prescribed arcs in F1 cannot be sorted as  ð2Þ ðik ,jk Þ,k ¼ 1,2,. . .9F 1 9 such that ik o jk and jk rik þ 1 ,

In the context of B&P, some arcs (i,j) in G may be fixed, posing opportunities to specialize TSA to gain effectiveness. Subsection 5.1 investigates properties related to fixed arcs on an acyclic graph and delineates MDFA to deal with them. Subsection 5.2 presents MDFA and a version of TSA for solving RCSP as a subproblem in B&P, and analyzes the complexity of MDFA. 5.1. Details of MDFA We start from graph GR which does not contain bottleneck arcs because stage 1 and 2 algorithms remove all bottleneck arcs (see Section 3) and make reduced graph GR and expanded graph GE available for use. We assume that binary variables corresponding to the bottleneck arcs have zero values and are eliminated permanently. Thus, the set of fixed arcs does not include the bottleneck arcs. Considering an arbitrary node in the B&B tree, let F0 (F1) be the set of arcs in GR that correspond to the binary variables fixed to 0 (1) by the branching rule. That is, F0 is the set of forbidden arcs which are not allowed on any path and F1 is the set of prescribed arcs that must be on the optimal path. At each B&B node, F0 and F1 are given and the corresponding RCSP sub-problem is subject to the constraints they induce. Because such an RCSP sub-problem at a B&B node must be solved repeatedly at each column-generation iteration, it is worth taking time to reduce GE using F0 and F1 before applying OA or ROA in stage 3 of TSA to solve the RCSP as a sub-problem in B&P. Define implied forbidden arcs in GR as the arcs that cannot be on v1  vn paths that contain all arcs in F1 and do not contain any arc in F0. Of course, forbidden arcs in F0 are also implied forbidden arcs. All arcs in GE corresponding to implied forbidden arcs in GR can be eliminated from GE. Theorem 1 identifies implied forbidden arcs that are induced by prescribed arcs in F1 Corollary 1 follows Theorem 1. Theorem 1. If arc (i,j) in acyclic graph g is prescribed (i.e., xij ¼ 1), then any arc (i0 ,j0 ) with (a) i0 oi and j0 4i, or (b) i0 ¼i and j0 aj, or (c) ioi0 oj is an implied forbidden arc.

then RCSP, subject to the set of prescribed arcs F1, is infeasible because prescribed arcs in F1 cannot be all on a same v1  vn path. Let F~ 0 be the set of all implied forbidden arcs in GR, including (i) ones in F0; (ii) ones induced by each prescribed arc (i,j)AF1: [ði,jÞ A F 1 ð[i0 o i,j0 4 i ði0 ,j0 ÞÞ [ ð[i0 ¼ i,j0 a j ði0 ,j0 ÞÞ [ ð[i o i0 o j ði0 ,j0 ÞÞ in which the three terms in parentheses are associated with cases (a), (b), and (c) of Theorem 1, respectively; and (iii) ones that are not on any v1 vn path after removing the forbidden arcs in (i) and (ii) from GR. Note that removing implied forbidden arcs in part (ii) is equivalent to fixing prescribed arcs. Given F0 and F1, the algorithm detailed in Fig. 4 generates F~ 0 . In Fig. 4, BS(vj) (FS(vj)) denotes the set of predecessors (successors) of vjAV(GR). Step 1 initializes F~ 0 with arcs in F0 and implied forbidden arcs in part (ii). Steps 2 and 3 identify implied forbidden arcs in part (iii). Step 4 checks if the remaining graph is connected. Now, using F~ 0 , the algorithm detailed in Fig. 5 identifies and removes all forbidden arcs in GE. Note that arc (i,j) A F~ 0 implies that arcs (shi ,skj ) with shi A Si and skj A Sj are forbidden in GE. We call the resulting graph the revised graph and denote it as GEF. Then, GEF is revised from GE by removing all forbidden arcs in GE, forming a subgraph of GE. Let BSEF ðskj Þ (FSEF ðskj Þ) denote the set of predecessors (successors) of node skj in GEF. BSEF ðskj Þ is a subset of BSE ðskj Þ and FSEF ðskj Þ a subset of FSE ðskj Þ. The algorithm in Fig. 5 applies to GE and first eliminates all forbidden arcs of GE that correspond to arcs in F~ 0 (step 1). Then, it identifies and removes the arcs in the remaining graph that cannot be on any s11 s1n path (steps 2 and 4). If BSEF ðs1n Þ ¼ | (step 3) or FSEF ðs11 Þ ¼ | (step 5), GEF is disconnected and so STOP; otherwise, return GEF, which is defined by BSEF ðskj Þ and FSEF ðskj Þ for skj A VðGE Þ (step 6). 5.2. Description of MDFA and its imbedding in TSA Fig. 6 details a version of TSA that incorporates MDFA for solving the dynamic RCSP in B&P. At the B&B root node (line 1), F1 ¼F0 ¼|, and RCSP is solved in lines 2–6, indeed the version of

390

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

Fig. 4. Generating F~ 0 algorithm.

Fig. 5. Revising GE to GEF algorithm.

Fig. 6. TSA incorporating MDFA (lines 8–10) within B&P.

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

TSA in Fig. 2. At each other node in the B&B tree (line 7), some arcs are fixed, and MDFA (specified in lines 8–10) first sorts arcs in F1 according to expression (2) if F1 a| (see Corollary 1), then generates F~ 0 , and finally revises GE to form GEF. After revising GE to GEF, either OA or ROA (lines 11–13) can be applied to GEF to find an optimal solution for each column-generation iteration. In ROA, set B should be replaced with B\ F~ 0 because arcs in F~ 0 do not have associated arcs in GEF. MDFA detects infeasibility caused by fixed arcs in three steps: sorting arcs in F1 detects infeasibility with respect to prescribed arcs (Corollary 1 and Fig. 6), generating F~ 0 detects infeasibility due to GR being disconnected (Fig. 4), and revising GE to form GEF detects infeasibility due to GE being disconnected (Fig. 5). Proposition 4. MDFA runs in O(9A(GE)9) ¼mPrAR(Tr þ1) time. Proof. In Fig. 6, MDFA operations are in lines 8–10. Line 8, which sorts arcs in F1 according to expression (2), can be done in O(9A(G)9)¼ O(m) time; line 9 generates F~ 0 and runs in O(9A(G)9)¼ O(m) time; and line 10 revises GE to form GEF and runs in O(9A(GE)9) time because each arc in GE is processed in constant time. Thus, MDFA runs O(9A(GE)9)¼mPrAR(Tr þ1) time according to Proposition 1. By Propositions 2–4, MDFA has the same complexity as OA and ROA. Thus, incorporating MDFA into TSA does not change its complexity. Note that stage 1 and 2 algorithms in TSA (lines 2 and 3) are invoked only once at the root node of B&B tree for the entire problem while MDFA is implemented once for each other B&B node. After using MDFA to deal with fixed arcs, OA or ROA (i.e., stage 3 of TSA) can obtain the solution at each columngeneration iteration in O(9A(GEF)9) ¼O(9A(GE)9) time. Remark 1. Define the span of arc (i,j) as j  i, assuming that the nodes are numbered topologically. Based on Theorem 1, it is appropriate to select a branching variable xij whose corresponding arc (i,j) has a large span, so that it will generate more implied forbidden arcs and cause GEF to be smaller than if the branching variable corresponds to an arc of a small span. 6. Computational evaluation This section describes the computational tests that we design to evaluate the effectiveness of ROA and MDFA imbed into TSA for solving the dynamic RCSP. The acyclic graphs are generated such that arc (i,j) appears in the graph with probability p for j defined as each integer on ½iþ 1,minðn,i þ qÞ. Specifying parameter q(1oqrn) to restrict the span of arc (i,j): j  irq, the expected number of arcs in a graph that is randomly generated in this manner with parameters of n, p, and q is pq(n (qþ1)/2). Note that we test instances of RCSP with one resource limitation constraint (i.e., 9R9¼1) and do not test instances with multiple constraints (i.e., 9R9 41) because ROA and MDFA are implemented after generating expanded graph GE and are thus not related to the number of resources. We classify test instances in two classes. For class s1, resource requirements on arcs are independently generated from discrete uniform distribution DU[1,10], and for class s2, resource requirement on arc (i,j) is a function of its span and equals [R(j  i)] where ½d denotes the nearest integer and R is a random number from continuous uniform distribution U(0.0,1.0). Resource limitation T is determined by T ¼ T min þðT max T min Þ  Z,

ð3Þ

where parameter Z on ð0:0,1:0 is used to control the tightness of a resource limitation. Resource requirement T min (T max ) denotes the minimum (maximum) amount of resource required by any v1  vn

391

path. These values can be obtained by setting cij ¼uij (cij ¼ uij) and implementing a classical SPP algorithm on input graph G. For all instances, we draw arc costs independently from U(  100.0, 100.0). Table 2 defines each of all four test instances (s1–500, s2–500, s1–1000, and s2–1000) in columns 1–9. Column 1 gives the code (a combination (class, n ¼9V(G)9)) that identifies the test instance. Columns 2 and 3 give parameters q and p, respectively, which we use to generate G. Column 4, 9A9, gives the number of arcs in input graph G. Columns 5–7 give T min , T max , and Z, respectively. Columns 8 and 9, 9VE9 and 9AE9, give, respectively, the numbers of nodes and arcs in expanded graph GE. We program all algorithms in C/Cþþ and conduct all experiments on a 3.2 GHz Pentium IV PC with 512 Mb of RAM. Sections 6.1 and 6.2 describe our computational tests of ROA and MDFA, respectively.

6.1. ROA Tests We evaluate ROA by comparing it with OA. Recall that B is the set of arcs in G that are assigned new costs. Since the stage 2 algorithm expands G (strictly speaking, GR), the number of arcs in expanded graph GE that are assigned new costs may be much larger than 9B9. Our experiment tests four values of 9B9 (1, 3, 5 and 10) on each of four instances, each with a specified Z ¼{0.1,0.5,0.9}. We run 100 replications for each test to simulate 100 column-generation iterations; each involves solving an RCSP with new arc costs specified by set B using both OA and ROA. For each replication, we randomly choose 9B9 arcs from G with equal likelihood for each arc. For a selected arc, we generate a random value D from U((  100.0, 1.0)[(1.0,100.0)) and add it to the current cost of the arc. If D is negative, the cost of the arc decreases; otherwise, it increases. Table 2 presents computational results in columns 10–16. The run times in Table 2 do not include times for stages 1 and 2 of TSA because OA and ROA are used only in stage 3 of TSA. Note that the run time for OA does not depend on set B and is almost constant over 100 replications; thus, for OA, column 11 presents the average run time per replication over 100 replications of each test. However, the run time for ROA depends on the number of arcs in B and their locations in the graph. Columns 12– 15 give the average, standard deviation, minimum, and maximum values of run times per replication over 100 replications. It shows that ROA run times differ significantly as a function of the arcs in set B, even for the same value of 9B9. For all tests, the minimum run time of ROA can be 0 s if the arcs in B do not alter the current shortest path tree. However, the maximum run time of ROA can be much larger than OA, if the arcs in B affect a large portion of GE (extremely if the arcs in B affect the entire graph). Consequently, standard deviations of ROA run times are large, even much larger than their average values. Furthermore, the last column in Table 2 gives the ratio of the average run time for one replication of ROA to that of OA. It shows that, on average, if 9B9 r3, ROA is faster than OA, except for instance s2-500 with Z ¼0.9, which has the smallest expanded graph (9484 nodes and 11,173 arcs) among all test instances. However, on average, if 9B9Z5, ROA is slower than OA for most instances except for instance s2-1000 with Z ¼0.1, 0.5, or 0.9. These results are consistent with expectations (see the discussion following Proposition 3). Basically, ROA may save time by resolving only portions of GE that are affected by the arcs in B. In general, the portion that is affected increases with 9B9, and may even encompass the entire graph. On the other hand, ROA must maintain an additional heap to store the portion of GE that has been affected and heap operations require additional run time.

392

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

Table 2 Comparison of OA and ROA. (0) test

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

(1) instance

s1-500

s2-500

s1-1000

s2-1000

(2) q

125

125

240

240

(3) p

0.1

0.1

0.05

0.05

(4) 9A9

5450

5450

10435

10435

(5) Tmin

9

38

14

49

(6) Tmax

494

491

615

962

(7)

Z

(8) 9VE9

(9) 9AE9

(10) 9B9

0.1

15004

152044

0.5

44860

526350

0.9

9484

11173

0.1

11604

57040

0.5

44531

469984

0.9

11104

125874

0.1

38268

378134

0.5

111871

1239797

0.9

21039

239205

0.1

37378

177380

0.5

173034

1761164

0.9

40906

447040

This trade-off balances at some value of 9B9, on average, around 3–5 for most of the tests. 6.2. MDFA Tests We evaluate MDFA by comparing it with a traditional method for dealing with fixed arcs, which we use as a benchmark. The traditional method [36] assigns a large cost to each forbidden arc and a small cost to each prescribed arc and then applies OA with respect to these adjusted arc costs. This subsection studies instance s1-1000 of 9A9¼10,435 arcs with Z ¼ 0.1, 0.5, and 0.9, using a factorial experiment design with four specified sets of F1 and two specified sets of F0, as shown in columns 4–7 in Table 3. Again, we solve each test (with a specified (Z,F1,F0) setting) for 100 randomly generated replications. At each replication, we generate a new set of arc costs randomly from U( 100.0, 100.0). Table 3 presents computational results. Columns 1–3 have the same meaning as those in Table 2.

1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10 1 3 5 10

(11) OA run time

(12) (13) (14) ROA run time (  10  3 s)

(15)

(  10  3 s)

mean

stdev

min

max

10.32 10.31 10.62 10.31 30.01 30.15 29.99 30.15 6.87 6.88 6.72 6.72 5.47 6.09 5.01 5.15 28.59 28.59 28.76 28.59 7.97 7.82 7.96 7.97 24.22 24.54 24.06 24.06 70.62 70.16 71.25 72.34 12.19 12.03 12.34 12.18 13.6 13.14 12.19 12.03 109.53 109.54 110.47 109.38 22.66 21.73 22.35 21.87

2.02 2.96 10.93 17.66 9.85 24.53 93.28 152.02 1.08 7.18 18.43 29.69 1.09 1.88 5.62 7.84 7.82 15.62 59.84 116.87 1.72 3.82 10.15 19.53 1.87 7.65 25.63 39.53 8.29 29.21 97.81 251.42 1.88 5.31 17.35 44.54 1.25 3.44 9.06 12.81 15 28.13 103.9 349.07 2.02 9.53 16.09 56.88

6.88 6.528 24.993 33.586 55.121 79.526 268.432 310.131 3.958 34.411 59.442 62.745 3.995 5.602 10.798 12.08 47.27 52.036 168.558 232.379 6.251 6.766 35.046 38.727 5.548 22.729 101.882 82.868 50.864 107.244 332.752 562.565 7.437 19.904 71.108 90.218 4.798 7.863 19.116 20.488 75.526 79.775 357.198 742.203 6.518 20.153 48.137 93.649

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

47 31 171 235 531 531 1656 1688 16 250 344 344 16 32 63 63 469 359 1078 1094 47 31 296 188 31 157 922 719 500 781 2282 2765 62 172 563 484 31 31 125 109 656 703 2844 3485 47 140 422 437

(16) Ratio

0.196 0.287 1.029 1.713 0.328 0.814 3.110 5.042 0.157 1.044 2.743 4.418 0.199 0.309 1.122 1.522 0.274 0.546 2.081 4.088 0.216 0.488 1.275 2.450 0.077 0.312 1.065 1.643 0.117 0.416 1.373 3.476 0.154 0.441 1.406 3.657 0.092 0.262 0.743 1.065 0.137 0.257 0.941 3.191 0.089 0.439 0.720 2.601

Column 8 denotes the cardinality of F~ 0 , which is calculated by the algorithm in Fig. 4. Columns 9–12 give the number of nodes and arcs in GEF, 9VEF9 and 9AEF9, and ratios 9VEF9/9VE9 and 9AEF9/9AE9, respectively. It shows that 9F~ 0 9 is large for all tests (9F~ 0 9 Z 4,050) although 9F19 and 9F09 are either 1 or 4, respectively, small numbers in comparison with 9A9¼10,435. A large 9F~ 0 9 results in a small GEF. As shown in Table 3, when F1 ¼ S1 ¼{(130,161), (369,382),(546,593),(830,865)} and F0 ¼ S3 ¼{(145,195), (390, 451), (582,707), (922,960)}, 9F~ 0 9 ¼ 9,496; that is, 91% of the arcs in G are forbidden (9A9 ¼10,435). Further, 9F~ 0 9 is much larger when F1 ¼{(273,504)} than when F1 ¼{(503,504)} (6870 versus 4050 for F0 ¼{(145,195)}; 6875 versus 4057 for F0 ¼S3), because the span of arcs (273,504) and (503,504) are 231 and 1, respectively; the former is obviously much larger than the latter. Tests with 9F19¼4 (i.e., F1 ¼S1 and F1 ¼S2 ¼{(382,503), (503,504), (504,512), (512,546)}) give similar results. Next, compare tests using F1 ¼{(273,504)} and F1 ¼ S2. Although the former prescribes only one arc and the latter prescribes four arcs, 9F~ 0 9 is larger for

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

393

Table 3 Results of instance s1-1000 with fixed arcs. (0) test (1)

(2) 9VE9 (3) 9AE9

Z

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

0.1

38,268

(4) 9F19

378,134 1

4

0.5

111,871 1239197 1

4

0.9

21039

239205

1

4

(5) F1

(6) 9F09

{(273,504)} 1 4 {(503,504)} 1 4 S1 1 4 S2 1 4 {(273,504)} 1 4 {(503,504)} 1 4 S1 1 4 S2 1 4 {(273,504)} 1 4 {(503,504)} 1 4 S1 1 4 S2 1 4

(7) F0

(8) 9F~ 0 9

(9) 9VEF9

(10) 9AEF9

(11) 9VEF9/ 9VE9

(12) 9AEF9/ 9AE9

(13)a run time (s)

(14)b (15)c (16)d (17) breakeven

{(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3 {(145,195)} S3

6870 6875 4050 4057 9495 9496 5504 5511 6870 6875 4050 4057 9495 9496 5504 5511 6870 6875 4050 4057 9495 9496 5504 5511

13427 13423 29381 29372 9157 9157 20183 20181 35048 35033 85662 85625 33280 33280 54910 54877 2947 2947 13905 13884 3504 3504 8442 8421

85310 85172 215113 214830 24590 24587 123131 122942 317905 317275 703631 702372 106433 106428 432095 431202 9708 9703 118796 118419 10391 10388 49620 49461

0.35 0.35 0.77 0.77 0.24 0.24 0.53 0.53 0.31 0.31 0.77 0.77 0.30 0.30 0.49 0.49 0.14 0.14 0.66 0.66 0.17 0.17 0.40 0.40

0.23 0.23 0.57 0.57 0.07 0.07 0.33 0.33 0.26 0.26 0.57 0.57 0.09 0.09 0.35 0.35 0.04 0.04 0.50 0.50 0.04 0.04 0.21 0.21

2.40 2.38 2.37 2.37 2.38 2.38 2.41 2.46 6.99 6.98 6.97 6.99 6.93 6.97 7.02 6.97 1.09 1.06 1.08 1.07 1.08 1.06 1.11 1.06

0.16 0.16 0.14 0.14 0.11 0.11 0.14 0.16 0.56 0.50 0.36 0.38 0.34 0.34 0.42 0.42 0.11 0.11 0.09 0.09 0.06 0.04 0.11 0.11

0.76 0.76 1.51 1.53 0.34 0.31 0.93 0.99 2.03 2.05 4.38 4.37 1.03 1.03 2.76 2.72 0.18 0.14 0.67 0.69 0.17 0.15 0.38 0.34

0.92 0.92 1.65 1.67 0.46 0.42 1.07 1.15 2.54 2.54 4.74 4.75 1.37 1.37 3.18 3.14 0.29 0.25 0.77 0.78 0.23 0.20 0.49 0.45

10 10 17 17 6 6 10 11 11 11 14 15 6 6 10 10 13 12 23 23 7 6 16 16

S1 ¼{(130,161), (369,382), (546,593), (830,865)}; S2 ¼ {(382,503), (503,504), (504,512), (512,546)}; S3 ¼ {(145,195), (390,451), (582,707), (922,960)}. a

Run time of traditional method for 100 replications. Run time for generating F~ 0 and revising GEF. Run time for solving SPP on GEF for 100 replications. d Sum of run times in columns 14 and 15. b c

{(273,504)} than for S2, because the span of arc (273,504) is 231, larger than the total span of path 382-503-504-512-546 formed by the arcs in S2, which is 546  382 ¼164. These observations verify Remark 1, which points out that many implied forbidden arcs can be identified if the spans of prescribed arcs are large. Note that the value of 9F~ 0 9 is not as sensitive to the set of prescribed arcs as it is to the set of forbidden arcs. Columns 13–16 are specified by the legend at the bottom of Table 3. Stage 3 of TSA uses OA to solve SPP on GEF. Column 16 gives the total time for MDFA and OA for solving SPP on GEF for 100 replications. MDFA requires run time to identify and remove forbidden arcs from GE (column 14), but this ‘‘set up’’ time is offset by reducing the size of GE significantly (the graph is reduced more than 43% up to 96% (test 21 and test 22) in terms of the number of arcs (column 12)), consequently, reducing run time that OA requires (column 15). This ‘‘setup’’ time is incurred once in revising GE but a time saving that results from the smaller graph accrues at each replication. For 100 replications, the total time for MDFA (column 16) is less than that for the traditional method (column 13) for all tests. The last column in Table 3 gives the breakeven number of replications, which is calculated using break even ¼



 100  cpuðMDFAÞ , cpuðtraditional methodÞcpuðOAÞ

ð4Þ

where dde denotes the ceiling function. cpu(traditional method), cpu(MDFA), and cpu(OA) are the run times in columns 13–15, respectively. If the number of replications is larger than breakeven, the total run time for MDFA (revising GE to GEF) plus solving SPP on GEF is less than that of the traditional method; otherwise, the traditional method is faster. breakeven (column 17) ranges from 4 to 23 for this set of tests, implying that MDFA is very effective.

7. Conclusions and discussion This paper proposes ROA for reoptimizing RCSP at each column-generation iteration (except the first-time solution) for which only a few arcs are assigned new costs, and devises MDFA to deal with fixed arcs when RCSP is used as a sub-problem in B&P. Our computational results show that incorporating ROA and MDFA in TSA can effectively reoptimize the dynamic RCSP even with fixed arcs. ROA and MDFA are designed with the basic goal of repeatedly solving the dynamic RCSP, for example, as a subproblem in B&P. At the same time, we demonstrate the flexibility of TSA in that it can incorporate specialized algorithms to deal with crucial issues (e.g., reoptimizing and fixed arcs). It is worth mentioning that implied forbidden arcs must receive correct treatment in selecting branching variables. Consider a node in the B&B tree with specified sets F1, F0, and computed F~ 0 . Decision variables corresponding to implied forbidden arcs in F~ 0 should not be selected as a branching variable because, if one were selected, the left child node (fixing the branching variable to 0) would be the same as its parent node, which already forbids this arc, and the right child node (fixing the branching variable to 1) would be infeasible since it prescribes the implied forbidden arc. In the implementation of B&P, F~ 0 should be stored in association with each active node in the B&B tree to avoid such ineffectiveness. Note that the construction of F~ 0 can be expedited at each B&B node by starting with the F~ 0 constructed at its parent node. The research in this paper shows that algorithms for solving sub-problems in column generation (e.g., RCSP) should explicitly consider issues related to repeated solutions. Moreover, if the algorithms are used to solve a sub-problem in B&P, they should include effective methods to deal with fixed variables. The branching strategy should be considered together with the algorithms for solving sub-problem(s) to obtain good overall

394

X. Zhu, W.E. Wilhelm / Computers & Operations Research 40 (2013) 385–394

performance. For example, if the sub-problem is RCSP and we use TSA incorporating MDFA to solve it, it is appropriate to choose a branching variable whose corresponding arc has a large span (see Remark 1) and it is not good to choose one that has an associated arc that can be relegated to set F~ 0 . Our literature review describes numerous applications that involve RCSP as a sub-problem in a B&P context so that the ideas and methods proposed in this paper can be applied to solve such real-world problems as crew scheduling, nurse scheduling, process planning for electronic assembly machines, and multi-commodity transportation in supply chains. We leave it to future research to test our methods in such applications.

Acknowledgment This material is based upon work supported by the Texas Advanced Technology Program under Grant no. 000512-02482001 and by the National Science Foundation under Grant no DMI-0217265. References [1] Wilhelm WE. A technical review of column generation in integer programming. Optimization Eng 2001;2(2):159–200. [2] Mingozzi A, Boschetti MA, Ricciardelli S, Bianco L. A set partitioning approach to the crew scheduling problem. Operations Research 1999;47:873–88. [3] Fahle T, Junker U, Karisch SE, Kohl N, Sellmann M, Vaaben B. Constraint programming based column generation for crew assignment. Journal of Heuristics 2002;8(1):59–81. [4] Wilhelm WE, Damodaran P, Li J. Prescribing the content and timing of product upgrades. IIE Transactions 2003;35:647–64. [5] Wilhelm WE, Arambula I, Choudhry NN. A model to optimize picking operations on dual-head placement machines. IEEE Transactions on Automation Science and Engineering 2006;3:1–15. [6] Wilhelm WE, Choudhry ND, Damodaran P. A model to optimize placement operations on dual-head placement machines. Discrete Optimization 2007;182:626–39. [7] Wilhelm WE, Zhu X. Enabling flexibility on a dual head placement machine by optimizing platform-tray-feeder picking operations. Flexible Services Manufacturing Journal 2009;21(1–2):1–30. [8] Holmberg K, Yuan D. A multicommodity network-flow problem with side constraints on paths solved by column generation. INFORMS Journal on Computing 2003;15:42–57. [9] Jaumard B, Semet F, Vovor T. A generalized linear programming model for nurse scheduling. European Journal of Operational Research 1998;107:1–18. [10] Desaulniers G, Desrosiers J, Dumas Y, Marc S, Rioux B, Solomon MM, et al. Crew pairing at Air France. European Journal of Operational Research 1997;97:245–59. [11] Stojkovic´ M, Soumis F. An optimization model for the simultaneous operational flight and pilot scheduling problem. Management Science 2001;47: 1290–305. [12] Dolgui A, Guschinsky N, Levin G. A special case of transfer lines balancing by graph approach. European Journal of Operational Research 2006;168: 732–46.

[13] Zhu X, Wilhelm WE. A three-stage approach for resource-constrained shortest-path problem on an acyclic graph. Computers and Operations Research 2011;38:1210–8. [14] Dumitrescu I, Boland N. Improved preprocessing, labeling and scaling algorithms for the weight-constrained shortest path problem. Networks 2003;42:135–53. [15] Handler GY, Zang I. A dual algorithm for the constrained shortest path problem. Networks 1980;10:293–310. [16] Beasley JE, Christofides N. An algorithm for the resource constrained shortest path problem. Networks 1989;19:379–94. [17] Garey MR, Johnson DS. Computers and intractability: a guide to the theory of NP-completeness. San Francisco: Freeman; 1979. [18] Mehlhorn K, Ziegelmann M. Resource constrained shortest paths. In: Proceedings of the 7th Annual European Symposium on Algorithms (ESA 2000), LNCS 1879; 2000. pp. 326–37. [19] Joksch HC. The shortest route problem with constraints. Journal of Mathematical Analysis and Applications 1966;14:191–7. [20] Hassin R. Approximation schemas for the restricted shortest path problems. Mathematics of Operations Research 1992;17:36–42. [21] Aneja YP, Aggarwal V, Nair KPK. Shortest chain subject to side constraints. Networks 1983;13:295–302. [22] Desrochers M, Soumis F. A reoptimization algorithm for the shortest path problem with time windows. European Journal of Operational Research 1988;35:242–54. [23] Desrochers M, Soumis F. A generalized permanent labeling algorithm for the shortest path problem with time windows. Information Systems and Operational Research 1988;26:193–214. [24] Dumitrescu I, Boland N. Algorithm for the weight constrained shortest path problem. International Transactions in Operational Research 2001;8:15–30. [25] Feillet D, Dejax P, Gendreau M, Gueguen C. An exact algorithm for the elementary shortest path problem with resource constraints: application to some vehicle routing problems. Networks 2004;44:216–29. [26] Ioachim I, Ge´linas S, Soumis F, Desrosiers J. A dynamic programming algorithm for the shortest path problem with time windows and linear node costs. Networks 1998;31:193–204. [27] Goto S, Sangiovanni-Vincentelli A. A new shortest path updating algorithm. Networks 1978;8:341–72. [28] Gallo G. Reoptimization procedures on shortest path problems. Rivista di Matematica per le Scienze Economiche e Sociali 1980;3:3–13. [29] Dijkstra EW. A note on two problems in connection with graphs. Numerische Mathematik 1959;1:269–71. [30] Fujishige S. A note on the problem of updating shortest paths. Networks 1981;11:317–9. [31] Buriol LS, Resende MGC, Thorup M. Speeding up dynamic shortest path algorithms. Technical report TD-5RJ8B. AT&T Labs Research. Florham Park, NJ; 2003. [32] Ramalingam G, Reps T. An incremental algorithm for a generalization of the shortest-path problem. Journal of Algorithms 1996;21:267–305. [33] King V, Thorup M. A space saving trick for directed dynamic transitive closure and shortest path algorithms. In: Proceedings of the 7th Annual International Computing and Combinatorial Conference COCOON, LNCS 2108; 2001. p. 268–77. [34] Demetrescu C. Fully dynamic algorithm for path problems on directed graphs. PhD thesis. Department of Computer and Systems Science, University of Rome ‘‘La Sapienza’’; 2001. [35] Pallottino S, Scutella MG. A new algorithm for reoptimizing shortest paths when the arc costs change. Operations Research Letters 2003;31:149–60. [36] Jaumard B, Semet F, Vovor T. A two-phase resource constrained shortest path algorithm for acyclic graphs, Les Cahiers du GERAD G-96-48. E´cole des Hautes E´tudes Commerciales, Montre´al, Canada, H3T 2A7; 1996.