Discrete Newton methods for the evacuation problem

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.1 (1-10)

Theoretical Computer Science ••• (••••) •••–•••

Contents lists available at ScienceDirect

Theoretical Computer Science www.elsevier.com/locate/tcs

Discrete Newton methods for the evacuation problem Naoyuki Kamiyama Institute of Mathematics for Industry, Kyushu University, Fukuoka, Japan

a r t i c l e

i n f o

Article history: Received 4 September 2017 Received in revised form 29 April 2018 Accepted 2 August 2019 Available online xxxx Communicated by L.A. Gasieniec Keywords: Dynamic network flow Evacuation problem Submodular function minimization

a b s t r a c t A dynamic network is a directed graph in which arcs have capacities and transit times. In this paper, we consider the evacuation problem in dynamic networks. In this problem, we are given a dynamic network with a single sink vertex in which each vertex except the sink vertex has a supply. Then the goal of this problem is to find the minimum time limit T such that we can send all the supplies to the sink vertex by T . In this paper, we propose a discrete Newton method for the evacuation problem. First, we prove that the number of iterations of this method is at most the number of vertices of the input dynamic network. Then we propose theoretical and practical implementation of this method. The theoretical implementation is based on submodular function minimization, and the practical implementation is based on maximum flow computations in timeexpanded networks. Finally, we compare the proposed practical implementation with an algorithm using a binary search in time-expanded networks in computational experiments. © 2019 Elsevier B.V. All rights reserved.

1. Introduction A dynamic network, which was introduced by Ford and Fulkerson [1,2], is a directed graph in which arcs have capacities and transit times (see, e.g., [3] for a detailed introduction to dynamic networks). A dynamic network is frequently used for modeling evacuation situations [4–7] (see, e.g., [8] for a survey of models of evacuation situations based on dynamic networks). One of the most fundamental problems in dynamic networks is the evacuation problem. In this problem, we are given a dynamic network with a single sink vertex and a supply for each vertex except the sink vertex. Then the goal of the evacuation problem is to find the minimum time limit T such that we can send all the supplies to the sink vertex by T . It is known [9] that we can solve the evacuation problem in polynomial time by using a submodular function minimization algorithm (e.g., [10,11]) as a subroutine in a binary search or a parametric search [12] (see also [13]). The aim of this paper is to propose a discrete Newton method [14,15] for the evacuation problem. In a binary search, we first estimate an upper bound of the optimal objective value, and then we iteratively reduce the interval from some lower bound to this upper bound. On the other hand, in a discrete Newton method, we start with some lower bound of the optimal objective value, and then we check whether the current lower bound is also an upper bound. If the answer is yes, then this lower bound is the optimal objective value. Otherwise, we update the current lower bound. In a discrete Newton method, we do not need to solve a problem with respect to a parameter that is larger than the optimal objective value. This is an important point in the evacuation problem because its optimal objective value may be very large. First, we propose a discrete Newton method for the evacuation problem (Section 3), and then we prove that the number of iterations of this method is at most the number of vertices of the input dynamic network (Section 4). Then we propose theoretical and practical implementation of this method (Sections 5 and 6). The theoretical implementation is based on

E-mail address: [email protected]. https://doi.org/10.1016/j.tcs.2019.08.004 0304-3975/© 2019 Elsevier B.V. All rights reserved.

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.2 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

2

submodular function minimization, and the practical implementation is based on maximum flow computations in timeexpanded networks. Lastly, we compare the proposed practical implementation with an algorithm using a binary search in time-expanded networks in computational experiments (Section 7). It should be noted that in [16, Chapter 4], Schlöter independently proposed an algorithm for the evacuation problem which is almost the same as the algorithm proposed in Section 5. 2. Preliminaries We denote by R, R+ , Z+ , and Z++ the sets of real numbers, non-negative real numbers, non-negative integers, and positive integers, respectively. For each finite directed graph H = ( N , L ), where N is the set of vertices and L is the set of arcs, and each subset X of N, we denote by δ( X | H ) (resp., ( X | H )) the set of arcs a = (u , v ) in L such that u ∈ X and v∈ / X (resp., u ∈ / X and v ∈ X ). For each vertex v in N, we write δ( v | H ) and ( v | H ) instead of δ({ v }| H ) and ({ v }| H ), respectively. 2.1. Dynamic networks and the evacuation problem The evacuation problem is defined as follows. In this problem, we are given a finite simple directed graph D = ( V , A ), a capacity function c : A → Z++ , a transit time function τ : A → Z++ , and a sink vertex r in V . Define n := | V |, m  := | A |, and U := V \ {r }. Furthermore, we are given a supply function b : U → Z+ . For each subset X of U , we define b( X ) := v ∈ X b( v ). In this paper, we assume that there exists at most one of arcs a = (u , v ) and a = ( v , u ) for each pair of distinct vertices u , v in V . Furthermore, we assume that the sink vertex is reachable from every vertex in V , and δ(r | D ) = ∅. A function f : A × Z+ → R+ is called a dynamic flow in D. For each arc a = (u , v ) in A and each integer θ in Z+ , f (a, θ) represents a flow rate entering u at the time θ which reaches v at the time θ + τ (a). For each dynamic flow f in D, each vertex v in V , and each integer θ in Z+ , we define

ex f ( v , θ) :=



θ− τ (a) 

