Universally maximum flow with piecewise-constant capacities
Universally Maximum Flow with Piecewise-Constant Capacities
Lisa K. Fleischer Graduate School of Industrial Administration, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, Pennsylvania 15213
A maximum flow over time generalizes standard maximum flow by introducing a time component. The object is to send as much flow from source to sink in T time units as possible, where capacities are interpreted as an upper bound on the rate of flow entering an arc. A related problem is the universally maximum flow, which is to send a flow from source to sink that maximizes the amount of flow arriving at the sink by time t simultaneously for all t ≤ T . We consider a further generalization of this problem that allows arc and node capacities to change over time. In particular, given a network with m arcs and n nodes, arc and node capacities that are piecewise constant functions of time with at most k breakpoints and a time bound T , we show how to compute a flow that maximizes the amount of flow reaching the sink in all time intervals (0,t ] simultaneously for all 0 < t ≤ T , in O (k 2mn log (kn 2/mm )) time. The best previous algorithm requires O (nk ) maximum-flow computations on a network with (m + n )k arcs and nk nodes. Our algorithm improves on this by a factor of O (nk ). © 2001 John Wiley & Sons, Inc. Keywords: dynamic flow; transshipment; parametric flow
1. INTRODUCTION In the 1960s, Ford and Fulkerson introduced flows over time to add time to the standard network flow model. Since then, flows over time have been used widely to model network-structured, decision-making problems over time: problems in electronic communication, production and distribution, economic planning, cash flow, job scheduling, and transportation. For examples, see the surveys of Aronson [4] and Powell et al. [26]. Flows over time have been previously referred to as dynamic network flows [4, 7, 12, 16, 17, 19, 20, 21, 22, 23, 26, 28, 30]. However, this term causes confusion about the problem being considered: “Dynamic” is more
Received November 1, 1998; accepted June 1, 2001 Correspondence to: L. K. Fleischer; e-mail:[email protected] Contract grant sponsor: AAUW; contract grant info: NSF; grant no. EIA-0049084 c 2001 John Wiley & Sons, Inc.
NETWORKS, Vol. 38(3), 115–125 2001
consistently used for a problem with input that changes online, over time. In many of the earlier papers, the network is fixed. While we allow for the network to change in this paper, all the changes are specified up front. In this paper and the earlier papers, it is the solution that changes over time. For this reason, we use the more descriptive, albeit slightly awkward, terminology of flows over time. A maximum flow over time generalizes standard maximum flow by the addition of time. A standard network consists of a set of nodes V with |V| = n and a set of arcs E with |E| = m, which is a subset of V × V. The capacity function from the arcs to the real numbers bounds the amount of flow allowed on each arc. In a flow-overtime network, the capacity function u limits the rate of flow into an arc. In addition, a flow-over-time network has a transit-time vector % associated with the arcs. The transit time of an arc is the amount of time it takes for flow to travel from one end of the arc to the other. Ford and Fulkerson [8] considered maximum flows over time: given a flow-over-time network with a specified source and sink, determine the maximum amount of flow that can be sent from source to sink in T time units. They show that the problem can be solved in polynomial time by using information obtained from a minimum-cost flow computation in a related network of comparable size. A generalization of this, a universally maximum flow, is a flow that simultaneously maximizes the amount of flow reaching a specified sink node d, sent from a specified source node s, by time t, for all 0 < t ≤ T. Wilkinson [30] and Minieka [19] showed that such a flow exists, and they both provided algorithms to solve this problem, but these do not run in polynomial time. There is no known polynomial time algorithm that solves the universally maximum-flow problem. It is not even known if the optimal solution is polynomial in the size of the input. Hoppe and Tardos [14, 15] presented an O(m−1 (m + n log n) log U) algorithm that computes a flow over time with the property that the quantity of
flow reaching the sink by time t is within (1 − ) of the maximum possible for all 0 ≤ t ≤ T. Another flow-over-time problem is the problem with multiple sources and sinks. A transshipment over time is a single-commodity flow that sends specified supplies located at source nodes to satisfy specified demands at sink nodes within a given time bound T. A universally quickest transshipment is a transshipment over time that satisfies as much demand as possible in every interval of time (0, t] for 0 < t ≤ T. If there is more than one sink, a universally quickest transshipment may not exist, so the universally quickest transshipment problem refers to a transshipment problem with one sink only. This paper considers the special case of flow problems over time when all transit times are zero—in other words, when flow is instantaneous. However, finite capacities still limit the rate of flow. In this setting, the universally maximum-flow problem is solved by sending flow at the rate of a maximum flow in the static network at every moment of time in (0, T]. (The static network is the interpretation of the flow-over-time network as a standard network. Transit times are ignored, and capacities are interpreted as bounds on the total flow allowed on an arc.) A universally quickest transshipment with zero transit times has a more complicated structure. This problem models the situation when the supplies and demands exceed network capacity, and it is necessary to send flow over time, in rounds. Unlike standard flow problems with multiple sources and sinks, this problem cannot be modeled as an equivalent s − t maximum flow over time. The reason such a transformation does not work here is because the arc capacities are upper bounds on the rate of flow per time period and do not correspond to the total amount of flow an arc may carry throughout an interval. Transshipment over time with zero transit times was discussed in [6, 12, 20, 29]. Hajek and Ogier [12] described the first polynomial time algorithm to find a universally quickest transshipment when all transit times are instantaneous. Their algorithm uses n maximum-flow computations on the underlying static network. In [6], this is improved upon by an algorithm that runs in the same asymptotic time as a preflow-push maximum-flow algorithm. All the above-mentioned problems may be generalized to allow storage of flow at the nodes. In this setting, there may be node-storage capacities that limit the total excess flow allowed at a node. Flow deficit is not allowed. In the above-mentioned problems, even when node excess is allowed, there always exists an optimal solution that does not use any node storage, and the algorithms mentioned above find such a solution. For a universally maximum flow with piecewise-constant capacity functions, the subject of this paper, this will not be the case: Node storage will be used in an optimal solution, and the amount of available storage will affect the solution.
116
NETWORKS–2001
Many other generalizations of this problem have been considered in the literature; however, few other generalizations are known to be polynomial-time solvable. One exception is an integer transshipment over time with general transit times, for which Hoppe and Tardos [16] gave the only polynomial-time algorithm known to solve this problem. Their algorithm repeatedly calls an oracle to minimize submodular functions. Orlin [23] provided polynomial- time algorithms for some infinite horizon problems, where the objective function is related to average throughput. Burkard et al. [5] found the minimum feasible time to send a specified flow from source to sink faster than incorporating the Ford–Fulkerson algorithm in a binary search framework. Most other work has focused on characterizing the structure of problems with more general capacity functions and determining when optimal solutions exist [1, 2, 3, 24, 25, 27, 28]. None of these problems are known to be polynomial-time solvable. Among the earlier work, Anderson et al. [2] considered the problem of finding a maximum flow in a network with zero transit times and with capacities and node storage that are much more general functions of time. They developed a duality theory for continuous network flows. We consider the generalization of the zero transit time, universally maximum-flow problem that allows both arc capacities and node-storage capacities to be piecewise-constant functions on [0, T] with at most k breakpoints. As with the general universally maximumflow problem, it is not clear that a universally maximum flow exists in these circumstances. Ogier [21] proved that it does and provided a polynomial-time algorithm that finds a universally maximum flow with at most nk breakpoints. These breakpoints can be computed with nk maximum-flow computations on a network with nk vertices and (m + n)k arcs. After the breakpoints are computed, additional maximum-flow computations on the same network are used to calculate the flow between two breakpoints. Thus, the total run time of Ogier’s algorithm is determined by the time to solve O(nk) maximum-flow problems on a network with nk vertices and (n + m)k arcs. As Ogier [21] demonstrated, this problem also generalizes the universally quickest transshipment problem. The main contribution of this paper is to (a) recognize that these problems can be solved by solving a parametric maximum-flow problem on a suitably defined graph and (b) to generalize the parametric maximum-flow algorithm of Gallo, Grigoriadis, and Tarjan [9] to fit the needs of this problem. The end result is that all the computations described in [21] can be performed in the same asymptotic time as one preflow-push maximumflow computation on a network with nk vertices and (n + m)k arcs. This improves the previous strongly polynomial run time by a factor of O(nk).
1.1.
The Model
The universally maximum-flow problem with piecewise-constant capacities (UMFP) is defined on a network with n nodes V and m arcs E, along with a nonnegative, left-continuous, piecewise-constant capacity function u : (E ∪ V) × (0, T] → Z+ ∪ {0}. When referring to the capacity of a particular arc e and a particular node i at a particular time t, we use the notation ue (t) and ui (t), respectively. The vector of all capacities at time t is simply denoted u(t), the vector of just arc capacities is denoted uE (t), and the vector of node capacities is denoted uV (t). For convenience of notation, we assume that the breakpoints of the capacity functions are exactly one unit apart. If this is not the case, we can scale time between each breakpoint to create an equivalent problem with this property; that is, if the time between two breakpoints is r, then we create an equivalent problem by changing this to 1 and multiplying all capacities in this interval by r. Any flow that solves this modified problem can be transformed back to a solution to the original problem by dividing the rate of flow sent in the new time interval by r and sending it over an interval that is r times longer. In such a manner, we scale T so that it equals k. Thus, for the remainder of the paper, we assume that our time bound is k, all k breakpoints are in the set {1, 2, . . . , k}, the domain of u is (E ∪ V) × (0, k], and u(t) = u(dte) for all 0 < t ≤ k. Let x be a flow over time. An arc e = (i, j) may be represented as ij in the context of a subscript: Here, xe and xij denote the same quantity, and similarly, for ue and uij . If (i, j) ∈ / E, then xij and uij are automatically set to zero. For A ⊂ V, we use the notation uA,A (t) to P denote u i∈A,j∈A ij (t) and the notation xA,A (t) to denote P [x i∈A,j∈A ij (t) − xji (t)]. We say that x is feasible if xe is a Lebesgue-measurable function on (0, k] for all e ∈ E, x obeys the arc capacity constraints xe (t) ≤ ue (t) for all e ∈ E, and all t ∈ (0, k], and the excess at each node, defined by Z τ xV,j (t)dt, pj (τ) := 0
obeys the node capacity constraint 0 ≤ pj (τ) ≤ uj (τ) for all j ∈ V − {s, d}, and all 0 < τ ≤ k. Denote the set of feasible flows over time by D. The value of flow x at time τ, denoted vτ (x), is the total flow reaching the sink by time τ, that is, Z τ xN,d (t)dt = pd (τ). vτ (x) :=
For each τ ∈ (0, k], we define a time-expanded graph Gτ on k copies of N, labeled N1 , . . . , Nk . Note that the vertex set is the same for all τ. Denote this vertex set by V∗ . Let iθ be the copy of vertex i in copy Nθ and eθ be the copy of arc e in Nθ . In Gτ , the arcs in Nθ for θ ∈ {1, . . . , dτe} are assigned capacities uE (θ); the arcs in Nθ for θ > dτe are given capacity 0; for θ = 1 to dτe−1, the holdover arcs leaving Nθ are given capacities equal to uV (θ); and for θ ≥ dτe, the holdover arcs leaving Nθ are assigned capacity 0. By this construction, we have that Gτ = Gσ if θ < τ ≤ σ ≤ θ + 1. Thus, there are exactly k unique graphs Gτ , one for each τ = θ ∈ {1, . . . , k}. The arcs with nonzero capacity in Gτ correspond to the unit intervals of time in which the corresponding arcs in the static network N may carry flow in any flow completing by time τ. Figure 1 depicts G2 and G2.5 for the simple network N on four nodes. The thin arcs correspond to arcs with 0 capacity. The thick arcs in Nθ correspond to arcs with capacity equal the capacity of the corresponding arc in N at time θ. The thick holdover arcs leaving Nθ have capacity equal to the capacity of the node they leave at time θ. The difference between G2 and G2.5 is that the arcs in N3 and the holdover arcs leaving N2 have nonzero capacity in G2.5 . This captures the fact that a flow completing by time τ = 2.5 may send some flow in N during the interval (2, 2.5], while a flow completing by time τ = 2 may not. 1.2.
Optimality Conditions
A τ-maximum flow is a feasible flow over time x that maximizes vτ (x). As in the static-flow setting, we can define a notion of time-cut that proves the optimality of a τ-maximum flow. A time-cut is a function C that maps (0, k] to subsets of vertices of the network N such that C(α) = C(β) for all α, β ∈ (θ − 1, θ], θ ∈ {1, . . . , k} and C(α) ∩ {s, d} = {d} for all α ∈ (0, k]. Note that this defines the cut as containing the sink. This is the complement of what is traditionally defined as a cut. Let C := {C|C is a time-cut}. Since the range of a timecut is a finite set, and over the domain (0, k], each timecut can change only k times (each is uniform over each interval (θ, θ + 1]), C is a finite set. The value wτ (C) of time-cut C at time τ is defined below: It is the sum of capacities of arcs with their head in C and tail in C for some interval of time, times the length of this interval (the first term), plus the sum of capacities of nodes at the last integer time they are in C (the second term).
0
In this paper, we will use a series of related timeexpanded graphs. A time-expanded network of N is a directed graph that contains a copy of N for each time unit and holdover arcs from a copy of a node at time θ to the copy of the same node at time θ + 1 for θ ∈ Z.
Z wτ (C) =
0
τ
uC(t)C(t) (t)dt +
dτe−1 X
X
θ=1 j∈C(θ+1)∩C(θ)
uj (θ). (1)
A time-cut has an alternate interpretation as a set of arcs entering a cut in an appropriately defined timeNETWORKS–2001
117
FIG. 1. Examples of time-expanded graphs Gτ and Hτ for τ = 2 and 2.5.
expanded network. This will be made more explicit in the corollaries and discussion following Theorem 1. A τ-minimum cut is a time-cut C that minimizes wτ (C). Theorem 1 states that the value of a τ-maximum flow equals the value of a τ-minimum cut. This is comparable to strong duality in standard network flows. The proof of even the weak duality case is a little more involved than is the corresponding static version and is omitted. A slightly weaker statement was proved by Anderson et al. [1], but for more general capacity functions. Philpott [24] extended this work to handle more general transit times. Theorem 1 (Ogier [21]). For all τ, minC∈C wτ (C) = maxx∈D vτ (x). The following corollary shows that for integral values of θ a θ-maximum flow corresponds to a maximum flow in Gθ . Corollary 2. If θ ∈ {1, 2, . . . , k}, then a maximum flow in the time-expanded network Gθ corresponds to a θmaximum flow xθ of the same value, where the flow on arc et in Gθ is precisely the flow on arc e in interval (t − 1, t], that is, for all t ∈ {1, . . . , θ}, the flow on et corresponds to xeθ (σ) for all σ ∈ (t − 1, t]. Proof. A maximum flow in Gθ gives rise to a minimum cut Rθ ⊂ V∗ that contains all copies of the sink node in Gθ and no copies of the source node. Letting C(σ) = C(dσe) := Ndσe ∩ Rθ for all σ ∈ (0, k], we have that C is a time-cut of a value equal to the maximum flow
118
NETWORKS–2001
in Gθ and, hence, equal to the value of the corresponding flow over time described in the corollary. We extend the series of time-expanded graphs Gτ to a more continuous version Hτ , so that we may extend Corollary 2 to apply more generally to τ ∈ (0, k]. For each real 0 < τ ≤ k, we define Hτ to be a graph on the same vertex set and arc set as Gdτe , but with a slightly different capacity function. In Hτ , the capacities of the arcs in Ndτe are multiplied by τ − bτc, which is always less than one. The capacities of all other arcs equal the capacities of the arcs in Gτ . Denote the capacity function of Hτ by w˜ τ . Note that Hθ = Gθ for all θ ∈ {1, . . . , k}. Examples of Hτ and Gτ are depicted in Figure 1 for a simple four-node network N with k = 4. Note that the capacity of arcs in N3 in H2.5 is half of the capacity of the corresponding arcs in G2.5 . This modification in H2.5 models the fact that a flow that completes by time τ = 2.5 uses arcs in N only for half of the unit interval (2, 3]. The value of a time-cut at time τ corresponds to the capacity of the arcs entering the sink side of a cut in Hτ . We make this explicit as follows: Since the vertex set is the same for all Hτ and equals V∗ , the vertex set of Gτ , we also use V∗ to denote the vertex set of Hτ . A cut in Hτ is a set R ⊂ V∗ such that R contains all copies of and no copies of the source, that Skthe sink nodeS k is, R ∩ θ=1 {sθ , dθ } = θ=1 {dθ }. ∗ For a set R ⊂ V , define w˜ τ (R) to be the sum of capacities of arcs in Hτ entering R. We can equate each time-cut C ∈ C with a cut R ⊂ V∗ of Hτ such that wτ (C) = w˜ τ (R) by setting R such that R ∩ Nθ = C(θ) for θ ∈ {1, . . . , k}. Let Dτ be the set of feasible flows over time that are obtained from flows in Hτ by setting the
rate of flow in arc e in interval (θ − 1, min{τ, θ}] to be the flow on arc eθ in Hτ . The following corollary follows easily from the above discussions: Corollary 3.
maxx∈Dτ vτ (x) = minR⊂V∗ w˜ τ (R).
By the submodularity of cut functions in a graph, the nonempty intersection of minimum cuts is a minimum cut. Thus, there is a minimum cut in Hτ with the smallest sink-side. Define Rτ to be the minimum cut such that Rτ ⊆ R for all R such that w˜ τ (R) = w˜ τ (Rτ ). The following lemma describes properties of Rτ as a function of τ:
Lemma 4 (Ogier [21]). For all 0 ≤ τ < σ ≤ k, Rτ ⊆ Rσ , and Rτ is a left continuous function of τ. We give some intuitive explanation for the behavior of Rτ between breakpoints of u. For τ close to bτc, the arcs in Ndτe have very little capacity in Hτ . Thus, a minimum cut in Hτ may have arcs entering the cut in Ndτe with large u(dτe) value to avoid holdover arcs in Hτ from Nbτc to Ndτe with large capacity. One way to avoid these holdover arcs is to exclude all nodes of Ndτe except ddτe . As τ increases, the contribution of the capacities of arcs in Ndτe toward the value of the
FIG. 2. (Top left) The time expanded graph G4 for k = 4 displaying all input capacities. (Second row, left) The graph H3.5 with capacities at half-value in N4 . (Top right) Cuts Rτ for τ = 1, 2, 3, 3.5, 4, demonstrating the behavior described in Corollary 4 and the subsequent discussion. (Second row, right) A τ-maximum flow xτ for τ = 3.5. Notice that this flow saturates the cut R3.5 in H3.5 . (Bottom left) the values of q. (Bottom right) the flows x0 and x∗ .
NETWORKS–2001
119
cut increases, and, thus, the minimum cut may shift to a smaller capacity cut of Ndτe , thus including more nodes of Ndτe and possibly increasing the total capacity of holdover arcs contributing to the value of the cut. Exactly this situation is depicted in Figure 2. Figure 2 shows Rτ for τ = 1, 2, 3, 3.5, 4, for the same original network N described in Figure 1. Here and elsewhere, Rτ+ indicates the limit of Rσ as σ ↓ τ. The figure also shows a flow in H3.5 . Note that this flow saturates the cut R3.5 and is thus a 3.5-maximum flow. Let W be the set of breakpoints of Rτ as a function of τ, that is, W is the set of values τ such that Rτ ≠ Rτ+ . Ogier’s main contribution is to prove that there is an optimal solution x∗ to UMFP that is uniform on the intervals between successive breakpoints in W and that such an optimal flow can be defined piecemeal by special τ-maximum flows [21]. This is a theorem of Ogier reprinted here as Theorem 6. Part of Theorem 6 can be explained by complementary slackness conditions. We first provide this explanation and then give the precise technical statement. We start with two easy observations: (1) For σ ≤ θ / Rσ . (2) An optimal flow x∗ for UMFP is a − 1, iθ ∈ τ-maximum flow for all τ ∈ (0, k] and thus saturates Rτ in Hτ for all τ. We now consider the three cases (A–C) that determine the flow on arc (i, j) in interval (θ − 1, θ]. / Rσ and jθ ∈ Case A. If there is a cut Rσ such that iθ ∈ Rσ , then arc (iθ , jθ ) is saturated in a σ-maximum flow and thus must be saturated in an optimal flow up until time ∗ (t) = uij (t) for t ∈ (θ − 1, min{θ, σ}]. A σ. Precisely, xij / Rσ and iθ ∈ Rσ symmetric argument shows that if jθ ∈ ∗ then xij (t) = 0 for t ∈ (θ − 1, min{θ, σ}]. Case B. If σ ≤ θ is the earliest time that iθ ∈ Rσ+ , then we know from the first observation that σ > θ − 1. Lemma 5 below implies that for all σ < τ ≤ θ there is a flow-saturating arc entering Rτ with excess pi (t) = 0 for all t ∈ (σ, θ]. Since the optimal flow x∗ must saturate arcs entering Rτ up until time τ, this implies that in an optimal flow x∗ we need not use any node capacity at i in (σ, θ]. Thus, between breakpoints σ ≤ σ1 < σ2 ≤ θ, the flow on all arcs in Nθ ∩ Rσ2 for t ∈ (σ1 , σ2 ] may be determined by the maximum flow in Nθ with sources Nθ \Rσ2 and sink dθ .
Lemma 5. For i ∈ V − {s, d}, if σ is the earliest time that iθ ∈ Rσ+ , then for all σ < τ ≤ θ, there is a flow in Dτ saturating arcs entering Rτ with excess pi (t) = 0 for all t ∈ (σ, θ]. Proof. Key to the proof of this lemma is that all flows in Dτ for τ ≤ θ correspond to flows in Hθ , which contains no holdover arcs for nodes in Nθ . The value of the cut function at Rσ and Rσ+ for τ ∈ (θ − 1, θ] is determined by arcs in Nθ which are parameterized by τ − θ + 1 and constant capacity arcs that are either
120
NETWORKS–2001
holdover arcs for times ψ ≤ θ − 1 or arcs in Nψ for ψ ≤ θ − 1. Thus, we can express these two values as (τ − θ + 1)aσ + bσ and (τ − θ + 1)aσ+ + bσ+ . For τ = σ, these two values are equal. By definition of Rσ and Rσ +, we have that aσ > aσ+ and bσ < bσ+ . Thus, 0 < bσ+ − bσ = (σ − θ + 1)(aσ − aσ+ ). In words, all flow that enters Rσ+ − Rσ along constant capacity arcs either leaves Rσ+ − Rσ along constant capacity arcs or drains from Rσ+ − Rσ by time σ. Since no positive, constant capacity arc leaving Rσ+ − Rσ is a holdover arc for any time after θ−1, this implies that no node capacity is used at any node in Nθ ∩ (Rσ+ − Rσ ) between time σ and time τ, where τ is the next breakpoint. In the interval (τ, θ], the only additional flow entering this set comes from the parametrized arcs in Nθ and thus is bounded from above by aσ . Since this flow saturates Rθ and drains by time θ, it drains in Nθ without requiring node storage. Case C. We are left with the case that σ > θ−1 is the last value such that both iθ and jθ are not in Rσ and for all τ > σ both iθ and jθ are in Rτ . In this case, since the flow crossing Rσ saturates all arcs crossing Rσ , the flow on all arcs on the source side of the cut should be enough to saturate Rσ . Thus, if iθ and jθ are not separated by some σ ∗ (t) = xij (t) for t ∈ (θ − 1, min{σ, θ}]. time-cut, then xij With this motivation in mind, we follow [21] and define xτ to be a τ-maximum flow that is constant on each interval (θ−1, θ] for θ ∈ {1, 2, . . . , bτc} and also constant on interval (bτc, τ]. In addition, let the excess function pτ of xτ satisfy pjτ (τ) = 0 for all j ∈ V − {s, d}. That such a flow exists is implied by Corollary 3: Compute a maximum flow in Hτ . This flow saturates the arcs entering Rτ and, hence, has the value of a τ-maximum flow. It does not use any node storage arc leaving N(dτe) and, hence, completes by time τ. This is depicted in Figure 2 for τ = 3.5. We also define x0 to be the flow over time such that for all t, f = x0 (t) is a maximum flow for static network N with capacities u(t) = u(dte), sources Ndte \Rt , and sink d. Figure 2 describes x0 for an example. For each node j of the network N, define max{τ|τ ≥ t, jdte ∈ Rτ }, if this set is nonempty, (2) qj (t) := 0, otherwise. In words, qj (t) is the largest value of τ for which there is no path from jdte to ddτe in the residual network of a τ-maximum flow. We use qj (t) to determine which of the above cases A–C applies. If qi (t) > qj (t): This implies that if idte has a path in the residual flow graph of xσ to ddσe then jdte has a residual path to ddσe and there is some τ for which idte does not have a residual path to ddτe when jdte does. For dte − 1 < τ, this is Case A for the interval (dte − 1, min{τ, dte}]. This is Case B for the interval (τ, dte].
Theorem 6 (Ogier [21]). Flow x∗ , defined by τ xij (t), if qi (t) = qj (t) = τ [ Case B if τ = 0; else Case C ] ∗ xij (t) = [ Case A ] uij (t), if qi (t) > qj (t) [ Case A ] 0, if qi (t) < qj (t) is a τ-maximum flow for all τ ∈ (0, T]. Notice that the only τ-maximum flows used to compute x∗ are for τ ∈ W, the set of breakpoints. Thus, we may compute a universally maximum flow by first computing W, the set of breakpoints, then computing x0 , Rτ , and xτ for all τ ∈ W, determining q, and then constructing x∗ as given by Theorem 6. For the example in Figure 2, the breakpoints are W = {1, 2, 3, 3.5, 4}. The corresponding flows, time-cuts, values of q, and final flow x∗ are all depicted in the figure. Once W, the set of all breakpoints of Rτ , along with the corresponding cuts Rτ , τ ∈ W, are known, then x0 can be computed with |W| maximum flows in a network the same size as N, and x∗ can be computed with an additional |W| maximum flows, one in each Hτ , for all τ ∈ W. Ogier described a method to compute W in a piecewise fashion, computing W ∩ (θ − 1, θ] for θ = 1, . . . , k. Define the minimum-cut function κ by setting κ(τ) equal to the value of the minimum cut in Hτ , that is, κ(τ) = w˜ τ (Rτ ). Figure 3 describes the relation between the Rτ and κ. Lemma 7. Within the domain (θ − 1, θ] for θ ∈ {1, . . . , k}, κ is concave and piecewise linear with at most n − 2 breakpoints. Proof. For τ ∈ (θ − 1, θ], θ ∈ {1, . . . , k}, κ(τ) is the lower envelope (the minimum function) of {w˜ τ (Rσ )}σ , a set of linear functions of τ. This implies that κ is concave. Since Rτ is increasing in τ by Lemma 4, changes in the slope of κ(τ) as τ increases imply a change in the cut Rτ defining κ and thus an increase in |Rτ ∩Nθ |. This / Rτ for all τ together with the fact that dθ ∈ Rτ and sθ ∈ implies that κ(τ) has at most n − 2 breakpoints in a unit interval. Ogier [21] computed each subset W ∩ (θ − 1, θ] by finding τ such that w˜ τ (R(θ−1)+ ) = w˜ τ (Rθ ), computing Rτ , and recursing on the subintervals (θ−1, τ] and (τ, θ] until each interval does not properly contain elements of W. The following lemma justifies this approach: Lemma 8. If τ is such that w˜ τ (Rσ1 ) = w˜ τ (Rσ2 ) for σ1 < σ2 , then σ1 ≤ τ ≤ σ2 and either τ is a breakpoint of κ or both (σ1 , τ) and (τ, σ2 ) contain breakpoints of κ. Proof. For a fixed cut Rσ and for τ in the interval between two breakpoints of u, w˜ τ (Rσ ) is a linear function
FIG. 3. The relationship of κ to w˜ τ and Rτ . Here, τ0 is such that w˜ τ0 (Rτ1 ) = w˜ τ0 (Rτ4 ).
with positive slope, that is, it has form τaσ + bσ . Since Rσ is the minimum cut when τ = σ, the slope of w˜ τ (Rσ ) as a function of τ decreases as σ increases. Thus, τ such that w˜ τ (Rσ1 ) = w˜ τ (Rσ2 ) satisfies σ1 ≤ τ ≤ σ2 . If, in addition, |(Rσ2 \Rσ1 ) ∩ Nθ | = 1, then, since as noted in the proof of Lemma 7, changes in the slope of the cut function imply an increase in |Rτ ∩ Nθ |, τ is a breakpoint, and is the only breakpoint contained in (σ1 , σ2 ). Otherwise, if τ is not a breakpoint, then w˜ τ (Rσ1 ) > w˜ τ (Rτ ) and w˜ τ (Rσ2 ) > w˜ τ (Rτ ). See Figure 3. Thus, there is at least one breakpoint in each of (σ1 , τ) and (τ, σ2 ). Given Rσ1 ⊂ Rσ2 , it is simple algebra to find τ with w˜ τ (Rσ1 ) = w˜ τ (Rσ2 ). With τ in hand, Rτ can be computed with one maximum-flow computation in Hτ . Thus, via the divide-and-conquer strategy described above, Ogier used at most 2n − 3 maximum-flow computations in the appropriate time-expanded graphs to find all breakpoints τ and corresponding cuts Rτ in (θ − 1, θ]. Over all intervals, this implies O(nk) maximum-flow computations on graphs with O(nk) vertices, yielding an O(n4 k4 ) algorithm to compute W. Substituting the fastest strongly polynomial maximum-flow algorithm—the preflow-push algorithm of Goldberg and Tarjan [10], Ogier’s approach leads to an O(k3 mn2 log(kn2 /m)) algorithm to compute x∗ . 1.3.
Parametric Maximum Flows
To speed up the computation of W and the xτ , τ ∈ W, we generalize the parametric maximum-flow algorithm of Gallo, Grigoriadis, and Tarjan [9]. A parametric flow problem is defined by a network whose capacities are parameters of one variable. Given a sequence of parameter values σ1 < σ2 < · · · < σr , the goal is to compute a maximum flow and minimum cut for the network defined by each parameter. For general-capacity functions, this may require r separate maximum-flow computations. However, Gallo et al. [9] NETWORKS–2001
121
showed that, for a special class of capacity functions, these computations can be done in the same time as a maximum- flow computation using the preflow-push algorithm of Goldberg and Tarjan [10]. The special class that Gallo et al. considered is the class of networks on n vertices and m arcs such that the capacities of arcs leaving the source are nondecreasing functions of σ, the capacities of arcs entering the sink are nonincreasing functions of σ, and all other capacities are constant. This class of networks has the property that the minimum-cardinality, minimum cuts for each parameter are nested. They used this fact to design an efficient algorithm that they called a parametric preflow algorithm to compute a maximum flow and a minimum cut for each value of σ in an increasing sequence in the same asymptotic time as the fastest known strongly polynomial time algorithm to compute one maximum flow [10]. The parametric preflow algorithm is a parameterized version of the Goldberg and Tarjan [10] preflow-push algorithm for computing maximum flows. The preflowpush algorithm maintains, at all times, a feasible preflow and a valid labeling. A feasible preflow is a flow f satisfying arc capacity constraints P fij ≤ uij and relaxed node conservation constraints i∈V [fij −fji ] ≥ 0 for all P j ∈ V − {s}. The quantity i∈V [fij − fji ], when strictly positive, is called the excess at node j. A valid labeling l for preflow f is a function from the vertices to the nonnegative integers satisfying l(s) = n, l(d) = 0, l(i) ≥ 0, and l(i) ≤ l(j) + 1 for all (i, j) such that fij < uij . The preflow-push algorithm typically starts with the preflow defined by saturating all arcs leaving the source and the labeling l(i) = 0 for all i ≠ s. Throughout the course of the algorithm, labels of nodes may increase but never decrease. The algorithm terminates when no nodes besides the source or sink contain excess. The entire analysis of the algorithm depends on the fact that the labels can only increase and that, at all times, there is a path of residual arcs from any vertex to either the source or the sink, implying that no node has label greater than 2n. The bound on the run time of the algorithm is determined by bounding the number of times each node is relabeled. The analysis of this algorithm remains valid as long as it starts with any feasible preflow and corresponding valid labeling. There is nothing that requires the algorithm to start with the stated preflow and labeling. This fact was exploited by Gallo et al. in their generalization of this algorithm to solve the parametric maximum-flow problem. The parametric preflow algorithm starts with capacities determined by the smallest value of the parameter σ, that is, σ1 , and computes a maximum flow. Arc capacities are then increased to the next largest value of σ. If there are arcs from s to nodes with label less than l(s), the flow on these arcs is increased to meet the new capacities, increasing the excess at nodes adjacent to the
122
NETWORKS–2001
source. The flow on arcs entering the sink is decreased, if necessary, to meet the new capacities, increasing the excess at the nodes adjacent to the sink. This results in a new, feasible preflow for which the previous labeling is still valid, since all saturated arcs in the previous preflow remain saturated in the new preflow. Thus, the preflowpush algorithm will compute a maximum flow in this network with updated capacities. The same bound on the run time of the preflow-push algorithm is also a bound on the run time of the parametric preflow push algorithm: O(nm log(n2 /m)) time [9, 10] (assuming the number of values of σ is not more than m log(n2 /m)). This is because the bound on the run-time of the preflow-push algorithm depends only on the number of times a node is relabeled, node labels never decrease in either algorithm, and all node labels are ≤ 2n. The parametric preflow algorithm can be extended to work efficiently for more general classes of graphs than the special class described by Gallo et al. This has been shown for specific examples by Gusfield and Martel [11], Hochbaum and Hong [13], and McCormick [18], namely, all the algorithm requires is the following two characteristics: Lemma 9 (Gallo, Grigoriadis, Tarjan [9]). A parametric maximum-flow problem has a parametric preflow algorithm if 1. The minimum-cardinality, minimum cuts of the network are nested, and 2. After computing the maximum flow for one value of σ, and increasing σ, it must be possible to adjust the flow so that it is a feasible preflow and to modify the labeling so that it remains valid for the new residual graph, without decreasing any labels.
Once a network has a parametric preflow algorithm specified by the flow and label adjustment, Gallo et al. showed that this can be used as a subroutine in more sophisticated algorithms with the same asymptotic runtime behavior, as long as the capacities are linear functions of the parameter [9]. Let κ be the minimum-cut function of the network, that is, κ(σ) equals the value of the minimum cut with capacities determined by σ; it is the lower envelope of the cut-value functions. Since capacity functions are linear and by the nestedness property of minimum cuts, κ is determined by at most n − 1 cuts and has at most n − 2 breakpoints. Gallo et al. showed how the parametric preflow algorithm can be used to compute all breakpoints of κ in the same asymptotic time as one parametric preflow computation. Once the breakpoints are found, the parametric preflow algorithm can be invoked again to compute the maximum flows and minimum cuts corresponding to these breakpoints. Lemma 10 (Gallo, Grigoriadis, Tarjan [9]). If a parametric flow problem has a parametric preflow algorithm, and all the capacities are linear functions of the parameter, then all breakpoints of the minimum-cut function of
the graph can be computed in the same asymptotic time as the time required by the parametric preflow algorithm. This use of the parametric preflow algorithm is independent of the graph class. As long as the class satisfies both characteristics above, the parametric preflow algorithm can be used to compute the breakpoints efficiently. The details of how this is done are described in [9]. For this reason, to demonstrate an efficient algorithm to compute the breakpoints of κ, it suffices to show the existence of an efficient parametric preflow algorithm via establishing the two characteristics described above. In the next section, we show that a parametric network defined by the Hτ , θ − 1 < τ ≤ θ satisfies these characteristics. This will imply that it is possible to determine all breakpoints of κ in an interval (θ − 1, θ]. We will then show how to use such a breakpoint algorithm to compute all breakpoints in (0, k] in the same asymptotic time. 2. SOLVING UMFP In this section, we discuss the main contribution of this paper, which is a generalization of the parametric maximum-flow algorithm of Gallo, Grigoriadis, and Tarjan [9] to solve a subproblem of UMFP. Our generalization of [9] enables us to reduce the time needed to compute the set of breakpoints of the universally maximum flow W and the τ-maximum flows xτ , τ ∈ W that are necessary to compute the optimal flow x∗ as detailed in Section 1.2. We reduce the time needed by a factor of O(nk). 2.1. A Generalized Parametric Preflow Algorithm We consider a parametric flow problem based on the graphs Hτ . Instead of considering the graphs Hτ for τ ∈ (0, k] as separate graphs, we consider one graph H on the same vertex set but with parameterized capacities, so that the capacities of arcs of H at time τ equal the capacities of Hτ , that is, H(τ) = Hτ . More precisely, using Nθ to denote the θth copy of network N in H, an arc in Nθ has capacity 0 for 0 ≤ t ≤ θ − 1, capacity (t − θ + 1)ue (θ) for θ − 1 < t ≤ θ and capacity ue (θ) for θ < t ≤ k. An arc from Nθ−1 to Nθ has capacity 0 for 0 ≤ t ≤ θ and capacity uj (θ) for θ < t ≤ k. Note that the capacity functions in H are piecewise linear. Thus, we may not invoke Lemma 10 directly. However, the capacities are linear functions within unit time intervals (θ−1, θ] for θ ∈ {1, 2, . . . , k}. For this reason, we start by restricting our attention to a unit time interval (θ−1, θ] and consider the reverse network HR of H restricted to this time interval. (The reverse network is the network with the directions of all arcs reversed.) In the reverse network, d is the source. A minimum-
cardinality source side cut in a reverse network is a minimum-cardinality sink side cut in the original network. Thus, to show that we can compute breakpoints in this interval, it suffices to show that, given a sequence of values θ − 1 = τ0 < τ1 < τ2 < . . . < τr = θ, we can compute the minimum cuts and maximum flows for each value τh via a parametric preflow algorithm. By Lemma 4, we have already established that the cuts in this network are nested. Thus, by Lemma 9, it remains to show that we can adjust the preflow and the labeling. After computing the first maximum flow, we increase the current value from τh−1 to τh , thus increasing the capacity of arcs in Nθ by a multiplicative factor of δh h −τ0 := ττh−1 . To adjust the preflow, we then increase the −τ0 current flow x to x0 by multiplying the flow on each arc (iθ , jθ ) by δh . Lemma 11. If x is a preflow in HR (τl−1 ) with a valid labeling, then x0 obtained by multiplying the xi,j (θ) by δh is a preflow in HR (τh ) for which the same labeling is valid. Proof. Since there are no arcs with positive capacity that enter Nθ in HR , there is no node in Nθ that receives flow from outside Nθ . Hence, this adjustment of the preflow maintains nonnegative excess within Nθ and the adjusted preflow is feasible. Since this update keeps previously saturated arcs saturated and empty arcs empty, the current labeling is still valid. Lemma 11 establishes the requirements for the existence of an efficient parametric preflow algorithm. Thus, by Lemma 10, we may assume that we have an efficient algorithm to compute the breakpoints of W ∩ (θ − 1, θ]. 2.2.
Solving UMFP with Parametric Preflow
To find an optimal solution to UMFP, we define and solve three different types of parametric maximum-flow problems: one to compute the minimum cuts Rθ and θ-maximum flows xθ in Hθ for θ = 1, . . . , k; one to compute the minimum cuts Rτ and corresponding τmaximum flows xτ for all τ ∈ (θ − 1, θ] ∩ W, for each θ = 1 . . . , k; and one to compute x0 . 2.2.1. Computing Rθ and xθ for θ ∈ {1, . . . , k}. To compute the θ-maximum flows xθ , θ ∈ {1, . . . , k}, we construct a parametric flow problem on a modified Gk (= Hk ). Recall that si is the copy of the source in Ni , and di is the corresponding sink. We introduce a supersource sS with infinite capacity arcs (sS , si ) for each i ∈ {1, . . . , k} and a supersink dS with arcs (di , dS ) for each i ∈ {1, . . . , k} with a capacity function that is infinite when i ≤ θ and zero otherwise and call this parametric network Gˆ k . (See Fig. 4.) We then solve the parametric flow problem in Gˆ k to find the maximum flows corresponding to each parameter θ ∈ {1, . . . , k}.
NETWORKS–2001
123
FIG. 4. Gˆ k for k = 4.
Lemma 12. The maximum flow in Gˆ k (θ) is the θmaximum flow xθ , for θ ∈ {1, . . . , k}. Proof. The maximum flow in Gˆ k (θ) is a maximum flow in Gθ , which is xθ by Corollary 2. To solve this parametric flow problem, we reverse Gˆ k by reversing the direction of every arc in Gˆ k and apply the algorithm in [9] to this network. This parametric flow problem is of the form that is solved in [9]: All arcs that vary with θ are leaving the source, and the capacities are all increasing functions of θ. We denote the algorithm that solves parametric flow problems of this form by GGT. Thus, we have established that xθ and Rθ for all θ ∈ {1, . . . , k} can be computed in the same asymptotic time as one maximum-flow computation in the graph Gˆ k : O(k2 nm log(kn2 /m)) time. This part of the algorithm corresponds to Step 1 in Figure 5. 2.2.2. Computing Rτ and xτ for All τ ∈ W. To find the elements of W, and the corresponding τ-maximum flows within interval (θ − 1, θ], we start with Rθ−1 and Rθ in graphs H(θ − 1) and H(θ) and invoke the breakpoint algorithm implied by Lemmas 9, 10, and 11. We denote this algorithm by ModifiedGGT. The breakpoint algorithm is repeated for each unit interval of the form (θ − 1, θ], θ = 1, . . . , k in [0, k]. This could take k parametric maximum-flow computakn2 tions in HR , that is, k times O(k2 nm log m ) time. To speed this up, we use the fact that we have already computed Rθ for all θ ∈ {1, . . . , k} in Section 2.2.1., and that, by Lemma 4, Rθ−1 ⊆ Rτ ⊆ Rθ , for all θ − 1 ≤ τ ≤ θ. Thus, in searching for breakpoints of κ(τ) in the interval (θ−1, θ], we can restrict the computations performed for each unit interval to the subgraph of HR on the vertex set (Rθ − Rθ−1 ) ∪ {v1 , v2 }, where v1 is the vertex resulting from the contraction of all vertices in Rθ and v2 is the vertex resulting from the contraction of Rθ−1 . For any two values of θ, these vertex sets are disjoint. Thus, the total time required to compute W∩(θ−1, θ] for all values of θ is bounded by the time needed to run ModifiedGGT on a graph of size equal to the sum of the sizes of the kn2 components, which is, at most, O(k2 nm log m ). Along with the breakpoints τ ∈ W ∩ (θ − 1, θ], ModifiedGGT
124
NETWORKS–2001
also returns the Rτ implicitly as Rτ − Rθ−1 . Thus, at this point, we have determined W and the corresponding Rτ for τ ∈ W. It remains to compute the τ-maximum flows xτ for all τ ∈ W. We do not actually need to know the flow on every arc at every time step for each xτ . By Theorem 6, it τ (t) for t and τ such is sufficient to determine the values xij that qi (t) = qj (t) = τ. Since xτ is constant on unit intervals by definition (see Section 1.2), the definition of qj (t) τ (t) only in (2) implies that it is necessary to compute xij for t ≤ τ, where i(dte), j(dte) ∈ Rτ and i(dte), j(dte) ∈ Rτ+1 . By Corollary 3, xτ corresponds to a maximum flow in Hτ . For τ ∈ {τ0 , τ1 , . . . , τr }, the particular flow values that we require are on arcs contained in the subgraph of HR induced by the vertices in Rτh+1 − Rτh . To compute xτh restricted to Rτh+1 − R, we must also specify the flow on all arcs entering and leaving this interval, that is, all arcs with exactly one endpoint in Rτh+1 − Rτh . Let ∆(Rτ ) be the set of arcs in Hτ with exactly one endpoint in Rτ . The flow on arcs in ∆(Rτh ) is determined by definition of Rτh : All arcs entering Rτh must be saturated in xτh , and all arcs leaving Rτh must be empty in xτh . The xτh flow on arcs in ∆(Rτh+1 ) is similarly determined by employing one additional lemma of Ogier’s. The following corollary along with the above discussion then implies that we can compute xτh on HR restricted to Rτh+1 \Rτh . This part of the algorithm is summarized in Step 2 of Figure 5. Lemma 13 (Ogier [21]). For each τ ∈ (0, T), Rτ+ is a τ-minimum cut. Corollary 14. The flow xτh saturates all arcs entering Rτh+1 and equals zero on all arcs leaving Rτh+1 .
FIG. 5. The algorithm to compute x0 and the τ-maximum flows needed for constructing a universally maximum flow.
Proof. For consecutive breakpoints τh < τh+1 in W, Rτh + = Rτh+1 , so Rτh+1 is also a τh -minimum cut. Thus, the above observations for arcs entering and leaving Rτh also hold for arcs in ∆(Rτh+1 ). 2.2.3. Computing x0 . Recall from Section 1.2 the definition of x0 : the flow over time such that for all t, x0 (t) is a maximum flow for a static network N with capacities u(t) = u(dte), sources Ndte \Rt , and sink d. Since Ndte \Rt is decreasing in t, we can compute x0 on Nθ with an application of the parametric preflow algorithm [9]. Consider Nθ with a supersource sˆθ attached to each node of Nθ except dθ . The capacity of arc (ˆsθ , iθ ) is ∞ until time min{θ, qi ((θ−1)+)} and afterward 0. The reverse network satisfies conditions for the applicability of the parametric preflow algorithm (see Section 1.3). Thus, we can compute x0 with k + 1 calls to GGT on a network with n nodes and m arcs. This time is dominated by the time required to compute W and xτ for all τ ∈ W. This part of the algorithm is summarized in Step 3 of Figure 5. The following theorem summarizes the above discussions of the steps in computing x∗ : Theorem 15. A universally maximum flow with piecewise-constant capacities can be constructed in O(k2 mn log(kn2 /m)) time. Acknowledgments I am grateful to Éva Tardos and Kevin Wayne for providing helpful comments on early drafts of this paper and to an anonymous referee for useful comments. REFERENCES [1] E.J. Anderson and P. Nash, Linear programming in infinitedimensional spaces, Wiley, New York, 1987. [2] E.J. Anderson, P. Nash, and A.B. Philpott, A class of continuous network flow problems. Math Oper Res 7 (1982), 501–514. [3] E.J. Anderson and A.B. Philpott, A continuous-time network simplex algorithm, Networks 19 (1989), 395–425. [4] J.E. Aronson, A survey of dynamic network flows, Ann Oper Res 20 (1989), 1–66. [5] R.E. Burkard, K. Dlaska, and B. Klinz, The quickest flow problem, ZOR Methods Models Oper Res 37 (1993), 31– 58. [6] L. Fleischer, Faster algorithms for the quickest transshipment problem, SIAM J Option 12 (2001), 18–35 (extended abstract appeared in Proc SODA 1998). [7] L. Fleischer and É. Tardos, Efficient continuous-time dynamic network flow algorithms, Oper Res Lett 23 (1998), 71–80. [8] L.R. Ford and D.R. Fulkerson. Flows in networks, Princeton University Press, Princeton, NJ, 1962. [9] G. Gallo, M.D. Grigoriadis, and R.E. Tarjan, A fast parametric maximum flow algorithm and applications, SIAM J Comput 18 (1989), 30–55.
[10] A.V. Goldberg and R.E. Tarjan, A new approach to the maximum flow problem, J ACM 35 (1988), 921–940. [11] D. Gusfield and C. Martel, A fast algorithm for the generalized parametric minimum cut problem and applications, Algorithmica 7 (1992), 499-519. [12] B. Hajek and R.G. Ogier, Optimal dynamic routing in communication networks with continuous traffic, Networks 14 (1984), 457–487. [13] D.S. Hochbaum and S.P. Hong, About strongly polynomial time algorithms for quadratic optimization over submodular constraints, Math Program 69 (1995), 269–309. [14] B. Hoppe, Efficient dynamic network flow algorithms, PhD Thesis, Cornell University, June 1995, Department of Computer Science Technical Report TR95-1524. [15] B. Hoppe and É. Tardos, Polynomial time algorithms for some evacuation problems, Proc 5th Annual ACM-SIAM Symp on Discrete Algorithms, 1994, pp. 433–441. [16] B. Hoppe and É. Tardos, The quickest transshipment problem, Math Oper Res 25 (2000), 36–62 (extended abstract appeared in Proc SODA 1995). [17] B. Klinz and G.J. Woeginger, “Minimum cost dynamic flows: The series-parallel case,” Integer programming and combinatorial optimization, E. Balas and J. Clausen (Editors), Proc 4th Int IPCO Conf, Lecture Notes in Computer Science 920, Springer-Verlag, 1995, pp. 329–343. [18] S.T. McCormick, Fast algorithms for parametric scheduling come from extensions to parametric maximum flow, Oper Res 47 (1999), 744–756 (extended abstract appeared in Proc STOC 96). [19] E. Minieka, Maximal, lexicographic, and dynamic network flows, Oper Res 21 (1973), 517–527. [20] F.H. Moss and A. Segall, An optimal control approach to dynamic routing in networks, IEEE Trans Autom Control 27 (1982), 329–339. [21] R.G. Ogier, Minimum-delay routing in continuous-time dynamic networks with piecewise-constant capacities, Networks 18 (1988), 303–318. [22] J.B. Orlin, Maximum-throughput dynamic network flows, Math Program 27 (1983), 214–231. [23] J.B. Orlin, Minimum convex cost dynamic network flows, Math Oper Res 9 (1984), 190–207. [24] A.B. Philpott, Continuous-time flows in networks, MOR 15 (1990), 640–661. [25] A.B. Philpott and M. Craddock, An adaptive discretization method for a class of continuous network programs, Networks 26 (1995), 1–11. [26] W.B. Powell, P. Jaillet, and A. Odoni, “Stochastic and dynamic networks and routing,” Handbooks in operations research and management science: Networks, M.O. Ball, T.L. Magnanti, C.L. Monma, and G.L. Nemhauser (editors), Elsevier, 1995. [27] M.C. Pullan, An algorithm for a class of continuous linear programs, SIAM J Control Optim 31 (1993), 1558–1577. [28] M.C. Pullan, A study of general dynamic network programs with arc time-delays, SIAM J Optim 7 (1997), 889– 912. [29] G.I. Stassinopoulos and P. Konstantopoulos, Optimal congestion control in single destination networks, IEEE Trans Commun 33 (1985), 792–800. [30] W.L. Wilkinson, An algorithm for universal maximal dynamic flows in a network, Oper Res 19 (1971), 1602–1612.
NETWORKS–2001
125
Lisa K. Fleischer Graduate School of Industrial Administration, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, Pennsylvania 15213
A maximum flow over time generalizes standard maximum flow by introducing a time component. The object is to send as much flow from source to sink in T time units as possible, where capacities are interpreted as an upper bound on the rate of flow entering an arc. A related problem is the universally maximum flow, which is to send a flow from source to sink that maximizes the amount of flow arriving at the sink by time t simultaneously for all t ≤ T . We consider a further generalization of this problem that allows arc and node capacities to change over time. In particular, given a network with m arcs and n nodes, arc and node capacities that are piecewise constant functions of time with at most k breakpoints and a time bound T , we show how to compute a flow that maximizes the amount of flow reaching the sink in all time intervals (0,t ] simultaneously for all 0 < t ≤ T , in O (k 2mn log (kn 2/mm )) time. The best previous algorithm requires O (nk ) maximum-flow computations on a network with (m + n )k arcs and nk nodes. Our algorithm improves on this by a factor of O (nk ). © 2001 John Wiley & Sons, Inc. Keywords: dynamic flow; transshipment; parametric flow
1. INTRODUCTION In the 1960s, Ford and Fulkerson introduced flows over time to add time to the standard network flow model. Since then, flows over time have been used widely to model network-structured, decision-making problems over time: problems in electronic communication, production and distribution, economic planning, cash flow, job scheduling, and transportation. For examples, see the surveys of Aronson [4] and Powell et al. [26]. Flows over time have been previously referred to as dynamic network flows [4, 7, 12, 16, 17, 19, 20, 21, 22, 23, 26, 28, 30]. However, this term causes confusion about the problem being considered: “Dynamic” is more
Received November 1, 1998; accepted June 1, 2001 Correspondence to: L. K. Fleischer; e-mail:[email protected] Contract grant sponsor: AAUW; contract grant info: NSF; grant no. EIA-0049084 c 2001 John Wiley & Sons, Inc.
NETWORKS, Vol. 38(3), 115–125 2001
consistently used for a problem with input that changes online, over time. In many of the earlier papers, the network is fixed. While we allow for the network to change in this paper, all the changes are specified up front. In this paper and the earlier papers, it is the solution that changes over time. For this reason, we use the more descriptive, albeit slightly awkward, terminology of flows over time. A maximum flow over time generalizes standard maximum flow by the addition of time. A standard network consists of a set of nodes V with |V| = n and a set of arcs E with |E| = m, which is a subset of V × V. The capacity function from the arcs to the real numbers bounds the amount of flow allowed on each arc. In a flow-overtime network, the capacity function u limits the rate of flow into an arc. In addition, a flow-over-time network has a transit-time vector % associated with the arcs. The transit time of an arc is the amount of time it takes for flow to travel from one end of the arc to the other. Ford and Fulkerson [8] considered maximum flows over time: given a flow-over-time network with a specified source and sink, determine the maximum amount of flow that can be sent from source to sink in T time units. They show that the problem can be solved in polynomial time by using information obtained from a minimum-cost flow computation in a related network of comparable size. A generalization of this, a universally maximum flow, is a flow that simultaneously maximizes the amount of flow reaching a specified sink node d, sent from a specified source node s, by time t, for all 0 < t ≤ T. Wilkinson [30] and Minieka [19] showed that such a flow exists, and they both provided algorithms to solve this problem, but these do not run in polynomial time. There is no known polynomial time algorithm that solves the universally maximum-flow problem. It is not even known if the optimal solution is polynomial in the size of the input. Hoppe and Tardos [14, 15] presented an O(m−1 (m + n log n) log U) algorithm that computes a flow over time with the property that the quantity of
flow reaching the sink by time t is within (1 − ) of the maximum possible for all 0 ≤ t ≤ T. Another flow-over-time problem is the problem with multiple sources and sinks. A transshipment over time is a single-commodity flow that sends specified supplies located at source nodes to satisfy specified demands at sink nodes within a given time bound T. A universally quickest transshipment is a transshipment over time that satisfies as much demand as possible in every interval of time (0, t] for 0 < t ≤ T. If there is more than one sink, a universally quickest transshipment may not exist, so the universally quickest transshipment problem refers to a transshipment problem with one sink only. This paper considers the special case of flow problems over time when all transit times are zero—in other words, when flow is instantaneous. However, finite capacities still limit the rate of flow. In this setting, the universally maximum-flow problem is solved by sending flow at the rate of a maximum flow in the static network at every moment of time in (0, T]. (The static network is the interpretation of the flow-over-time network as a standard network. Transit times are ignored, and capacities are interpreted as bounds on the total flow allowed on an arc.) A universally quickest transshipment with zero transit times has a more complicated structure. This problem models the situation when the supplies and demands exceed network capacity, and it is necessary to send flow over time, in rounds. Unlike standard flow problems with multiple sources and sinks, this problem cannot be modeled as an equivalent s − t maximum flow over time. The reason such a transformation does not work here is because the arc capacities are upper bounds on the rate of flow per time period and do not correspond to the total amount of flow an arc may carry throughout an interval. Transshipment over time with zero transit times was discussed in [6, 12, 20, 29]. Hajek and Ogier [12] described the first polynomial time algorithm to find a universally quickest transshipment when all transit times are instantaneous. Their algorithm uses n maximum-flow computations on the underlying static network. In [6], this is improved upon by an algorithm that runs in the same asymptotic time as a preflow-push maximum-flow algorithm. All the above-mentioned problems may be generalized to allow storage of flow at the nodes. In this setting, there may be node-storage capacities that limit the total excess flow allowed at a node. Flow deficit is not allowed. In the above-mentioned problems, even when node excess is allowed, there always exists an optimal solution that does not use any node storage, and the algorithms mentioned above find such a solution. For a universally maximum flow with piecewise-constant capacity functions, the subject of this paper, this will not be the case: Node storage will be used in an optimal solution, and the amount of available storage will affect the solution.
116
NETWORKS–2001
Many other generalizations of this problem have been considered in the literature; however, few other generalizations are known to be polynomial-time solvable. One exception is an integer transshipment over time with general transit times, for which Hoppe and Tardos [16] gave the only polynomial-time algorithm known to solve this problem. Their algorithm repeatedly calls an oracle to minimize submodular functions. Orlin [23] provided polynomial- time algorithms for some infinite horizon problems, where the objective function is related to average throughput. Burkard et al. [5] found the minimum feasible time to send a specified flow from source to sink faster than incorporating the Ford–Fulkerson algorithm in a binary search framework. Most other work has focused on characterizing the structure of problems with more general capacity functions and determining when optimal solutions exist [1, 2, 3, 24, 25, 27, 28]. None of these problems are known to be polynomial-time solvable. Among the earlier work, Anderson et al. [2] considered the problem of finding a maximum flow in a network with zero transit times and with capacities and node storage that are much more general functions of time. They developed a duality theory for continuous network flows. We consider the generalization of the zero transit time, universally maximum-flow problem that allows both arc capacities and node-storage capacities to be piecewise-constant functions on [0, T] with at most k breakpoints. As with the general universally maximumflow problem, it is not clear that a universally maximum flow exists in these circumstances. Ogier [21] proved that it does and provided a polynomial-time algorithm that finds a universally maximum flow with at most nk breakpoints. These breakpoints can be computed with nk maximum-flow computations on a network with nk vertices and (m + n)k arcs. After the breakpoints are computed, additional maximum-flow computations on the same network are used to calculate the flow between two breakpoints. Thus, the total run time of Ogier’s algorithm is determined by the time to solve O(nk) maximum-flow problems on a network with nk vertices and (n + m)k arcs. As Ogier [21] demonstrated, this problem also generalizes the universally quickest transshipment problem. The main contribution of this paper is to (a) recognize that these problems can be solved by solving a parametric maximum-flow problem on a suitably defined graph and (b) to generalize the parametric maximum-flow algorithm of Gallo, Grigoriadis, and Tarjan [9] to fit the needs of this problem. The end result is that all the computations described in [21] can be performed in the same asymptotic time as one preflow-push maximumflow computation on a network with nk vertices and (n + m)k arcs. This improves the previous strongly polynomial run time by a factor of O(nk).
1.1.
The Model
The universally maximum-flow problem with piecewise-constant capacities (UMFP) is defined on a network with n nodes V and m arcs E, along with a nonnegative, left-continuous, piecewise-constant capacity function u : (E ∪ V) × (0, T] → Z+ ∪ {0}. When referring to the capacity of a particular arc e and a particular node i at a particular time t, we use the notation ue (t) and ui (t), respectively. The vector of all capacities at time t is simply denoted u(t), the vector of just arc capacities is denoted uE (t), and the vector of node capacities is denoted uV (t). For convenience of notation, we assume that the breakpoints of the capacity functions are exactly one unit apart. If this is not the case, we can scale time between each breakpoint to create an equivalent problem with this property; that is, if the time between two breakpoints is r, then we create an equivalent problem by changing this to 1 and multiplying all capacities in this interval by r. Any flow that solves this modified problem can be transformed back to a solution to the original problem by dividing the rate of flow sent in the new time interval by r and sending it over an interval that is r times longer. In such a manner, we scale T so that it equals k. Thus, for the remainder of the paper, we assume that our time bound is k, all k breakpoints are in the set {1, 2, . . . , k}, the domain of u is (E ∪ V) × (0, k], and u(t) = u(dte) for all 0 < t ≤ k. Let x be a flow over time. An arc e = (i, j) may be represented as ij in the context of a subscript: Here, xe and xij denote the same quantity, and similarly, for ue and uij . If (i, j) ∈ / E, then xij and uij are automatically set to zero. For A ⊂ V, we use the notation uA,A (t) to P denote u i∈A,j∈A ij (t) and the notation xA,A (t) to denote P [x i∈A,j∈A ij (t) − xji (t)]. We say that x is feasible if xe is a Lebesgue-measurable function on (0, k] for all e ∈ E, x obeys the arc capacity constraints xe (t) ≤ ue (t) for all e ∈ E, and all t ∈ (0, k], and the excess at each node, defined by Z τ xV,j (t)dt, pj (τ) := 0
obeys the node capacity constraint 0 ≤ pj (τ) ≤ uj (τ) for all j ∈ V − {s, d}, and all 0 < τ ≤ k. Denote the set of feasible flows over time by D. The value of flow x at time τ, denoted vτ (x), is the total flow reaching the sink by time τ, that is, Z τ xN,d (t)dt = pd (τ). vτ (x) :=
For each τ ∈ (0, k], we define a time-expanded graph Gτ on k copies of N, labeled N1 , . . . , Nk . Note that the vertex set is the same for all τ. Denote this vertex set by V∗ . Let iθ be the copy of vertex i in copy Nθ and eθ be the copy of arc e in Nθ . In Gτ , the arcs in Nθ for θ ∈ {1, . . . , dτe} are assigned capacities uE (θ); the arcs in Nθ for θ > dτe are given capacity 0; for θ = 1 to dτe−1, the holdover arcs leaving Nθ are given capacities equal to uV (θ); and for θ ≥ dτe, the holdover arcs leaving Nθ are assigned capacity 0. By this construction, we have that Gτ = Gσ if θ < τ ≤ σ ≤ θ + 1. Thus, there are exactly k unique graphs Gτ , one for each τ = θ ∈ {1, . . . , k}. The arcs with nonzero capacity in Gτ correspond to the unit intervals of time in which the corresponding arcs in the static network N may carry flow in any flow completing by time τ. Figure 1 depicts G2 and G2.5 for the simple network N on four nodes. The thin arcs correspond to arcs with 0 capacity. The thick arcs in Nθ correspond to arcs with capacity equal the capacity of the corresponding arc in N at time θ. The thick holdover arcs leaving Nθ have capacity equal to the capacity of the node they leave at time θ. The difference between G2 and G2.5 is that the arcs in N3 and the holdover arcs leaving N2 have nonzero capacity in G2.5 . This captures the fact that a flow completing by time τ = 2.5 may send some flow in N during the interval (2, 2.5], while a flow completing by time τ = 2 may not. 1.2.
Optimality Conditions
A τ-maximum flow is a feasible flow over time x that maximizes vτ (x). As in the static-flow setting, we can define a notion of time-cut that proves the optimality of a τ-maximum flow. A time-cut is a function C that maps (0, k] to subsets of vertices of the network N such that C(α) = C(β) for all α, β ∈ (θ − 1, θ], θ ∈ {1, . . . , k} and C(α) ∩ {s, d} = {d} for all α ∈ (0, k]. Note that this defines the cut as containing the sink. This is the complement of what is traditionally defined as a cut. Let C := {C|C is a time-cut}. Since the range of a timecut is a finite set, and over the domain (0, k], each timecut can change only k times (each is uniform over each interval (θ, θ + 1]), C is a finite set. The value wτ (C) of time-cut C at time τ is defined below: It is the sum of capacities of arcs with their head in C and tail in C for some interval of time, times the length of this interval (the first term), plus the sum of capacities of nodes at the last integer time they are in C (the second term).
0
In this paper, we will use a series of related timeexpanded graphs. A time-expanded network of N is a directed graph that contains a copy of N for each time unit and holdover arcs from a copy of a node at time θ to the copy of the same node at time θ + 1 for θ ∈ Z.
Z wτ (C) =
0
τ
uC(t)C(t) (t)dt +
dτe−1 X
X
θ=1 j∈C(θ+1)∩C(θ)
uj (θ). (1)
A time-cut has an alternate interpretation as a set of arcs entering a cut in an appropriately defined timeNETWORKS–2001
117
FIG. 1. Examples of time-expanded graphs Gτ and Hτ for τ = 2 and 2.5.
expanded network. This will be made more explicit in the corollaries and discussion following Theorem 1. A τ-minimum cut is a time-cut C that minimizes wτ (C). Theorem 1 states that the value of a τ-maximum flow equals the value of a τ-minimum cut. This is comparable to strong duality in standard network flows. The proof of even the weak duality case is a little more involved than is the corresponding static version and is omitted. A slightly weaker statement was proved by Anderson et al. [1], but for more general capacity functions. Philpott [24] extended this work to handle more general transit times. Theorem 1 (Ogier [21]). For all τ, minC∈C wτ (C) = maxx∈D vτ (x). The following corollary shows that for integral values of θ a θ-maximum flow corresponds to a maximum flow in Gθ . Corollary 2. If θ ∈ {1, 2, . . . , k}, then a maximum flow in the time-expanded network Gθ corresponds to a θmaximum flow xθ of the same value, where the flow on arc et in Gθ is precisely the flow on arc e in interval (t − 1, t], that is, for all t ∈ {1, . . . , θ}, the flow on et corresponds to xeθ (σ) for all σ ∈ (t − 1, t]. Proof. A maximum flow in Gθ gives rise to a minimum cut Rθ ⊂ V∗ that contains all copies of the sink node in Gθ and no copies of the source node. Letting C(σ) = C(dσe) := Ndσe ∩ Rθ for all σ ∈ (0, k], we have that C is a time-cut of a value equal to the maximum flow
118
NETWORKS–2001
in Gθ and, hence, equal to the value of the corresponding flow over time described in the corollary. We extend the series of time-expanded graphs Gτ to a more continuous version Hτ , so that we may extend Corollary 2 to apply more generally to τ ∈ (0, k]. For each real 0 < τ ≤ k, we define Hτ to be a graph on the same vertex set and arc set as Gdτe , but with a slightly different capacity function. In Hτ , the capacities of the arcs in Ndτe are multiplied by τ − bτc, which is always less than one. The capacities of all other arcs equal the capacities of the arcs in Gτ . Denote the capacity function of Hτ by w˜ τ . Note that Hθ = Gθ for all θ ∈ {1, . . . , k}. Examples of Hτ and Gτ are depicted in Figure 1 for a simple four-node network N with k = 4. Note that the capacity of arcs in N3 in H2.5 is half of the capacity of the corresponding arcs in G2.5 . This modification in H2.5 models the fact that a flow that completes by time τ = 2.5 uses arcs in N only for half of the unit interval (2, 3]. The value of a time-cut at time τ corresponds to the capacity of the arcs entering the sink side of a cut in Hτ . We make this explicit as follows: Since the vertex set is the same for all Hτ and equals V∗ , the vertex set of Gτ , we also use V∗ to denote the vertex set of Hτ . A cut in Hτ is a set R ⊂ V∗ such that R contains all copies of and no copies of the source, that Skthe sink nodeS k is, R ∩ θ=1 {sθ , dθ } = θ=1 {dθ }. ∗ For a set R ⊂ V , define w˜ τ (R) to be the sum of capacities of arcs in Hτ entering R. We can equate each time-cut C ∈ C with a cut R ⊂ V∗ of Hτ such that wτ (C) = w˜ τ (R) by setting R such that R ∩ Nθ = C(θ) for θ ∈ {1, . . . , k}. Let Dτ be the set of feasible flows over time that are obtained from flows in Hτ by setting the
rate of flow in arc e in interval (θ − 1, min{τ, θ}] to be the flow on arc eθ in Hτ . The following corollary follows easily from the above discussions: Corollary 3.
maxx∈Dτ vτ (x) = minR⊂V∗ w˜ τ (R).
By the submodularity of cut functions in a graph, the nonempty intersection of minimum cuts is a minimum cut. Thus, there is a minimum cut in Hτ with the smallest sink-side. Define Rτ to be the minimum cut such that Rτ ⊆ R for all R such that w˜ τ (R) = w˜ τ (Rτ ). The following lemma describes properties of Rτ as a function of τ:
Lemma 4 (Ogier [21]). For all 0 ≤ τ < σ ≤ k, Rτ ⊆ Rσ , and Rτ is a left continuous function of τ. We give some intuitive explanation for the behavior of Rτ between breakpoints of u. For τ close to bτc, the arcs in Ndτe have very little capacity in Hτ . Thus, a minimum cut in Hτ may have arcs entering the cut in Ndτe with large u(dτe) value to avoid holdover arcs in Hτ from Nbτc to Ndτe with large capacity. One way to avoid these holdover arcs is to exclude all nodes of Ndτe except ddτe . As τ increases, the contribution of the capacities of arcs in Ndτe toward the value of the
FIG. 2. (Top left) The time expanded graph G4 for k = 4 displaying all input capacities. (Second row, left) The graph H3.5 with capacities at half-value in N4 . (Top right) Cuts Rτ for τ = 1, 2, 3, 3.5, 4, demonstrating the behavior described in Corollary 4 and the subsequent discussion. (Second row, right) A τ-maximum flow xτ for τ = 3.5. Notice that this flow saturates the cut R3.5 in H3.5 . (Bottom left) the values of q. (Bottom right) the flows x0 and x∗ .
NETWORKS–2001
119
cut increases, and, thus, the minimum cut may shift to a smaller capacity cut of Ndτe , thus including more nodes of Ndτe and possibly increasing the total capacity of holdover arcs contributing to the value of the cut. Exactly this situation is depicted in Figure 2. Figure 2 shows Rτ for τ = 1, 2, 3, 3.5, 4, for the same original network N described in Figure 1. Here and elsewhere, Rτ+ indicates the limit of Rσ as σ ↓ τ. The figure also shows a flow in H3.5 . Note that this flow saturates the cut R3.5 and is thus a 3.5-maximum flow. Let W be the set of breakpoints of Rτ as a function of τ, that is, W is the set of values τ such that Rτ ≠ Rτ+ . Ogier’s main contribution is to prove that there is an optimal solution x∗ to UMFP that is uniform on the intervals between successive breakpoints in W and that such an optimal flow can be defined piecemeal by special τ-maximum flows [21]. This is a theorem of Ogier reprinted here as Theorem 6. Part of Theorem 6 can be explained by complementary slackness conditions. We first provide this explanation and then give the precise technical statement. We start with two easy observations: (1) For σ ≤ θ / Rσ . (2) An optimal flow x∗ for UMFP is a − 1, iθ ∈ τ-maximum flow for all τ ∈ (0, k] and thus saturates Rτ in Hτ for all τ. We now consider the three cases (A–C) that determine the flow on arc (i, j) in interval (θ − 1, θ]. / Rσ and jθ ∈ Case A. If there is a cut Rσ such that iθ ∈ Rσ , then arc (iθ , jθ ) is saturated in a σ-maximum flow and thus must be saturated in an optimal flow up until time ∗ (t) = uij (t) for t ∈ (θ − 1, min{θ, σ}]. A σ. Precisely, xij / Rσ and iθ ∈ Rσ symmetric argument shows that if jθ ∈ ∗ then xij (t) = 0 for t ∈ (θ − 1, min{θ, σ}]. Case B. If σ ≤ θ is the earliest time that iθ ∈ Rσ+ , then we know from the first observation that σ > θ − 1. Lemma 5 below implies that for all σ < τ ≤ θ there is a flow-saturating arc entering Rτ with excess pi (t) = 0 for all t ∈ (σ, θ]. Since the optimal flow x∗ must saturate arcs entering Rτ up until time τ, this implies that in an optimal flow x∗ we need not use any node capacity at i in (σ, θ]. Thus, between breakpoints σ ≤ σ1 < σ2 ≤ θ, the flow on all arcs in Nθ ∩ Rσ2 for t ∈ (σ1 , σ2 ] may be determined by the maximum flow in Nθ with sources Nθ \Rσ2 and sink dθ .
Lemma 5. For i ∈ V − {s, d}, if σ is the earliest time that iθ ∈ Rσ+ , then for all σ < τ ≤ θ, there is a flow in Dτ saturating arcs entering Rτ with excess pi (t) = 0 for all t ∈ (σ, θ]. Proof. Key to the proof of this lemma is that all flows in Dτ for τ ≤ θ correspond to flows in Hθ , which contains no holdover arcs for nodes in Nθ . The value of the cut function at Rσ and Rσ+ for τ ∈ (θ − 1, θ] is determined by arcs in Nθ which are parameterized by τ − θ + 1 and constant capacity arcs that are either
120
NETWORKS–2001
holdover arcs for times ψ ≤ θ − 1 or arcs in Nψ for ψ ≤ θ − 1. Thus, we can express these two values as (τ − θ + 1)aσ + bσ and (τ − θ + 1)aσ+ + bσ+ . For τ = σ, these two values are equal. By definition of Rσ and Rσ +, we have that aσ > aσ+ and bσ < bσ+ . Thus, 0 < bσ+ − bσ = (σ − θ + 1)(aσ − aσ+ ). In words, all flow that enters Rσ+ − Rσ along constant capacity arcs either leaves Rσ+ − Rσ along constant capacity arcs or drains from Rσ+ − Rσ by time σ. Since no positive, constant capacity arc leaving Rσ+ − Rσ is a holdover arc for any time after θ−1, this implies that no node capacity is used at any node in Nθ ∩ (Rσ+ − Rσ ) between time σ and time τ, where τ is the next breakpoint. In the interval (τ, θ], the only additional flow entering this set comes from the parametrized arcs in Nθ and thus is bounded from above by aσ . Since this flow saturates Rθ and drains by time θ, it drains in Nθ without requiring node storage. Case C. We are left with the case that σ > θ−1 is the last value such that both iθ and jθ are not in Rσ and for all τ > σ both iθ and jθ are in Rτ . In this case, since the flow crossing Rσ saturates all arcs crossing Rσ , the flow on all arcs on the source side of the cut should be enough to saturate Rσ . Thus, if iθ and jθ are not separated by some σ ∗ (t) = xij (t) for t ∈ (θ − 1, min{σ, θ}]. time-cut, then xij With this motivation in mind, we follow [21] and define xτ to be a τ-maximum flow that is constant on each interval (θ−1, θ] for θ ∈ {1, 2, . . . , bτc} and also constant on interval (bτc, τ]. In addition, let the excess function pτ of xτ satisfy pjτ (τ) = 0 for all j ∈ V − {s, d}. That such a flow exists is implied by Corollary 3: Compute a maximum flow in Hτ . This flow saturates the arcs entering Rτ and, hence, has the value of a τ-maximum flow. It does not use any node storage arc leaving N(dτe) and, hence, completes by time τ. This is depicted in Figure 2 for τ = 3.5. We also define x0 to be the flow over time such that for all t, f = x0 (t) is a maximum flow for static network N with capacities u(t) = u(dte), sources Ndte \Rt , and sink d. Figure 2 describes x0 for an example. For each node j of the network N, define max{τ|τ ≥ t, jdte ∈ Rτ }, if this set is nonempty, (2) qj (t) := 0, otherwise. In words, qj (t) is the largest value of τ for which there is no path from jdte to ddτe in the residual network of a τ-maximum flow. We use qj (t) to determine which of the above cases A–C applies. If qi (t) > qj (t): This implies that if idte has a path in the residual flow graph of xσ to ddσe then jdte has a residual path to ddσe and there is some τ for which idte does not have a residual path to ddτe when jdte does. For dte − 1 < τ, this is Case A for the interval (dte − 1, min{τ, dte}]. This is Case B for the interval (τ, dte].
Theorem 6 (Ogier [21]). Flow x∗ , defined by τ xij (t), if qi (t) = qj (t) = τ [ Case B if τ = 0; else Case C ] ∗ xij (t) = [ Case A ] uij (t), if qi (t) > qj (t) [ Case A ] 0, if qi (t) < qj (t) is a τ-maximum flow for all τ ∈ (0, T]. Notice that the only τ-maximum flows used to compute x∗ are for τ ∈ W, the set of breakpoints. Thus, we may compute a universally maximum flow by first computing W, the set of breakpoints, then computing x0 , Rτ , and xτ for all τ ∈ W, determining q, and then constructing x∗ as given by Theorem 6. For the example in Figure 2, the breakpoints are W = {1, 2, 3, 3.5, 4}. The corresponding flows, time-cuts, values of q, and final flow x∗ are all depicted in the figure. Once W, the set of all breakpoints of Rτ , along with the corresponding cuts Rτ , τ ∈ W, are known, then x0 can be computed with |W| maximum flows in a network the same size as N, and x∗ can be computed with an additional |W| maximum flows, one in each Hτ , for all τ ∈ W. Ogier described a method to compute W in a piecewise fashion, computing W ∩ (θ − 1, θ] for θ = 1, . . . , k. Define the minimum-cut function κ by setting κ(τ) equal to the value of the minimum cut in Hτ , that is, κ(τ) = w˜ τ (Rτ ). Figure 3 describes the relation between the Rτ and κ. Lemma 7. Within the domain (θ − 1, θ] for θ ∈ {1, . . . , k}, κ is concave and piecewise linear with at most n − 2 breakpoints. Proof. For τ ∈ (θ − 1, θ], θ ∈ {1, . . . , k}, κ(τ) is the lower envelope (the minimum function) of {w˜ τ (Rσ )}σ , a set of linear functions of τ. This implies that κ is concave. Since Rτ is increasing in τ by Lemma 4, changes in the slope of κ(τ) as τ increases imply a change in the cut Rτ defining κ and thus an increase in |Rτ ∩Nθ |. This / Rτ for all τ together with the fact that dθ ∈ Rτ and sθ ∈ implies that κ(τ) has at most n − 2 breakpoints in a unit interval. Ogier [21] computed each subset W ∩ (θ − 1, θ] by finding τ such that w˜ τ (R(θ−1)+ ) = w˜ τ (Rθ ), computing Rτ , and recursing on the subintervals (θ−1, τ] and (τ, θ] until each interval does not properly contain elements of W. The following lemma justifies this approach: Lemma 8. If τ is such that w˜ τ (Rσ1 ) = w˜ τ (Rσ2 ) for σ1 < σ2 , then σ1 ≤ τ ≤ σ2 and either τ is a breakpoint of κ or both (σ1 , τ) and (τ, σ2 ) contain breakpoints of κ. Proof. For a fixed cut Rσ and for τ in the interval between two breakpoints of u, w˜ τ (Rσ ) is a linear function
FIG. 3. The relationship of κ to w˜ τ and Rτ . Here, τ0 is such that w˜ τ0 (Rτ1 ) = w˜ τ0 (Rτ4 ).
with positive slope, that is, it has form τaσ + bσ . Since Rσ is the minimum cut when τ = σ, the slope of w˜ τ (Rσ ) as a function of τ decreases as σ increases. Thus, τ such that w˜ τ (Rσ1 ) = w˜ τ (Rσ2 ) satisfies σ1 ≤ τ ≤ σ2 . If, in addition, |(Rσ2 \Rσ1 ) ∩ Nθ | = 1, then, since as noted in the proof of Lemma 7, changes in the slope of the cut function imply an increase in |Rτ ∩ Nθ |, τ is a breakpoint, and is the only breakpoint contained in (σ1 , σ2 ). Otherwise, if τ is not a breakpoint, then w˜ τ (Rσ1 ) > w˜ τ (Rτ ) and w˜ τ (Rσ2 ) > w˜ τ (Rτ ). See Figure 3. Thus, there is at least one breakpoint in each of (σ1 , τ) and (τ, σ2 ). Given Rσ1 ⊂ Rσ2 , it is simple algebra to find τ with w˜ τ (Rσ1 ) = w˜ τ (Rσ2 ). With τ in hand, Rτ can be computed with one maximum-flow computation in Hτ . Thus, via the divide-and-conquer strategy described above, Ogier used at most 2n − 3 maximum-flow computations in the appropriate time-expanded graphs to find all breakpoints τ and corresponding cuts Rτ in (θ − 1, θ]. Over all intervals, this implies O(nk) maximum-flow computations on graphs with O(nk) vertices, yielding an O(n4 k4 ) algorithm to compute W. Substituting the fastest strongly polynomial maximum-flow algorithm—the preflow-push algorithm of Goldberg and Tarjan [10], Ogier’s approach leads to an O(k3 mn2 log(kn2 /m)) algorithm to compute x∗ . 1.3.
Parametric Maximum Flows
To speed up the computation of W and the xτ , τ ∈ W, we generalize the parametric maximum-flow algorithm of Gallo, Grigoriadis, and Tarjan [9]. A parametric flow problem is defined by a network whose capacities are parameters of one variable. Given a sequence of parameter values σ1 < σ2 < · · · < σr , the goal is to compute a maximum flow and minimum cut for the network defined by each parameter. For general-capacity functions, this may require r separate maximum-flow computations. However, Gallo et al. [9] NETWORKS–2001
121
showed that, for a special class of capacity functions, these computations can be done in the same time as a maximum- flow computation using the preflow-push algorithm of Goldberg and Tarjan [10]. The special class that Gallo et al. considered is the class of networks on n vertices and m arcs such that the capacities of arcs leaving the source are nondecreasing functions of σ, the capacities of arcs entering the sink are nonincreasing functions of σ, and all other capacities are constant. This class of networks has the property that the minimum-cardinality, minimum cuts for each parameter are nested. They used this fact to design an efficient algorithm that they called a parametric preflow algorithm to compute a maximum flow and a minimum cut for each value of σ in an increasing sequence in the same asymptotic time as the fastest known strongly polynomial time algorithm to compute one maximum flow [10]. The parametric preflow algorithm is a parameterized version of the Goldberg and Tarjan [10] preflow-push algorithm for computing maximum flows. The preflowpush algorithm maintains, at all times, a feasible preflow and a valid labeling. A feasible preflow is a flow f satisfying arc capacity constraints P fij ≤ uij and relaxed node conservation constraints i∈V [fij −fji ] ≥ 0 for all P j ∈ V − {s}. The quantity i∈V [fij − fji ], when strictly positive, is called the excess at node j. A valid labeling l for preflow f is a function from the vertices to the nonnegative integers satisfying l(s) = n, l(d) = 0, l(i) ≥ 0, and l(i) ≤ l(j) + 1 for all (i, j) such that fij < uij . The preflow-push algorithm typically starts with the preflow defined by saturating all arcs leaving the source and the labeling l(i) = 0 for all i ≠ s. Throughout the course of the algorithm, labels of nodes may increase but never decrease. The algorithm terminates when no nodes besides the source or sink contain excess. The entire analysis of the algorithm depends on the fact that the labels can only increase and that, at all times, there is a path of residual arcs from any vertex to either the source or the sink, implying that no node has label greater than 2n. The bound on the run time of the algorithm is determined by bounding the number of times each node is relabeled. The analysis of this algorithm remains valid as long as it starts with any feasible preflow and corresponding valid labeling. There is nothing that requires the algorithm to start with the stated preflow and labeling. This fact was exploited by Gallo et al. in their generalization of this algorithm to solve the parametric maximum-flow problem. The parametric preflow algorithm starts with capacities determined by the smallest value of the parameter σ, that is, σ1 , and computes a maximum flow. Arc capacities are then increased to the next largest value of σ. If there are arcs from s to nodes with label less than l(s), the flow on these arcs is increased to meet the new capacities, increasing the excess at nodes adjacent to the
122
NETWORKS–2001
source. The flow on arcs entering the sink is decreased, if necessary, to meet the new capacities, increasing the excess at the nodes adjacent to the sink. This results in a new, feasible preflow for which the previous labeling is still valid, since all saturated arcs in the previous preflow remain saturated in the new preflow. Thus, the preflowpush algorithm will compute a maximum flow in this network with updated capacities. The same bound on the run time of the preflow-push algorithm is also a bound on the run time of the parametric preflow push algorithm: O(nm log(n2 /m)) time [9, 10] (assuming the number of values of σ is not more than m log(n2 /m)). This is because the bound on the run-time of the preflow-push algorithm depends only on the number of times a node is relabeled, node labels never decrease in either algorithm, and all node labels are ≤ 2n. The parametric preflow algorithm can be extended to work efficiently for more general classes of graphs than the special class described by Gallo et al. This has been shown for specific examples by Gusfield and Martel [11], Hochbaum and Hong [13], and McCormick [18], namely, all the algorithm requires is the following two characteristics: Lemma 9 (Gallo, Grigoriadis, Tarjan [9]). A parametric maximum-flow problem has a parametric preflow algorithm if 1. The minimum-cardinality, minimum cuts of the network are nested, and 2. After computing the maximum flow for one value of σ, and increasing σ, it must be possible to adjust the flow so that it is a feasible preflow and to modify the labeling so that it remains valid for the new residual graph, without decreasing any labels.
Once a network has a parametric preflow algorithm specified by the flow and label adjustment, Gallo et al. showed that this can be used as a subroutine in more sophisticated algorithms with the same asymptotic runtime behavior, as long as the capacities are linear functions of the parameter [9]. Let κ be the minimum-cut function of the network, that is, κ(σ) equals the value of the minimum cut with capacities determined by σ; it is the lower envelope of the cut-value functions. Since capacity functions are linear and by the nestedness property of minimum cuts, κ is determined by at most n − 1 cuts and has at most n − 2 breakpoints. Gallo et al. showed how the parametric preflow algorithm can be used to compute all breakpoints of κ in the same asymptotic time as one parametric preflow computation. Once the breakpoints are found, the parametric preflow algorithm can be invoked again to compute the maximum flows and minimum cuts corresponding to these breakpoints. Lemma 10 (Gallo, Grigoriadis, Tarjan [9]). If a parametric flow problem has a parametric preflow algorithm, and all the capacities are linear functions of the parameter, then all breakpoints of the minimum-cut function of
the graph can be computed in the same asymptotic time as the time required by the parametric preflow algorithm. This use of the parametric preflow algorithm is independent of the graph class. As long as the class satisfies both characteristics above, the parametric preflow algorithm can be used to compute the breakpoints efficiently. The details of how this is done are described in [9]. For this reason, to demonstrate an efficient algorithm to compute the breakpoints of κ, it suffices to show the existence of an efficient parametric preflow algorithm via establishing the two characteristics described above. In the next section, we show that a parametric network defined by the Hτ , θ − 1 < τ ≤ θ satisfies these characteristics. This will imply that it is possible to determine all breakpoints of κ in an interval (θ − 1, θ]. We will then show how to use such a breakpoint algorithm to compute all breakpoints in (0, k] in the same asymptotic time. 2. SOLVING UMFP In this section, we discuss the main contribution of this paper, which is a generalization of the parametric maximum-flow algorithm of Gallo, Grigoriadis, and Tarjan [9] to solve a subproblem of UMFP. Our generalization of [9] enables us to reduce the time needed to compute the set of breakpoints of the universally maximum flow W and the τ-maximum flows xτ , τ ∈ W that are necessary to compute the optimal flow x∗ as detailed in Section 1.2. We reduce the time needed by a factor of O(nk). 2.1. A Generalized Parametric Preflow Algorithm We consider a parametric flow problem based on the graphs Hτ . Instead of considering the graphs Hτ for τ ∈ (0, k] as separate graphs, we consider one graph H on the same vertex set but with parameterized capacities, so that the capacities of arcs of H at time τ equal the capacities of Hτ , that is, H(τ) = Hτ . More precisely, using Nθ to denote the θth copy of network N in H, an arc in Nθ has capacity 0 for 0 ≤ t ≤ θ − 1, capacity (t − θ + 1)ue (θ) for θ − 1 < t ≤ θ and capacity ue (θ) for θ < t ≤ k. An arc from Nθ−1 to Nθ has capacity 0 for 0 ≤ t ≤ θ and capacity uj (θ) for θ < t ≤ k. Note that the capacity functions in H are piecewise linear. Thus, we may not invoke Lemma 10 directly. However, the capacities are linear functions within unit time intervals (θ−1, θ] for θ ∈ {1, 2, . . . , k}. For this reason, we start by restricting our attention to a unit time interval (θ−1, θ] and consider the reverse network HR of H restricted to this time interval. (The reverse network is the network with the directions of all arcs reversed.) In the reverse network, d is the source. A minimum-
cardinality source side cut in a reverse network is a minimum-cardinality sink side cut in the original network. Thus, to show that we can compute breakpoints in this interval, it suffices to show that, given a sequence of values θ − 1 = τ0 < τ1 < τ2 < . . . < τr = θ, we can compute the minimum cuts and maximum flows for each value τh via a parametric preflow algorithm. By Lemma 4, we have already established that the cuts in this network are nested. Thus, by Lemma 9, it remains to show that we can adjust the preflow and the labeling. After computing the first maximum flow, we increase the current value from τh−1 to τh , thus increasing the capacity of arcs in Nθ by a multiplicative factor of δh h −τ0 := ττh−1 . To adjust the preflow, we then increase the −τ0 current flow x to x0 by multiplying the flow on each arc (iθ , jθ ) by δh . Lemma 11. If x is a preflow in HR (τl−1 ) with a valid labeling, then x0 obtained by multiplying the xi,j (θ) by δh is a preflow in HR (τh ) for which the same labeling is valid. Proof. Since there are no arcs with positive capacity that enter Nθ in HR , there is no node in Nθ that receives flow from outside Nθ . Hence, this adjustment of the preflow maintains nonnegative excess within Nθ and the adjusted preflow is feasible. Since this update keeps previously saturated arcs saturated and empty arcs empty, the current labeling is still valid. Lemma 11 establishes the requirements for the existence of an efficient parametric preflow algorithm. Thus, by Lemma 10, we may assume that we have an efficient algorithm to compute the breakpoints of W ∩ (θ − 1, θ]. 2.2.
Solving UMFP with Parametric Preflow
To find an optimal solution to UMFP, we define and solve three different types of parametric maximum-flow problems: one to compute the minimum cuts Rθ and θ-maximum flows xθ in Hθ for θ = 1, . . . , k; one to compute the minimum cuts Rτ and corresponding τmaximum flows xτ for all τ ∈ (θ − 1, θ] ∩ W, for each θ = 1 . . . , k; and one to compute x0 . 2.2.1. Computing Rθ and xθ for θ ∈ {1, . . . , k}. To compute the θ-maximum flows xθ , θ ∈ {1, . . . , k}, we construct a parametric flow problem on a modified Gk (= Hk ). Recall that si is the copy of the source in Ni , and di is the corresponding sink. We introduce a supersource sS with infinite capacity arcs (sS , si ) for each i ∈ {1, . . . , k} and a supersink dS with arcs (di , dS ) for each i ∈ {1, . . . , k} with a capacity function that is infinite when i ≤ θ and zero otherwise and call this parametric network Gˆ k . (See Fig. 4.) We then solve the parametric flow problem in Gˆ k to find the maximum flows corresponding to each parameter θ ∈ {1, . . . , k}.
NETWORKS–2001
123
FIG. 4. Gˆ k for k = 4.
Lemma 12. The maximum flow in Gˆ k (θ) is the θmaximum flow xθ , for θ ∈ {1, . . . , k}. Proof. The maximum flow in Gˆ k (θ) is a maximum flow in Gθ , which is xθ by Corollary 2. To solve this parametric flow problem, we reverse Gˆ k by reversing the direction of every arc in Gˆ k and apply the algorithm in [9] to this network. This parametric flow problem is of the form that is solved in [9]: All arcs that vary with θ are leaving the source, and the capacities are all increasing functions of θ. We denote the algorithm that solves parametric flow problems of this form by GGT. Thus, we have established that xθ and Rθ for all θ ∈ {1, . . . , k} can be computed in the same asymptotic time as one maximum-flow computation in the graph Gˆ k : O(k2 nm log(kn2 /m)) time. This part of the algorithm corresponds to Step 1 in Figure 5. 2.2.2. Computing Rτ and xτ for All τ ∈ W. To find the elements of W, and the corresponding τ-maximum flows within interval (θ − 1, θ], we start with Rθ−1 and Rθ in graphs H(θ − 1) and H(θ) and invoke the breakpoint algorithm implied by Lemmas 9, 10, and 11. We denote this algorithm by ModifiedGGT. The breakpoint algorithm is repeated for each unit interval of the form (θ − 1, θ], θ = 1, . . . , k in [0, k]. This could take k parametric maximum-flow computakn2 tions in HR , that is, k times O(k2 nm log m ) time. To speed this up, we use the fact that we have already computed Rθ for all θ ∈ {1, . . . , k} in Section 2.2.1., and that, by Lemma 4, Rθ−1 ⊆ Rτ ⊆ Rθ , for all θ − 1 ≤ τ ≤ θ. Thus, in searching for breakpoints of κ(τ) in the interval (θ−1, θ], we can restrict the computations performed for each unit interval to the subgraph of HR on the vertex set (Rθ − Rθ−1 ) ∪ {v1 , v2 }, where v1 is the vertex resulting from the contraction of all vertices in Rθ and v2 is the vertex resulting from the contraction of Rθ−1 . For any two values of θ, these vertex sets are disjoint. Thus, the total time required to compute W∩(θ−1, θ] for all values of θ is bounded by the time needed to run ModifiedGGT on a graph of size equal to the sum of the sizes of the kn2 components, which is, at most, O(k2 nm log m ). Along with the breakpoints τ ∈ W ∩ (θ − 1, θ], ModifiedGGT
124
NETWORKS–2001
also returns the Rτ implicitly as Rτ − Rθ−1 . Thus, at this point, we have determined W and the corresponding Rτ for τ ∈ W. It remains to compute the τ-maximum flows xτ for all τ ∈ W. We do not actually need to know the flow on every arc at every time step for each xτ . By Theorem 6, it τ (t) for t and τ such is sufficient to determine the values xij that qi (t) = qj (t) = τ. Since xτ is constant on unit intervals by definition (see Section 1.2), the definition of qj (t) τ (t) only in (2) implies that it is necessary to compute xij for t ≤ τ, where i(dte), j(dte) ∈ Rτ and i(dte), j(dte) ∈ Rτ+1 . By Corollary 3, xτ corresponds to a maximum flow in Hτ . For τ ∈ {τ0 , τ1 , . . . , τr }, the particular flow values that we require are on arcs contained in the subgraph of HR induced by the vertices in Rτh+1 − Rτh . To compute xτh restricted to Rτh+1 − R, we must also specify the flow on all arcs entering and leaving this interval, that is, all arcs with exactly one endpoint in Rτh+1 − Rτh . Let ∆(Rτ ) be the set of arcs in Hτ with exactly one endpoint in Rτ . The flow on arcs in ∆(Rτh ) is determined by definition of Rτh : All arcs entering Rτh must be saturated in xτh , and all arcs leaving Rτh must be empty in xτh . The xτh flow on arcs in ∆(Rτh+1 ) is similarly determined by employing one additional lemma of Ogier’s. The following corollary along with the above discussion then implies that we can compute xτh on HR restricted to Rτh+1 \Rτh . This part of the algorithm is summarized in Step 2 of Figure 5. Lemma 13 (Ogier [21]). For each τ ∈ (0, T), Rτ+ is a τ-minimum cut. Corollary 14. The flow xτh saturates all arcs entering Rτh+1 and equals zero on all arcs leaving Rτh+1 .
FIG. 5. The algorithm to compute x0 and the τ-maximum flows needed for constructing a universally maximum flow.
Proof. For consecutive breakpoints τh < τh+1 in W, Rτh + = Rτh+1 , so Rτh+1 is also a τh -minimum cut. Thus, the above observations for arcs entering and leaving Rτh also hold for arcs in ∆(Rτh+1 ). 2.2.3. Computing x0 . Recall from Section 1.2 the definition of x0 : the flow over time such that for all t, x0 (t) is a maximum flow for a static network N with capacities u(t) = u(dte), sources Ndte \Rt , and sink d. Since Ndte \Rt is decreasing in t, we can compute x0 on Nθ with an application of the parametric preflow algorithm [9]. Consider Nθ with a supersource sˆθ attached to each node of Nθ except dθ . The capacity of arc (ˆsθ , iθ ) is ∞ until time min{θ, qi ((θ−1)+)} and afterward 0. The reverse network satisfies conditions for the applicability of the parametric preflow algorithm (see Section 1.3). Thus, we can compute x0 with k + 1 calls to GGT on a network with n nodes and m arcs. This time is dominated by the time required to compute W and xτ for all τ ∈ W. This part of the algorithm is summarized in Step 3 of Figure 5. The following theorem summarizes the above discussions of the steps in computing x∗ : Theorem 15. A universally maximum flow with piecewise-constant capacities can be constructed in O(k2 mn log(kn2 /m)) time. Acknowledgments I am grateful to Éva Tardos and Kevin Wayne for providing helpful comments on early drafts of this paper and to an anonymous referee for useful comments. REFERENCES [1] E.J. Anderson and P. Nash, Linear programming in infinitedimensional spaces, Wiley, New York, 1987. [2] E.J. Anderson, P. Nash, and A.B. Philpott, A class of continuous network flow problems. Math Oper Res 7 (1982), 501–514. [3] E.J. Anderson and A.B. Philpott, A continuous-time network simplex algorithm, Networks 19 (1989), 395–425. [4] J.E. Aronson, A survey of dynamic network flows, Ann Oper Res 20 (1989), 1–66. [5] R.E. Burkard, K. Dlaska, and B. Klinz, The quickest flow problem, ZOR Methods Models Oper Res 37 (1993), 31– 58. [6] L. Fleischer, Faster algorithms for the quickest transshipment problem, SIAM J Option 12 (2001), 18–35 (extended abstract appeared in Proc SODA 1998). [7] L. Fleischer and É. Tardos, Efficient continuous-time dynamic network flow algorithms, Oper Res Lett 23 (1998), 71–80. [8] L.R. Ford and D.R. Fulkerson. Flows in networks, Princeton University Press, Princeton, NJ, 1962. [9] G. Gallo, M.D. Grigoriadis, and R.E. Tarjan, A fast parametric maximum flow algorithm and applications, SIAM J Comput 18 (1989), 30–55.
[10] A.V. Goldberg and R.E. Tarjan, A new approach to the maximum flow problem, J ACM 35 (1988), 921–940. [11] D. Gusfield and C. Martel, A fast algorithm for the generalized parametric minimum cut problem and applications, Algorithmica 7 (1992), 499-519. [12] B. Hajek and R.G. Ogier, Optimal dynamic routing in communication networks with continuous traffic, Networks 14 (1984), 457–487. [13] D.S. Hochbaum and S.P. Hong, About strongly polynomial time algorithms for quadratic optimization over submodular constraints, Math Program 69 (1995), 269–309. [14] B. Hoppe, Efficient dynamic network flow algorithms, PhD Thesis, Cornell University, June 1995, Department of Computer Science Technical Report TR95-1524. [15] B. Hoppe and É. Tardos, Polynomial time algorithms for some evacuation problems, Proc 5th Annual ACM-SIAM Symp on Discrete Algorithms, 1994, pp. 433–441. [16] B. Hoppe and É. Tardos, The quickest transshipment problem, Math Oper Res 25 (2000), 36–62 (extended abstract appeared in Proc SODA 1995). [17] B. Klinz and G.J. Woeginger, “Minimum cost dynamic flows: The series-parallel case,” Integer programming and combinatorial optimization, E. Balas and J. Clausen (Editors), Proc 4th Int IPCO Conf, Lecture Notes in Computer Science 920, Springer-Verlag, 1995, pp. 329–343. [18] S.T. McCormick, Fast algorithms for parametric scheduling come from extensions to parametric maximum flow, Oper Res 47 (1999), 744–756 (extended abstract appeared in Proc STOC 96). [19] E. Minieka, Maximal, lexicographic, and dynamic network flows, Oper Res 21 (1973), 517–527. [20] F.H. Moss and A. Segall, An optimal control approach to dynamic routing in networks, IEEE Trans Autom Control 27 (1982), 329–339. [21] R.G. Ogier, Minimum-delay routing in continuous-time dynamic networks with piecewise-constant capacities, Networks 18 (1988), 303–318. [22] J.B. Orlin, Maximum-throughput dynamic network flows, Math Program 27 (1983), 214–231. [23] J.B. Orlin, Minimum convex cost dynamic network flows, Math Oper Res 9 (1984), 190–207. [24] A.B. Philpott, Continuous-time flows in networks, MOR 15 (1990), 640–661. [25] A.B. Philpott and M. Craddock, An adaptive discretization method for a class of continuous network programs, Networks 26 (1995), 1–11. [26] W.B. Powell, P. Jaillet, and A. Odoni, “Stochastic and dynamic networks and routing,” Handbooks in operations research and management science: Networks, M.O. Ball, T.L. Magnanti, C.L. Monma, and G.L. Nemhauser (editors), Elsevier, 1995. [27] M.C. Pullan, An algorithm for a class of continuous linear programs, SIAM J Control Optim 31 (1993), 1558–1577. [28] M.C. Pullan, A study of general dynamic network programs with arc time-delays, SIAM J Optim 7 (1997), 889– 912. [29] G.I. Stassinopoulos and P. Konstantopoulos, Optimal congestion control in single destination networks, IEEE Trans Commun 33 (1985), 792–800. [30] W.L. Wilkinson, An algorithm for universal maximal dynamic flows in a network, Oper Res 19 (1971), 1602–1612.
NETWORKS–2001
125