a∈( v | D )

t =0

f (a, t ) −

θ  

f (a, t ).

a∈δ( v | D ) t =0

Intuitively speaking, ex f ( v , θ) represents the excess (or deficiency) of flow at the vertex v at the time θ . For each integer θ in Z+ , a dynamic flow f in D is said to be feasible with respect to θ if the following conditions are satisfied. (D1) For every arc a in A, the following conditions are satisfied. a) For every integer t in Z+ such that t ≤ θ − τ (a), we have f (a, t ) ≤ c (a). b) For every integer t in Z+ such that t > θ − τ (a), we have f (a, t ) = 0. (D2) For every vertex v in U , the following conditions are satisfied. a) For every integer t in Z+ such that t < θ , we have ex f ( v , t ) ≥ −b( v ). b) ex f ( v , θ) = −b( v ). The goal of the evacuation problem is to compute the minimum integer θ in Z+ such that there exists a feasible dynamic flow with respect to θ in D. Define UB := n · maxa∈ A τ (a) + b(U ). Then it is known [9] that the optimal solution is at most UB. Although in this paper we do not consider the problem of finding a feasible dynamic flow with respect to the optimal solution of the evacuation problem, it should be noted that we can find such a dynamic flow in polynomial time by using the algorithm in [9, Section 6]. 2.2. Time-expanded networks In this subsection, we define a time-expanded network, which was introduced by Ford and Fulkerson [1,2]. For each integer θ in Z+ , we define the time-expanded network D θ = ( V θ , A θ ) with respect to θ as follows. The vertex set V θ contains a super-source vertex s∗ and a vertex v t for each vertex v in V and each integer t in {0, 1, . . . , θ}. The arc set A θ is the union of A θ,1 , A θ,2 , and A θ,3 defined as follows. The arc set A θ,1 contains an arc at from ut to v t +τ (a) for each arc a = (u , v ) in A and each integer t in {0, 1, . . . , θ − τ (a)}. The arc set A θ,2 contains an arc htv from v t to v t +1 for each vertex v in V and each integer t in {0, 1, . . . , θ − 1}. The arc set A θ,3 contains an arc a v from s∗ to v 0 for each vertex v in U . Furthermore, we define the capacity function c θ : A θ → R+ ∪ {∞} as follows. For each arc at in A θ,1 , we define c θ (at ) := c (a). For each arc htv in A θ,2 , we define c θ (htv ) := ∞. For each arc a v in A θ,3 , we define c θ (a v ) := b( v ). Assume that we are given an integer θ in Z+ . A function ξ : A θ → R+ is called a feasible static flow in D θ if the following conditions are satisfied. (S1) For every arc  in A θ , we have ξ() ≤ c θ ().

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.3 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

3

(S2) For every vertex v in V θ \ {s∗ , rθ }, we have



ξ() =

∈δ( v | D θ )



ξ().

∈( v | D θ )

The value val(ξ ) of a feasible static flow ξ in D θ is defined by



val(ξ ) :=

ξ().

∈δ(s∗ | D θ )

Notice that val(ξ ) ≤ b(U ) for every feasible static flow ξ in D θ . Then a feasible static flow ξ in D θ such that val(ξ ) is maximum among all the feasible static flows in D θ is called a maximum static flow in D θ . If we view the time-expanded network as a flow network with the source s∗ and the sink rθ , then the notions of feasible and maximum static flows and the value of a feasible static flow are the standard notions from the context of the (static) maximum flow problem, where (S1) are capacity constraints and (S2) are flow conservation constraints. See, e.g., [17]. Lemma 1 (Ford and Fulkerson [1,2]). Assume that we are given an integer θ in Z+ . Then there exists a feasible dynamic flow with respect to θ in D if and only if there exists a feasible static flow ξ in D θ such that val(ξ ) = b(U ). Assume that we are given an integer θ in Z+ . Then Lemma 1 implies that by finding a maximum static flow in D θ , we can check whether there exists a feasible dynamic flow with respect to θ in D. This implies that by computing a maximum static flow in a time-expended network O (log UB) times, we can solve the evacuation problem. However, the size of D θ is not bounded by a polynomial in the input size of the evacuation problem. Thus, this algorithm is not a polynomial-time algorithm. 3. Discrete Newton methods for the evacuation problem In this section, we propose a discrete Newton method for the evacuation problem. First, we give an intuitive explanation about it. To propose a discrete Newton method, we introduce a function dθ : 2U → R for each integer θ in Z+ . The functions dθ have the following properties. (i) For every subset X of U , the value dθ ( X ) is a non-decreasing function with respect to θ . (ii) For every integer θ in Z+ , the optimal solution of the evacuation problem is at most θ if and only if min X ⊆U dθ ( X ) ≥ 0. Then our method works as follows. First, we compute the minimum integer θ in Z+ such that dθ (U ) ≥ 0. The property (ii) implies that θ is a lower bound of the optimal solution of the evacuation problem. If min X ⊆U dθ ( X ) ≥ 0, then the optimal objective value of the evacuation problem is at most θ , i.e., θ is the optimal solution. Otherwise, there exists a subset X of U such that dθ ( X ) < 0. In this case, we compute the minimum integer θ  in Z+ such that dθ  ( X ) ≥ 0. Then θ  is a lower bound of the optimal solution of the evacuation problem. If min X ⊆U dθ  ( X ) ≥ 0, then θ  is the optimal solution of the evacuation problem. We repeat this procedure. In Section 4, we prove that the number of iterations of this method can be bounded by a polynomial in the input size. In order to define the functions dθ , we first introduce a new directed graph E X for each subset X of U . This directed graph was introduced in [1,2] to reduce the problem of finding a maximum dynamic flow in a dynamic network to the minimum-cost circulation problem in a static network. Assume that we are given a subset X of U . Then the directed graph E X := ( V X , A X ) is defined as follows. Define V X := V ∪ { s}, where  s is a new vertex. The arc set A X is defined by

A X := A ∪ {e v = ( s, v ) | v ∈ X } ∪ {er = (r , s)}. The capacity function c X : A X → R+ ∪ {∞} is defined by



X

c (a) :=

c (a)

if a ∈ A

∞

if a ∈ A X \ A .

For each integer θ in Z+ , we define the cost function

πθX : A X → R by

⎧ ⎪ if a ∈ A ⎨τ (a) πθX (a) := −θ − 1 if a = er ⎪ ⎩ 0 if a = e v for some vertex v in X .

A function

ϕ : A X → R+ is called a circulation in E X if the following conditions are satisfied.

(C1) For every arc a in A X , we have

ϕ (a) ≤ c X (a).

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.4 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

4

(C2) For every vertex v in V X , we have





ϕ (a) =

a∈δ( v | E X )

The cost costθX (ϕ ) of a circulation

costθX (ϕ ) :=



ϕ (a).

a∈( v | E X )

ϕ in E X with respect to πθX is defined by

πθX (a) · ϕ (a).

a∈ A X

For each integer θ in Z+ , we define the function oθ : 2U → R+ by

oθ ( X ) := − min{costθX (ϕ ) | ϕ is a circulation in E X }. Notice that the notions here are the standard notions from the context of the (static) minimum-cost circulation problem. See, e.g., [17]. It is known [18] that for every integer θ in Z+ and every subset X of U , we can compute o θ ( X ) in O (m2 log2 n) time.1 For each integer θ in Z+ , we define the function dθ : 2U → R by

dθ ( X ) := oθ ( X ) − b( X ). We are now ready to propose a discrete Newton method for the evacuation problem. See Algorithm 1. Algorithm 1: A discrete Newton method for the evacuation problem. 1 2 3 4 5 6 7

Compute the minimum integer θ1 in Z+ such that dθ1 (U ) ≥ 0. Find a minimal minimizer X 1 of dθ1 . Set i := 1. while dθi ( X i ) < 0 do Compute the minimum integer θi +1 in Z+ such that dθi+1 ( X i ) ≥ 0. Find a minimal minimizer X i +1 of dθi+1 . Set i := i + 1.

8 end 9 Output θi (i.e., θi is the optimal solution), and halt.

The correctness of Algorithm 1 follows from the following lemmas. Lemma 2 (Klinz [9, Theorem 5.1]). Assume that we are given an integer θ in Z+ . Then there exists a feasible dynamic flow with respect to θ in D if and only if oθ ( X ) ≥ b( X ) for every subset X of U . Lemma 3. Assume that Algorithm 1 halts when i = k. Then θk is the optimal solution of the evacuation problem. Proof. Lemma 2 and the definition of Algorithm 1 imply that θi is a lower bound of the optimal solution of the evacuation problem for every integer i in {1, 2, . . . , k}. Thus, since dθk ( X ) ≥ 0 for every subset X of U , Lemma 2 implies that θk is the optimal solution of the evacuation problem. This completes the proof. 2 4. The number of iterations In this section, we evaluate the number of iterations of Algorithm 1. Assume that we are given a finite set S and a function g : 2 S → R. Then g is said to be submodular if

g ( X ) + g (Y ) ≥ g ( X ∪ Y ) + g ( X ∩ Y ) for every pair of subsets X , Y of S. It is not difficult to see that if g is submodular, then the family of minimizers of g is closed under union and intersection, and thus there exists the unique minimal minimizer of g. Assume that we are given a finite set S and submodular functions g 1 , g 2 : 2 S → R. Then we write g 1 ← g 2 if for every pair of subsets X , Y of S such that X ⊆ Y ,

g 1 (Y ) − g 1 ( X ) ≤ g 2 (Y ) − g 2 ( X ). 1 Precisely speaking, the original evaluation of the time complexity of the algorithm of [18] is O ((m log n)(m + n log n)). In this paper, for ease of presentation, we use this weaker bound.

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.5 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

5

Lemma 4 (Topkis [19]). For every finite set S and every pair of submodular functions g 1 , g 2 : 2 S → R such that g 1 ← g 2 , the minimal minimizer of g 2 is a subset of the minimal minimizer of g 1 . Lemma 5 (Hoppe and Tardos [9, p. 51]). For every integer θ in Z+ , oθ is a submodular function. Lemma 6 (Baumann and Skutella [20, Lemma 3.1]). For every pair of integers θ, θ  in Z+ such that θ ≤ θ  , we have oθ ← oθ  .2 Lemma 5 implies that dθ is a submodular function. Furthermore, Lemma 6 implies that for every pair of integers θ, θ  in Z+ such that θ ≤ θ  , we have dθ ← dθ  . We are now ready to prove an upper bound of the number of iterations of Algorithm 1. Lemma 7. Algorithm 1 halts when i = k for some positive integer k such that k ≤ n. Proof. It is not difficult to see that for every subset X of U and every pair of integers θ, θ  in Z+ , if oθ ( X ) < oθ  ( X ), then θ < θ  . This implies that during the course of Algorithm 1, since oθi ( X i ) < b( X i ) ≤ oθi+1 ( X i ), θi < θi +1 . Thus, Lemmas 4, 5, and 6 imply that during the course of Algorithm 1, X i +1 ⊆ X i . Next, we prove that if Algorithm 1 does not halt when i = t for some positive integer t such that t ≥ 2, then X t −1 = X t . The definition of Algorithm 1 implies that dθt ( X t ) < 0. Furthermore, the definition of θt implies that dθt ( X t −1 ) ≥ 0. These imply that X t −1 = X t . Thus, the number of iterations of Algorithm 1 is finite. If Algorithm 1 halts when i = k for some positive integer k such that k ≥ 2, then X k = ∅ and X k−1 = ∅. This implies that X i = X i +1 during the course of Algorithm 1. Thus, since | X 1 | ≤ n − 1, this completes the proof. 2 5. Theoretical implementation In this section, we propose theoretical implementation of Algorithm 1, which has a similar flavor as the algorithm of [20] for the earliest arrival flow problem. In fact, Baumann and Skutella [20] suggested using parametric submodular function minimization in their algorithm. In the rest of this section, we assume that Algorithm 1 halts when i = k. Before considering the whole part of Algorithm 1, we first consider the time required to compute θi for each integer i in {1, 2, . . . , k}. Based on the algorithm of Lin and Jaillet [21], Saho and Shigeno [22] proposed a strongly polynomial-time algorithm for computing θi . Lemma 8 (Saho and Shigeno [22, Theorem 11]). For every subset X of U , we can compute the minimum integer θ in Z+ such that dθ ( X ) ≥ 0 in O (nm2 log2 n) time.3 First, we propose naive implementation of Algorithm 1 based on the following lemma. Lemma 9 (Lee, Sidford, and Wong [23]). For every finite set S and every submodular function g : 2 S → R such that g (∅) = 0, we can find a minimizer of g in O (| S |3 (log2 | S |)EO + | S |4 log O (1) | S |) time, where EO is the time required to evaluate g ( X ) for each subset X of S. Theorem 1. Algorithm 1 can be implemented in O (n5m2 log4 n + n6 log O (1) n) time. Proof. First, Lemma 9 implies that we can find a minimizer of dθi in O (n3 m2 log4 n + n4 log O (1) n) time. Furthermore, it is known (see, e.g., [24, Note 10.12]) that the minimal minimizer of a submodular function g : 2 S → R such that g (∅) = 0 can be computed by using any submodular function minimization algorithm O (| S |) times. Thus, Lemmas 7 and 8 imply this theorem. 2 Here we give a remark about the time complexity in Theorem 1. If we use the algorithm of Orlin [25] for the submodular function minimization problem, we can find a minimizer X ∗ of dθi in O (n5 m2 log2 n) time. Although this algorithm is not the current fastest one, it is known (see, e.g., [24, Note 10.11] and [26, Lemma 1]) that by using the information about X ∗ in the algorithm of Orlin [25], we can find the minimal minimizer X i of dθi in O (n3 m2 log2 n) time. Thus, we do not need to use a submodular function minimization algorithm O (n) times. This implies that if we can do the same thing by using the algorithm of [23], then the time complexity in Theorem 1 can be improved. Since this topic is beyond the scope of this paper, we do not go into the details.

2 Precisely speaking, Baumann and Skutella [20] considered the continuous version of dynamic flows. However, it is not difficult to see that this lemma is also valid for the discrete version, i.e., the setting that is considered in this paper. 3 As in the paper [20], Saho and Shigeno [22] considered the continuous version of dynamic flows. However, it is not difficult to see that this lemma is also valid for the discrete version.

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.6 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

6

Next, we consider implementation based on a parametric submodular function minimization algorithm. Nagano [26] proposed a parametric submodular function minimization algorithm based on the algorithm of Orlin [25]. (See also [27] for another parametric submodular function minimization algorithm, which is not the current fastest one.) Lemma 10 (Nagano [26, Corollary 14]). Assume that we are given a finite set S and submodular functions g 1 , g 2 , . . . , g z : 2 S → R such that (i) z = O (| S |2 ), (ii) gt (∅) = 0 for every integer t in {1, 2, . . . , z}, and (iii) g 1 ← g 2 ← · · · ← g z . Furthermore, we assume that the minimal minimizers of the submodular functions g 1 , g 2 , . . . , g z are Y 1 , Y 2 , . . . , Y z , respectively. Then we can find the minimal minimizers Y 1 , Y 2 , . . . , Y z in O (| S |5 EO + | S |6 ) time, where EO is the time required to evaluate gt ( X ) for each integer t in {1, 2, . . . , z} and each subset X of S. This algorithm can be run even if g 1 , g 2 , . . . , g z are obtained in this order in an on-line manner during the course of the algorithm. We are now ready to give the main result of this section. Theorem 2. Algorithm 1 can be implemented in O (n5m2 log2 n) time. Proof. Lemmas 6, 7, and 10 imply that the total time complexity of Steps 2 and 6 of Algorithm 1 is O (n5 m2 log2 n). Lemmas 7 and 8 imply that the total time complexity of Steps 1 and 5 is O (n2 m2 log2 n). Thus, the time complexity of Algorithm 1 is O (n5 m2 log2 n). 2 Notice that the algorithm of [25] is not the current fastest submodular function minimization algorithm. The current fastest algorithm is due to [23]. If we can extend Lemma 10 by using the algorithm of [23], then Algorithm 1 will become much faster. Since this topic is beyond the scope of this paper, we do not go into the details. The time complexity of an algorithm using a binary search and a submodular function minimization algorithm depends on maxa∈ A τ (a) and b(U ). Thus, since the time complexities of the proposed algorithms depend on only n and m, we cannot compare the proposed algorithms with this algorithm. Furthermore, in the papers [9,13], algorithms based on a parametric search were proposed. However, since the precise time complexities of these algorithms were not given, we cannot compare the proposed algorithms with these algorithms. 6. Practical implementation In this section, we propose practical implementation of Algorithm 1. In practical applications, algorithms based on timeexpanded networks are frequently used instead of algorithms based on submodular function minimization algorithms. This is because algorithms based on time-expanded networks can be easily implemented, and there exists a good package for computing a maximum static flow, e.g., NetworkX package [28] for Python. Furthermore, algorithms based on time-expanded networks have the merit that we can easily find not only the minimum evacuation time but also an optimal flow itself. In this section, we propose practical implementation of Algorithm 1 based on time-expanded networks and maximum static flow computations. In our practical implementation, for each subset X of U , we use a binary search instead of the algorithm of [22] to compute the minimum integer θ in Z+ such that dθ ( X ) ≥ 0 (i.e., Steps 1 and 5 of Algorithm 1). Notice that each iteration of this binary search computes a minimum-cost circulation in E X , which is of a similar size as the input dynamic network. Since the input dynamic network is much smaller than the time-expanded network, this algorithm is not time-consuming. Furthermore, since the naive upper bound UB can be much larger than the optimal solution of the evacuation problem, we use the following well-known Algorithm 2 to compute the minimum integer θ in Z+ such that dθ ( X ) ≥ 0. Notice that the number of iterations of Algorithm 2 is O (log T ), where T is the optimal solution of the evacuation problem. Algorithm 2: An algorithm for computing the minimum integer θ in Z+ such that dθ ( X ) ≥ 0. 1 2 3 4 5

If b( X ) = 0, then output 0 and halt. Otherwise, set t := 1. while dt ( X ) < 0 do Set t := 2t. end Compute the minimum integer θ in Z+ such that dθ ( X ) ≥ 0 by a binary search using t and t /2 as initial upper and lower bounds, respectively.

Notice that Lemmas 1 and 2 imply that for every integer θ in Z+ and the minimal minimizer X of dθ , dθ ( X ) ≥ 0 holds if and only if val(ξ ) = b(U ) for a maximum static flow ξ in D θ . Thus, Step 4 of Algorithm 1 can be implemented by finding a maximum static flow in a time-expanded network. In the rest of this section, we are given an integer θ in Z+ . Then we consider how to compute the minimal minimizer of dθ by using a maximum static flow algorithm in a time-expanded network. Assume that we are given a feasible static flow ξ in D θ . Define the directed graph D θ (ξ ) = ( V θ , A θ (ξ )) as follows. For each arc  = (u , v ) in A θ , we define the arc −1 by −1 := ( v , u ). The arc set A θ (ξ ) is the union of A ◦θ (ξ ) and A •θ (ξ ) defined

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.7 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

7

as follows. Define A ◦θ (ξ ) as the set of arcs  in A θ such that ξ() < c θ (). Define A •θ (ξ ) as the set of arcs −1 such that  ∈ A θ and ξ() > 0. The proposed algorithm for computing the minimal minimizer of dθ by using a maximum static flow algorithm in D θ is described in Algorithm 3. Algorithm 3: An algorithm for finding the minimal minimizer of dθ in D θ . 1 2 3 4

Find a maximum static flow ξθ in D θ . Define Z θ as the set of vertices in V θ reachable from s∗ in D θ (ξθ ). Define X θ := { v ∈ U | v 0 ∈ Z θ }. Output X θ , and halt.

From here, we prove the correctness of Algorithm 3. A subset Z of V θ is called a cut in D θ if s∗ ∈ Z and rθ ∈ / Z . For each cut Z in D θ , we define the capacity cap( Z ) of Z by



cap( Z ) :=

c θ ().

∈δ( Z | D θ )

Then a cut Z in D θ is called a minimum cut if cap( Z ) is minimum among all the cuts in D θ . For each subset X of U , we define D X as the family of cuts Z in D θ such that v 0 ∈ Z for every vertex v in X . For each subset X of U , we define C X as the family of members Z of D X such that v 0 ∈ / Z holds for every vertex v in U \ X . For each subset X of U , we define mθ ( X ) := min{cap( Z ) | Z ∈ C X }. Lemma 11. For every subset X of U ,

oθ ( X ) = mθ ( X ) − b(U \ X ). Proof. Assume that we are given a subset X of U . Then we define the function c θ : A θ → R+ ∪ {∞} by



c θ () :=

0

if  = a v for some vertex v in U \ X

c θ ()

otherwise.

Furthermore, for each cut Z in D θ , we define

cap( Z ) :=



c θ ().

∈δ( Z | D θ )

It is well known [1,2] that

oθ ( X ) = min{cap( Z ) | Z ∈ D X }. Precisely speaking, this result follows from the so-called max-flow min-cut theorem and the result on the maximum dynamic flow problem in [1,2]. If there exists a member Z ∗ of C X such that

cap( Z ∗ ) = min{cap( Z ) | Z ∈ D X },

(1)

then we have

oθ ( X ) = min{cap( Z ) | Z ∈ D X } = min{cap( Z ) | Z ∈ C X } = min{cap( Z ) − b(U \ X ) | Z ∈ C X }

= min{cap( Z ) | Z ∈ C X } − b(U \ X ) = mθ ( X ) − b(U \ X ). This completes the proof. Here we prove that there exists a member Z ∗ of C X satisfying (1). Define M := min{cap( Z ) | Z ∈ D X }. Let Z ∗ be a member of D X such that cap( Z ∗ ) = M. Assume that it minimizes the number of vertices v in U \ X such that v 0 ∈ Z ∗ among all the members Z of D X satisfying the condition that cap( Z ) = M. If v 0 ∈ / Z ∗ holds for every vertex v in U \ X , ∗ ∗ then since Z ∈ C X , the proof is done. Assume that v 0 ∈ Z for some vertex v in U \ X . Since τ (a) > 0 for every arc a in A, ( v 0 | D θ ) = {a v }. Thus, since c θ (a v ) = 0,

cap( Z ∗ ) ≥ cap( Z ∗ \ { v 0 }). Since Z ∗ \ { v 0 } ∈ D X , this implies that cap( Z ∗ \ { v 0 }) = M. This contradicts the definition of Z ∗ . This completes the proof.

2

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.8 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

8

Lemma 12. For every subset X of U and every cut Z in D θ such that Z ∈ C X and cap( Z ) = mθ ( X ), we have

dθ ( X ) = cap( Z ) − b(U ). Proof. Lemma 11 implies that

dθ ( X ) = oθ ( X ) − b( X ) = mθ ( X ) − b(U ) = cap( Z ) − b(U ). This completes the proof.

2

Lemma 13. For every subset X of U , there exists a minimum cut Z in D θ such that Z ∈ C X if and only if X is a minimizer of dθ . Proof. Assume that we are given a subset X of U . First, we prove the “if” part. Assume that X is a minimizer of dθ . Let Z be a cut in D θ such that Z ∈ C X and cap( Z ) = mθ ( X ). We prove that Z is a minimum cut in D θ . Let K be a cut in D θ . Define Y as the set of vertices v in U such that v 0 ∈ K . Let L be a cut in D θ such that L ∈ CY and cap( L ) = mθ (Y ). The definition of L implies that cap( L ) ≤ cap( K ). Furthermore, since X is a minimizer of dθ , Lemma 12 implies that cap( Z ) ≤ cap( L ). These imply that cap( Z ) ≤ cap( K ). This completes the proof. Next, we prove the “only if” part. Assume that there exists a minimum cut Z in D θ such that Z ∈ C X . Notice that cap( Z ) = mθ ( X ). Let Y be a subset of U . We prove that dθ ( X ) ≤ dθ (Y ). Let K be a cut in D θ such that K ∈ CY and cap( K ) = mθ (Y ). Then since Z is a minimum cut in D θ , we have cap( Z ) ≤ cap( K ). Thus, Lemma 12 implies that dθ ( X ) ≤ dθ (Y ). This completes the proof. 2 Lemma 14 (Picard and Queyranne [29, p. 11]). Z θ is a minimal minimum cut in D θ .4 Since it is not difficult to see that the family of minimum cuts in D θ is closed under union and intersection, Lemma 14 implies that Z θ ⊆ K for every minimum cut K in D θ . Lemma 15. X θ is the minimal minimizer of dθ . Proof. Since Lemmas 13 and 14 imply that X θ is a minimizer of dθ . Assume that there exists a minimizer Y of dθ such that Y  X θ . Lemma 13 implies that there exists a minimum cut K in D θ such that K ∈ CY . Thus, since Y  X θ , there exists a vertex v in X θ \ Y such that v 0 ∈ Z θ \ K . This implies that Z θ  K , which contradicts Lemma 14. This completes the proof. 2 The proposed practical implementation of Algorithm 1 is described in Algorithm 4. Algorithm 4: Practical implementation of Algorithm 1. 1 2 3 4 5 6 7 8 9

Compute the minimum integer θ1 in Z+ such that dθ1 (U ) ≥ 0 by Algorithm 2. Find a maximum static flow ξ1 in D θ1 . Define Z 1 as the set of vertices in V θ1 reachable from s∗ in D θ1 (ξ1 ). Define X 1 := { v ∈ U | v 0 ∈ Z 1 }, and set i := 1. while val(ξi ) < b(U ) do Compute the minimum integer θi +1 in Z+ such that dθi+1 ( X i ) ≥ 0 by Algorithm 2. Find a maximum static flow ξi +1 in D θi+1 . Define Z i +1 as the set of vertices in V θi+1 reachable from s∗ in D θi+1 (ξi +1 ). Define X i +1 := { v ∈ U | v 0 ∈ Z i +1 }, and set i := i + 1.

10 end 11 Output θi (i.e., θi is the optimal solution), and halt.

7. Computational experiments In this section, we compare Algorithm 4 with an algorithm based on a binary search in computational experiments. First, we give an algorithm for the evacuation problem based on a binary search as the baseline of our experiments. Since UB can be very large, we first estimate a smaller upper bound of the optimal solution. Then we do a binary search. This algorithm is described in Algorithm 5. In this algorithm, the number of computations of a maximum static flow in a time-expanded network is (log T ), where T is the optimal solution of the evacuation problem. Notice that in each of

4

This statement is rephrased to fit this paper. Picard and Queyranne [29] proved this result for a general network.

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

[m3G; v1.261; Prn:16/09/2019; 10:58] P.9 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

9

Table 1 The results of the computational experiments.



0.25

US

25

50

100

25

0.5 50

100

25

1 50

100

N = 10

Average Min Max

2.41 1.32 5.66

3.63 1.24 8.59

4.64 1.98 11.41

3.31 1.28 11.69

4.60 1.73 11.20

5.11 2.28 11.89

4.39 1.78 10.96

5.75 2.37 12.02

6.27 2.67 13.23

N = 15

Average Min Max

3.13 1.39 11.19

4.42 1.75 10.44

5.64 2.70 13.34

4.59 2.32 8.27

5.62 2.61 11.77

5.92 2.78 12.02

5.69 2.78 11.15

6.16 2.49 12.45

6.93 2.55 13.41

the (log T ) iterations of the binary search in Step 5, we have to find a maximum static flow in a time-expanded network for θ = ( T ). Algorithm 5: An algorithm using a binary search. 1 2 3 4 5

If b(U ) = 0, then output 0 and halt. Otherwise, set t := 1. while the maximum value of a static flow in D t is less than b(U ) do Set t := 2t. end Compute the minimum integer θ in Z+ such that the maximum value of a static flow in D θ is b(U ) by a binary search using t and t /2 as initial upper and lower bounds, respectively.

The algorithms were implemented by Python 3.5.1 with NetworkX 1.11 [28]. The experimental results were obtained on a 64-bit machine with Windows 7, Intel Core i7 (3.6 GHz), and 32 GB of RAM, and the running time was measured in seconds. In the implementation of Algorithms 4 and 5, we used preflow_push command in NetworkX package to compute a maximum static flow. Furthermore, in the implementation of Algorithm 4, we used min_cost_flow_cost command to compute the minimum cost of a circulation. To the best of our knowledge, there exists no benchmark instance of the evacuation problem. Thus, we randomly made instances in the following setting.

• The underlying graph is an N × N grid. Thus, the number of vertices is N 2 , and the number of arcs is 2N ( N − 1). Network topologies appearing in real applications are frequently close to grids (see, e.g., [7]). In the following experiments, we set N := 10, 15. For example, Zheng, Chiu, and Mirchandani [30] used a network with 180 vertices. Thus, the sizes of our instances seem to be moderate for practical applications. • We randomly choose some vertex as a sink vertex. Then we direct all the arcs toward this sink vertex. • For each arc, its capacity and transit time are randomly chosen from {1, 2, . . . , 10} by using random.uniform command in Python. • The supply of each vertex except the sink vertex is determined as follows. First, we choose a parameter , and then we give positive supplies to  · N 2  vertices. In the following experiments, we set  := 0.25, 0.5, 1. Furthermore, we choose an upper bound US, and randomly choose each supply from {1, 2, . . . , US} by using random.uniform command in Python. In the following experiments, we set US := 25, 50, 100. In our experiments, we made 100 instances for each combination of N, , and US. Then we solved them by using Algorithms 4 and 5. We consider the ratio

the seconds required to solve the instance by Algorithm 5 the seconds required to solve the instance byAlgorithm 4

(2)

for each instance. Table 1 shows the average, the minimum, and the maximum of the ratio (2) for each combination of N, , and US. This table shows that the proposed algorithm is faster than the algorithm based on a binary search. This ratio increases in proportion to the increase in N, , and US. 8. Conclusion In this paper, we propose a discrete Newton method for the evacuation problem. First, we prove that the number of iterations of the proposed method is at most the number of vertices of the input dynamic network. Then we propose theoretical and practical implementation of the proposed method. It is important to speed up the proposed algorithms. For example, in Algorithm 4, the time-expanded network with respect to θi is a subgraph of the time-expanded network with respect to θi +1 . Thus, we can use the information obtained in the time-expanded network with respect to θi to compute a maximum static flow in the time-expanded network with respect to θi +1 . We can hope that this fact can speed up the proposed algorithms.

JID:TCS AID:12123 /FLA

Doctopic: Algorithms, automata, complexity and games

10

[m3G; v1.261; Prn:16/09/2019; 10:58] P.10 (1-10)

N. Kamiyama / Theoretical Computer Science ••• (••••) •••–•••

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The author would like to thank anonymous referees for helpful comments. This research was supported by JST, PRESTO Grant Number JPMJPR14E1, Japan. References [1] L.R. Ford Jr., D.R. Fulkerson, Constructing maximal dynamic flows from static flows, Oper. Res. 6 (3) (1958) 419–433. [2] L.R. Ford Jr., D.R. Fulkerson, Flows in Networks, Princeton University Press, 1962. [3] M. Skutella, An introduction to network flows over time, in: W. Cook, L. Lovasz, J. Vygen (Eds.), Research Trends in Combinatorial Optimization, Springer, 2009, pp. 451–482. [4] D. Dressler, G. Flötteröd, G. Lämmel, K. Nagel, M. Skutella, Optimal evacuation solutions for large-scale scenarios, in: Operations Research Proceedings 2010, 2010, pp. 239–244. [5] D. Dressler, M. Groß, J.-P. Kappmeier, T. Kelter, J. Kulbatzki, D. Plümpe, G. Schlechter, M. Schmidt, M. Skutella, S. Temme, On the use of network flow techniques for assigning evacuees to exits, Proc. Eng. 3 (2010) 205–215. [6] C. Even, V. Pillac, P.V. Hentenryck, Convergent plans for large-scale evacuations, in: Proceedings of the 29th AAAI Conference on Artificial Intelligence, 2015, pp. 1121–1127. [7] A. Takizawa, M. Inoue, N. Katoh, An emergency evacuation planning model using the universally quickest flow, Rev. Socionetwork Strategies 6 (1) (2012) 15–28. [8] H.W. Hamacher, S.A. Tjandra, Mathematical modelling of evacuation problems: a state of the art, in: M. Schreckenberg, S.D. Sharma (Eds.), Pedestrian and Evacuation Dynamics, Springer, 2002, pp. 227–266. [9] B. Hoppe, É. Tardos, The quickest transshipment problem, Math. Oper. Res. 25 (1) (2000) 36–62. [10] S. Iwata, L. Fleischer, S. Fujishige, A combinatorial strongly polynomial algorithm for minimizing submodular functions, J. ACM 48 (4) (2001) 761–777. [11] A. Schrijver, A combinatorial algorithm minimizing submodular functions in strongly polynomial time, J. Comb. Theory, Ser. B 80 (2) (2000) 346–355. [12] N. Megiddo, Combinatorial optimization with rational objective functions, Math. Oper. Res. 4 (4) (1979) 414–424. [13] M. Schlöter, M. Skutella, Fast and memory-efficient algorithms for evacuation problems, in: Proceedings of the 28th Annual ACM-SIAM Symposium on Discrete Algorithms, 2017, pp. 821–840. [14] W. Dinkelbach, On nonlinear fractional programming, Manag. Sci. 13 (1967) 492–498. [15] T. Radzik, Newton’s method for fractional combinatorial optimization, in: Proceedings of the 33rd Annual Symposium on Foundations of Computer Science, 1992, pp. 659–669. [16] M. Schlöter, Flows over Time and Submodular Function Minimization, Ph.D. thesis, Technische Universität Berlin, 2018. [17] R.K. Ahuja, T.L. Magnanti, J.B. Orlin, Network Flows: Theory, Algorithms, and Applications, Prentice-Hall, Inc., 1993. [18] J.B. Orlin, A faster strongly polynomial minimum cost flow algorithm, Oper. Res. 41 (2) (1993) 338–350. [19] D.M. Topkis, Minimizing a submodular function on a lattice, Oper. Res. 26 (2) (1978) 305–321. [20] N. Baumann, M. Skutella, Earliest arrival flows with multiple sources, Math. Oper. Res. 34 (2) (2009) 499–512. [21] M. Lin, P. Jaillet, On the quickest flow problem in dynamic networks – a parametric min-cost flow approach, in: Proceedings of the 26th Annual ACM-SIAM Symposium on Discrete Algorithms, 2015, pp. 1343–1356. [22] M. Saho, M. Shigeno, Cancel-and-tighten algorithm for quickest flow problems, Networks 69 (2) (2017) 179–188. [23] Y.T. Lee, A. Sidford, S.C. Wong, A faster cutting plane method and its implications for combinatorial and convex optimization, in: Proceedings of the 56th Annual Symposium on Foundations of Computer Science, 2015, pp. 1049–1065. [24] K. Murota, Discrete Convex Analysis, SIAM Monographs on Discrete Mathematics and Applications, vol. 10, Society for Industrial and Applied Mathematics, 2003. [25] J.B. Orlin, A faster strongly polynomial time algorithm for submodular function minimization, Math. Program. 118 (2) (2009) 237–251. [26] K. Nagano, A Faster Parametric Submodular Function Minimization Algorithm and Applications, Tech. Rep. METR 2007-43, The University of Tokyo, 2007. [27] L. Fleischer, S. Iwata, A push-relabel framework for submodular function minimization and applications to parametric optimization, Discrete Appl. Math. 131 (2) (2003) 311–322. [28] A.A. Hagberg, D.A. Schult, P.J. Swart, Exploring network structure, dynamics, and function using networkX, in: Proceedings of the 7th Python in Science Conference, 2008, pp. 11–15. [29] J.-C. Picard, M. Queyranne, On the structure of all minimum cuts in a network and applications, in: V.J. Rayward-Smith (Ed.), Combinatorial Optimization II, in: Mathematical Programming Studies, vol. 13, Springer Berlin Heidelberg, 1980, pp. 8–16. [30] H. Zheng, Y.-C. Chiu, P.B. Mirchandani, On the system optimum dynamic traffic assignment and earliest arrival flow problems, Transp. Sci. 49 (1) (2015) 13–27.