11institutetext:  Graduate School of Information Science, University of Hyogo,
Japan
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
#  Sink Location Problems in Dynamic Flow Grid Networks††thanks: Supported by
JSPS KAKENHI Grant Numbers 19H04068, 23H03349.

Yuya Higashikawa  Ayano Nishii  Junichi Teruyama  Yuki Tokuni

###### Abstract

A dynamic flow network consists of a directed graph, where nodes called
sources represent locations of evacuees, and nodes called sinks represent
locations of evacuation facilities. Each source and each sink are given supply
representing the number of evacuees and demand representing the maximum number
of acceptable evacuees, respectively. Each edge is given capacity and transit
time. Here, the capacity of an edge bounds the rate at which evacuees can
enter the edge per unit time, and the transit time represents the time which
evacuees take to travel across the edge. The evacuation completion time is the
minimum time at which each evacuees can arrive at one of the evacuation
facilities. Given a dynamic flow network without sinks, once sinks are located
on some nodes or edges, the evacuation completion time for this sink location
is determined. We then consider the problem of locating sinks to minimize the
evacuation completion time, called the sink location problem. The problems
have been given polynomial-time algorithms only for limited networks such as
paths [11, 2, 1], cycles [1], and trees [18, 10, 4], but no polynomial-time
algorithms are known for more complex network classes. In this paper, we prove
that the 1-sink location problem can be solved in polynomial-time when an
input network is a grid with uniform edge capacity and transit time.

###### Keywords:

facility location problem dynamic flow quickest transshipment problem
evacuation problem polynomial-time algorithm.

##  1 Introduction

In recent years, natural disasters such as earthquakes, tsunamis and
hurricanes have been occurring more frequently, and in Japan, the revision of
the Basic Act on Disaster Management [5] has accelerated the development of
evacuation facilities, i.e., locating tsunami evacuation towers and setting
tsunami evacuation buildings, which had been delayed. On the other hand, when
considering evacuation planning for urgent large-scale disasters such as
tsunamis, floods, and nuclear power plant accidents, one of the most important
issues to be considered is the delay in evacuation time due to traffic
congestion. In fact, it is known that many people were killed in the Great
East Japan Earthquake due to delayed evacuation caused by traffic congestion
[17].

To attack this issue, the dynamic flow network model proposed by Ford and
Fulkerson [7] can be applied. Since the dynamic flow network model can handle
the movement of people or objects over time, it is possible to develop an
evacuation plan that quantitatively takes traffic congestion into account. A
dynamic flow network consists of a directed graph, where nodes called sources
represent locations of evacuees, and nodes called sinks represent locations of
evacuation facilities. Each source and each sink are given supply representing
the number of evacuees and demand representing the maximum number of
acceptable evacuees, respectively. Each edge is given capacity and transit
time. The capacity of an edge bounds the rate at which evacuees can enter the
edge per unit time, and the transit time represents the time which evacuees
take to travel across the edge.

One of the most fundamental problems in dynamic flow networks is the quickest
transshipment problem. The objective of this problem is to compute the minimum
time by which each evacuees can arrive at one of sinks, i.e., the evacuation
completion time, and to find the optimal flow of evacuees that achieves the
evacuation completion time. For the quickest transshipment problem, several
strongly polynomial-time algorithms have been developed so far [16, 13, 19].

Given a dynamic flow network without sinks, once sinks are located on some
nodes or edges, the evacuation completion time for this sink location is
determined. We then consider the problem of locating sinks to minimize the
evacuation completion time, called the sink location problem. The problems
have been given polynomial-time algorithms only for limited networks such as
paths [11, 2, 1], cycles [1], and trees [18, 10, 4], but no polynomial-time
algorithms are known for more complex network classes. In this paper, we
address the sink location problem on grid networks. A grid network can model
an actual road network better than the networks studied so far, in that it
consists of a number of cycles. We present the first polynomial-time algorithm
for the 1-sink location problem in dynamic flow grid networks with uniform
edge capacity and transit time.

We describe the basic ideas of this study. In our model, a sink can be located
on a node or an edge. When a sink is located at a node, the evacuation
completion time for that sink can be computed in polynomial-time using the
algorithms [16, 13, 19] for the quickest transshipment problem. Therefore, if
one can compute in polynomial-time the optimal sink location for each edge,
i.e., the sink location that minimizes the evacuation completion time over the
points on the edge, the optimal sink location for the entire network can be
obtained in polynomial-time. In the following, we deal with the 1-sink
location problem on a particular edge of the input network.

This paper is organized as follows. In Section 2, we introduce notations,
models and known related results used throughout the paper. In Section 3, we
will give a polynomial-time algorithm for the 1-sink location problem on a
particular edge. Combining this algorithm and previous results, we will give a
polynomial-time algorithm for the 1-sink location problem on a grid network.
We conclude this paper in Section 4 with some discussions.

##  2 Preliminaries

###  2.1 Models

Let $\mathbb{R}$ and $\mathbb{R}_{+}$ denote the sets of real values and non-
negative real values. A dynamic flow network $\mathcal{N}$ is given as a
6-tuple $\mathcal{N}=(G=(V,E),S^{+},S^{-},w,c,\tau)$, where $G=(V,E)$ is a
directed graph with node set $V$ and edge set $E$, $S^{+}\subseteq{V}$ and
$S^{-}\subseteq V$ are sets of sources and sinks, respectively, function
$w\colon S^{+}\cup S^{-}\rightarrow\mathbb{R}$ represents supply of a source
or demand of a sink, function $c\colon E\rightarrow\mathbb{R}_{+}$ is capacity
of an edge, and function $\tau\colon E\rightarrow\mathbb{R}_{+}$ is transit
time of an edge. As for function $w$, for a source $s\in S^{+}$,
$w(s)(\geq{0})$ represents supply of source $s$, and for a sink $s\in S^{-}$,
$w(s)(<0)$ represents demand of sink $s$. We define
$w(X)\coloneqq\sum_{s\in{X}}w(s)$ for a subset $X\subseteq{S^{+}\cup S^{-}}$.

In this paper, we consider the problem of locating a sink in a dynamic flow
grid network. Here, let us describe an input network $\mathcal{N}$ in our
problem. First, a graph $G=(V,E)$ is a grid, where $V$ consists of $n=N\times
N$ nodes and $E$ consists of bidirected edges in between each adjacent nodes
of $V$ (see Fig. 2). Furthermore, we assume that all edge capacities and
transit times are uniform, thus $c$ and $\tau$ are constant functions. In the
sink location problem, input network $\mathcal{N}$ has no sinks and all nodes
are treated as sources. Therefore, $\mathcal{N}$ is represented as
$(G=(V,E),S^{+}=V,S^{-}=\emptyset,w,c,\tau)$. In the following, we abuse $c$
and $\tau$ as non-negative real constants.

Figure 1: The input network $\mathcal{N}$ |  Figure 2: $\mathcal{N}(y)$
---|---

Let us consider the problem of locating a sink $s^{-}$ in between two adjacent
nodes $p,q\in V$ at a distance of $y\ (0<y<\tau)$ from node $p$. For a given
dynamic flow grid network $\mathcal{N}$ without sinks, we apply the following
operations to $\mathcal{N}$:

  1. 1.

Remove directed edges $(p,q)$ and $(q,p)$.

  2. 2.

Add a sink $s^{-}$ with demand $w(s^{-})\\!=\\!-\\!\sum_{v\in{V}}\\!w(v)$.

  3. 3.

Add directed edges $(p,s^{-})$ with capacity $c$ and transit time $y$, and
$(q,s^{-})$ with capacity $c$ and transit time $\tau-y$.

The dynamic flow network obtained from the above operations is denoted as
$\mathcal{N}(y)$ (see Fig. 2). Hereafter, we denote the sets of nodes and
edges in $\mathcal{N}(y)$ as $V(y)$ and $E(y)$, respectively. In the
following, we consider embedding network $\mathcal{N}(y)$ in a Cartesian
coordinate system so that $p$ and $q$ are mapped at $(0,0)$ and $(\tau,0)$,
respectively, and the other nodes are mapped at points $(i\tau,j\tau)$ with
some integers $i,j$. For simplicity, a node mapped at $(i\tau,j\tau)$ is
called node $(i,j)$. We refer to the directions from $(0,0)$ to $(0,1)$, from
$(0,1)$ to $(0,0)$, from $(0,0)$ to $(1,0)$, and from $(1,0)$ to $(0,0)$ as
upward, downward, rightward, and leftward, respectively.

###  2.2 Evacuation completion time

A dynamic flow $f$ is defined as a function $f\colon
E(y)\times\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$, where $f(e,\theta)$
represents the flow rate entering edge $e\in E(y)$ at time
$\theta\in\mathbb{R}_{+}$. In this paper, we deal with dynamic flows in
continuous-time model [6, 12]. Let us consider the following constraints for a
dynamic flow $f$:

$\begin{split}0\leq f(e,\theta)\leq c\quad\text{for any}\ e\in E(y),\
\text{and for any}\ \theta\in\mathbb{R}_{+},\end{split}$ (1)
$\begin{split}\sum_{(v,u)\in E(y)}\int_{0}^{\theta}f((v,u),t)dt-\sum_{(u,v)\in
E(y)}\int_{0}^{\theta-\tau}f((u,v),t)dt\leq w(v)&\\\ \text{for any}\ v\in
S^{+}(=V(y)\setminus\\{s^{-}\\}),\ \text{and for any}\
\theta\in\mathbb{R}_{+}.&\end{split}$ (2)

The constraints (1) and (LABEL:eq:flow_conservation) are called the capacity
constraint and the conserve constraint. The conserve constraint
(LABEL:eq:flow_conservation) means that for any time $\theta$ and any source
$v$, the amount of flow out of $v$ within time $\theta$ is at most the amount
of flow entering $v$ within $\theta$ plus the amount of the supply at $v$.
Furthermore, for a time horizon $T\in\mathbb{R}_{+}$, consider the following
constraint:

$\begin{split}\sum_{(v,u)\in E(y)}\int_{0}^{T}f((v,u),t)dt-\sum_{(u,v)\in
E(y)}\int_{0}^{T-\tau}f((u,v),t)dt&=w(v)\\\ \text{for any}\ v\in
S^{+}(=V(y)\setminus\\{s^{-}\\}),&\\\
\int_{0}^{T-y}f((p,s^{-}),t)dt+\int_{0}^{T-(\tau-y)}f((q,s^{-}),t)dt&=-w(s^{-}).\end{split}$
(3)

The constraint (3) implies that for each node $v$, the net amount of flow
accumulated at $v$ within time $T$ equals its supply or demand. For a time
horizon $T$, if a dynamic flow $f$ on $\mathcal{N}(y)$ satisfies the above
constraints (1), (LABEL:eq:flow_conservation), and (3), then $f$ is said to be
feasible for $T$. Moreover, a time horizon $T$ for which a dynamic flow is
feasible is said to be feasible. The evacuation completion time denotes the
minimum time for which a feasible dynamic flow exists. Letting $\Theta^{*}(y)$
denote the evacuation completion time in $\mathcal{N}(y)$, the 1-sink location
problem between $p$ and $q$ is formulated as follows:

minimize $\displaystyle\Theta^{*}(y)$ subject to $\displaystyle 0<y<\tau,$

which we call Sink-Location-on-Edge (for short, SLE). Next, we describe known
properties of the evacuation completion time. Given a dynamic flow network
$\mathcal{N}(y)$, for a subset of the sources $X\subseteq S^{+}$ and a time
horizon $\theta\in\mathbb{R}_{+}$, let $o_{\theta}(X,y)$ be the maximum amount
of dynamic flow that can reach the sink $s^{-}$ from $X$ within time horizon
$\theta$ (more formal definition in our model will be given in (8)). The
following theorem by Hoppe and Tardos [13] says a property of a feasible time
horizon.

######  Theorem 2.1 ([13])

Given a dynamic flow network $\mathcal{N}(y)$ for a real value $y$ with
$0<y<\tau$, and a time horizon $\theta\in\mathbb{R}_{+}$, there exists a
feasible dynamic flow for $\theta$ on $\mathcal{N}(y)$ if and only if

$\min\\{o_{\theta}(X,y)-w(X)\mid X\subseteq{S^{+}}\\}\geq 0$ (4)

holds.

Therefore, the evacuation completion time $\Theta^{*}(y)$ in $\mathcal{N}(y)$
is the minimum $\theta$ that satisfies (4). Here, we define

$\Theta(X,y)\coloneqq\min\\{\theta\mid o_{\theta}(X,y)-w(X)\geq{0}\\}$ (5)

and then, because $o_{\theta}(X,y)-w(X)$ is monotonically non-decreasing in
$\theta\in\mathbb{R}_{+}$ for any $X\subseteq S^{+}$ and any $y$, we have

$\Theta^{*}(y)=\max\\{\Theta(X,y)\mid{X}\subseteq{S^{+}}\\}.$ (6)

###  2.3 Residual networks

The dynamic flow network $\mathcal{N}(y)$ can be treated as a static flow
network $\overline{\mathcal{N}}(y)$ by considering the transit time as the
cost of each edge. We say that a static flow $\bar{f}\colon
E(y)\rightarrow\mathbb{R}_{+}$ is feasible with respect to $X\subseteq S^{+}$
if $\bar{f}$ satisfies the following conditions:

$\displaystyle 0\leq\bar{f}(e)$ $\displaystyle\leq c$ $\displaystyle\text{for
any}\ e\in E(y),$ $\displaystyle\sum_{(v,u)\in
E(y)}\bar{f}(v,u)-\sum_{(u,v)\in E(y)}\bar{f}(u,v)$ $\displaystyle=0$
$\displaystyle\text{for any}\ v\in V(y)\setminus(X\cup\\{s^{-}\\}),$
$\displaystyle\sum_{v\in X}\sum_{(v,u)\in
E(y)}\bar{f}(v,u)-\sum_{u\in\\{p,q\\}}\bar{f}(u,s^{-})$ $\displaystyle=0.$

When a feasible static flow $\bar{f}\colon E(y)\rightarrow\mathbb{R}_{+}$ is
given to $\overline{\mathcal{N}}(y)$, the residual network with respect to
$\bar{f}$, denoted by $\overline{\mathcal{N}}(y)_{\bar{f}}$, is constructed as
follows. For each edge $e$ with $\bar{f}(e)>0$, a reverse edge
$\overset{\leftarrow}{e}$ is added, and the edge set obtained by adding such
reverse edges to $E(y)$ is denoted by $E(y)_{\bar{f}}$. The capacity of
residual edges $c_{\bar{f}}\colon E(y)_{\bar{f}}\rightarrow{\mathbb{R}_{+}}$
is defined so that $c_{\bar{f}}(e)=c-{\bar{f}}(e)$ for original edges $e\in
E(y)$, and $c_{\bar{f}}(\overset{\leftarrow}{e})={\bar{f}}(e)$ for reverse
edges $\overset{\leftarrow}{e}\in E(y)_{\bar{f}}\setminus E(y)$. The cost of
residual edges $\tau_{\bar{f}}\colon E(y)_{\bar{f}}\rightarrow{\mathbb{R}}$ is
defined as follows: For edges $(p,s^{-})$, $(s^{-},p)$, $(q,s^{-})$, and
$(s^{-},q)$, $\tau_{\bar{f}}(p,s^{-})=y$, $\tau_{\bar{f}}(s^{-},p)=-y$,
$\tau_{\bar{f}}(q,s^{-})=\tau-y$, and $\tau_{\bar{f}}(s^{-},q)=-\tau+y$
respectively. For other original edges $e\in E(y)$, $\tau_{\bar{f}}(e)=\tau$,
and for reverse edges $\overset{\leftarrow}{e}\in E(y)_{\bar{f}}\setminus
E(y)$, $\tau_{\bar{f}}(\overset{\leftarrow}{e})=-\tau$.

###  2.4 Envelope of two-dimensional line segments

Let $\mathcal{F}$ be a family of one-variable linear functions
$f_{1},\ldots,f_{n}$ defined on closed intervals $[a_{1},b_{1}],\ldots,$
$[a_{n},b_{n}]$, respectively. We simply refer to such $\mathcal{F}$ as a set
of two-dimensional line segments. Moreover, for each function
$f\in\mathcal{F}$ defined on a closed interval $[a,b]$, we define
$\underline{f}(x)$ as follows:

$\underline{f}(x)\coloneqq\left\\{\begin{array}[]{ll}\\!f(x)&(x\in[a,b]),\\\
\\!-\infty&(x\notin[a,b]).\end{array}\right.$

The upper envelope $U_{\mathcal{F}}$ of the two-dimensional line segments
$\mathcal{F}$ is defined as

$U_{\mathcal{F}}(x)\coloneqq\max_{1\leq{i}\leq{n}}\underline{f_{i}}(x)$

that is defined on $\bigcup_{1\leq i\leq n}[a_{i},b_{i}]$. A point
$(x^{\prime},U_{\mathcal{F}}(x^{\prime}))$ in the upper envelope
$U_{\mathcal{F}}$ is called a break point if a function constituting
$U_{\mathcal{F}}$ switches at $x^{\prime}$, that is, if

$\lim_{x\to
x^{\prime}+0}\mathrm{argmax}_{1\leq{i}\leq{n}}\underline{f_{i}}(x)\neq\lim_{x\to
x^{\prime}-0}\mathrm{argmax}_{1\leq{i}\leq{n}}\underline{f_{i}}(x).$

It is known that the following theorem holds for the upper envelope
$U_{\mathcal{F}}$ of the set of two-dimensional line segments $\mathcal{F}$.

######  Theorem 2.2 ([8, 9])

Given a set $\mathcal{F}$ containing $n$ two-dimensional line segments, the
upper envelope $U_{\mathcal{F}}$ has at most $O(n\alpha(n))$ break points, and
they can be computed in $O(n\log n)$ time. Here, the function $\alpha$ denotes
the inverse function of the Ackermann function.

In the upper envelope $U_{\mathcal{F}}$ of a two-dimensional line segments
$\mathcal{F}$, a line segment connects each pair of adjacent break points.
Then, this implies that some of break points gives the minimum value of
$U_{\mathcal{F}}$, and we have the following corollary.

###### Corollary 1

Given a two-dimensional line segments $\mathcal{F}$ containing $n$ line
segments, the break point that minimizes the upper envelope $U_{\mathcal{F}}$
can be computed in $O(n\log n)$ time.

##  3 Sink location on an edge

In this section, we propose a polynomial-time algorithm for problem SLE.
According to (6), if we can compute the function $\Theta(X,y)$ in $y$, for all
$X\subseteq S^{+}$, then the objective function $\Theta^{*}(y)$ can be
represented as the upper envelope of $\Theta(X,y)$ over $0<y<\tau$. However,
this brute-force method does not provide a polynomial-time algorithm since
there are exponentially many subsets $X$. On the contrary, we give a family of
$O(\sqrt{n})$ source sets that must contain $X$ maximizing $\Theta(X,y)$ for
any $y$ with $0<y<\tau$, and based on this property, provide a polynomial-time
algorithm. In Section 3.1, we show properties of $\Theta(X,y)$, in Section
3.2, we show properties of $X$ that maximizes $\Theta(X,y)$, and in Section
3.3, a polynomial-time algorithm is provided.

###  3.1 Properties of $\Theta(X,y)$

First, we introduce the notation used in the following. For any two sources
$u,v\in S^{+}$, we define $d(u,v)$ as the length of the shortest path from $u$
to $v$ in the static flow network $\overline{\mathcal{N}}(y)$. Note that the
shortest path from $u$ to $v$ does not use edges $(p,s^{-})$ and $(q,s^{-})$
in $E(y)$. For any source set $X\subseteq S^{+}$ and source $v\in S^{+}$, we
define $d(X,v)$ as the length of the shortest path from $X$ to $v$ in
$\overline{\mathcal{N}}(y)$, that is,

$d(X,v)\coloneqq\min\\{d(x,v)\mid x\in{X}\\}.$ (7)

Here, we discuss the properties of $o_{\theta}(X,y)$ which represents the
maximum amount of dynamic flow that can reach the sink $s^{-}$ from $X$ within
time horizon $\theta$. It is known that $o_{\theta}(X,y)$ is characterized by
a minimum-cost flow from $X$ to $S^{-}$ in the static flow network
$\overline{\mathcal{N}}(y)$ [7]. Specifically, given $X\subseteq S^{+}\cup
S^{-}$, $o_{\theta}(X,y)$ as a function in $\theta$ can be computed by
applying the successively shortest path algorithm [15, 3, 14] for the minimum-
cost flow problem in the following manner. Initially, set $\bar{f}$ as a zero
static flow, i.e., $\bar{f}:E(y)\rightarrow 0$. At each step $i\ (\geq 1)$,
execute the following two procedures.

  1. 1.

Find the shortest (i.e., minimum-cost) path $P_{i}(X,y)$ from $X$ to $s^{-}$
in the current residual network $\overline{\mathcal{N}}(y)_{\bar{f}}$. If
there is no such path, then break the iteration.

  2. 2.

Update $\bar{f}$ by adding a static flow of amount $c$ along path
$P_{i}(X,y)$.

Note that the number of paths obtained when the above iteration halts is
exactly two. Letting $|P_{1}(X,y)|$ and $|P_{2}(X,y)|$ be the lengths of paths
$P_{1}(X,y)$ and $P_{2}(X,y)$, respectively, $o_{\theta}(X,y)$ is represented
by the following equation.

$\begin{split}o_{\theta}(X,y)=\left\\{\begin{array}[]{ll}0&0\leq\theta<|P_{1}(X,y)|,\\\
c(\theta-|P_{1}(X,y)|)&|P_{1}(X,y)|\leq\theta<|P_{2}(X,y)|,\\\
c(\theta-|P_{1}(X,y)|)+c(\theta-|P_{2}(X,y)|)&|P_{2}(X,y)|\leq\theta.\end{array}\right.\end{split}$
(8)

Since (8) implies that $o_{\theta}(X,y)$ is a piecewise linear convex function
in $\theta$, it can be transformed into the following:

$\begin{split}o_{\theta}(X,y)=\max\left\\{0,c(\theta\\!-\\!|P_{1}(X,y)|),c(\theta\\!-\\!|P_{1}(X,y)|)+c(\theta\\!-\\!|P_{2}(X,y)|)\right\\}.\end{split}$
(9)

Note that $P_{1}(X,y)$ clearly does not contain any reverse edge since
$P_{1}(X,y)$ is a path in $\overline{\mathcal{N}}(y)$. In fact, not only that,
there exists $P_{2}(X,y)$ not containing any reverse edge, which follows Lemma
1 below. Thus, it is not necessary to consider the residual network when
calculating $o_{\theta}(X,y)$. Here, we define $d_{1}(X,y)$ as the minimum of
the shortest path lengths that pass through $p$ or $q$, respectively, from $X$
to $s^{-}$ on $\overline{\mathcal{N}}(y)$. We define $d_{2}(X,y)$ as the
maximum of these shortest path lengths. We thus have

$\begin{split}d_{1}(X,y)&\coloneqq\min\\{d(X,p)+y,d(X,q)+\tau-y\\},\text{
and}\\\ d_{2}(X,y)&\coloneqq\max\\{d(X,p)+y,d(X,q)+\tau-y\\}.\end{split}$ (10)

As for paths in $\overline{\mathcal{N}}(y)$ corresponding to $d_{1}(X,y)$ and
$d_{2}(X,y)$, we present the following lemma.

###### Lemma 1

For any source set $X\subseteq{S^{+}}$, there exists on
$\overline{\mathcal{N}}(y)$ a pair of edge-disjoint paths from $X$ to $y$
whose lengths are $d_{1}(X,y)$ and $d_{2}(X,y)$, respectively.

###### Proof

Figure 3: (a) The shortest path in $\mathcal{N}(y)$ from $p^{\prime}$ to $p$
that passes through $(i,0)$.
(b) Positions where reverse edges exist in the residual network
$\overline{\mathcal{N}}(y)_{\bar{f}}$.

First, for any pair of paths corresponding to $d_{1}(X,y)$ and $d_{2}(X,y)$,
if one path passes through $p$, then the other path passes through $q$. In
other words, if one path includes the directed edge $(p,s^{-})$, then the
other path includes the directed edge $(q,s^{-})$. Therefore, from equation
(10), it is sufficient to show the existence of a pair of shortest paths from
$X$ to $p$ and $X$ to $q$ in $\overline{\mathcal{N}}(y)$ that are edge-
disjoint. In the following discussion, we only consider the case where
$d(X,p)\leq d(X,q)$. Note that a similar argument holds even when
$d(X,p)>d(X,q)$.

Let $p^{\prime}=(i,j)\in X$ be the starting node of the shortest path from $X$
to $p$. Note that due to the assumption $d(X,p)\leq d(X,q)$, we have $i\leq
0$. Among shortest paths from $X$ to $p$, we choose the unique path passing
through node $(i,0)$, and let $\bar{f}$ be the static flow of amount $c$ along
this path (see Fig. 3(a)).

In the residual network $\overline{\mathcal{N}}(y)_{\bar{f}}$ (see Fig. 3(b)),
let $q^{\prime}=(i^{\prime},j^{\prime})\in X$ be the starting node of the
shortest path from $X$ to $q$. Let $P^{*}_{q}$ be the shortest path in
$\overline{\mathcal{N}}(y)_{\bar{f}}$ from
$q^{\prime}=(i^{\prime},j^{\prime})$ to $q=(1,0)$ via $(1,j^{\prime})$ that is
guaranteed not to contain any reverse edge, we thus have

$|P^{*}_{q}|=(|i^{\prime}-1|+|j^{\prime}|)\tau.$ (11)

In the following, we show that $P^{*}_{q}$ is one of shortest paths in
$\overline{\mathcal{N}}(y)_{\bar{f}}$ from $q^{\prime}$ to $q$, i.e., the
length of any other path is at least $(|i^{\prime}-1|+|j^{\prime}|)\tau$.

Let $P_{q}$ be any path from $q^{\prime}$ to $q$ in
$\overline{\mathcal{N}}(y)_{\bar{f}}$. Let us consider vertical and horizontal
lines passing through node $q$. We then call the half plane to the left (resp.
right) of the vertical line L-region (resp. R-region), and similarly the half
plane to above (resp. below) the horizontal line U-region (resp. D-region).
Here notice that $P_{q}$ consists of four types of edges: leftward, rightward,
upward and downward edges. Let $r_{1}$ (resp. $l_{1}$) be the number of
rightward (resp. leftward) edges of $P_{q}$ lying in R-region. Let $r_{2}$
(resp. $l_{2}$) be the number of rightward (resp. leftward) edges of $P_{q}$
lying in L-region. Let $u_{1}$ (resp. $d_{1}$) be the number of upward (resp.
downward) edges of $P_{q}$ lying in U-region. Let $u_{2}$ (resp. $d_{2}$) be
the number of upward (resp. downward) edges of $P_{q}$ lying in D-region. In
the residual network $\overline{\mathcal{N}}(y)_{\bar{f}}$, there are leftward
reverse edges from $p=(0,0)$ to $(i,0)$ in L-region and upward reverse edges
from $(i,0)$ to $p^{\prime}=(i,j)$ in U-region. Therefore, the number of
reverse edges in $P_{q}$ is at most $u_{1}+l_{2}$, we thus have

$|P_{q}|\geq(l_{1}-l_{2}+r_{1}+r_{2}-u_{1}+u_{2}+d_{1}+d_{2})\tau.$ (12)

First, we show that

$l_{1}-l_{2}+r_{1}+r_{2}=2r_{1}+|i^{\prime}-1|.$ (13)

When $i^{\prime}\leq 1$, we have

$l_{1}=r_{1},\quad l_{2}=r_{2}-|i^{\prime}-1|,$ (14)

and when $i^{\prime}>1$, we have

$\begin{split}l_{1}=r_{1}+|i^{\prime}-1|,\quad l_{2}=r_{2}.\end{split}$ (15)

In both cases of (14) and (15), equation (13) holds. Similarly, we show that

$-u_{1}+u_{2}+d_{1}+d_{2}=2d_{2}+|j^{\prime}|.$ (16)

When $j^{\prime}\leq 0$, we have

$u_{1}=d_{1},\quad u_{2}=d_{2}+|j^{\prime}|,$ (17)

and when $j^{\prime}>0$, we have

$u_{1}=d_{1}-|j^{\prime}|,\quad u_{2}=d_{2}.$ (18)

In both cases of (17) and (18), equation (16) holds. Therefore, by equations
(12), (13) and (16), we have

$\begin{split}|P_{q}|&\geq(2r_{1}+2d_{2}+|i^{\prime}-1|+|j^{\prime}|)\tau\\\
&\geq(|i^{\prime}-1|+|j^{\prime}|)\tau.\end{split}$ (19)

The last inequality is led by the fact that $r_{1}$ and $d_{2}$ are non-
negative. Therefore, by equations (19) and (11), we can conclude that
$P^{*}_{q}$ is the shortest path in $\overline{\mathcal{N}}(y)_{\bar{f}}$ from
$q^{\prime}$ to $q$. ∎

As mentioned above, by Lemma 1, we have $|P_{1}(X,y)|=d_{1}(X,y)$ and
$|P_{2}(X,y)|$ $=d_{2}(X,y)$. Substituting these into (9), $o_{\theta}(X,y)$
is represented by

$\begin{split}o_{\theta}(X,y)=\max\left\\{0,c\left(\theta-
d_{1}(X,y)\right),c\left(\theta-d_{1}(X,y)\right)+c\left(\theta-
d_{2}(X,y)\right)\right\\}.\end{split}$ (20)

By the definition of $d_{1}(X,y)$, we obtain

$\begin{split}c\left(\theta-
d_{1}(X,y)\right)&=c\left(\theta-\min\\{d(X,p)+y,d(X,q)+\tau-y\\}\right)\\\
&=\max\left\\{c(\theta-d(X,p)-\\!y),c(\theta-d(X,q)-\tau+y)\right\\}.\end{split}$
(21)

Since

$d_{1}(X,y)+d_{2}(X,y)=d(X,p)+d(X,q)+\tau$

holds by (10), we have

$\begin{split}c\left(\theta-d_{1}(X,y)\right)+c\left(\theta-
d_{2}(X,y)\right)=c\left(2\theta-d(X,p)-d(X,q)-\tau\right).\end{split}$ (22)

Substituting (21) and (22) into (20), we obtain

$\begin{split}o_{\theta}(X,y)=\max\bigl{\\{}&0,c(\theta-d(X,p)-y),c(\theta-d(X,q)-\tau+y),\\\
&c(2\theta-d(X,p)-d(X,q)-\tau)\bigr{\\}}.\end{split}$ (23)

Here, we describe the properties of $\Theta(X,y)$. According to its definition
(5), $\Theta(X,y)$ is the value of $\theta$ that satisfies
$o_{\theta}(X,y)=w(X)$. Regarding $\Theta(X,y)$ as a function in $y$, we have
the following theorem (see Fig. 4 for the shape of $\Theta(X,y)$ for
$0<y<\tau$).

Figure 4: The shape of $\Theta(X,y)$ for $0<y<\tau$ in the cases where (a)
$d(X,p)<d(X,q)$, (b) $d(X,p)>d(X,q)$, and (c) $d(X,p)=d(X,q)$, respectively.

###### Theorem 3.1

For any source set $X\subseteq{S^{+}}$, $\Theta(X,y)$ is a piecewise linear
function in $y$ with $0<y<\tau$ of at most three line segments, represented by

$\begin{split}\Theta(X,y)=\min\biggl{\\{}&y+\frac{w(X)}{c}+d(X,p),-y+\tau+\frac{w(X)}{c}+d(X,q),\\\
&\frac{w(X)}{2c}+\frac{d(X,p)+d(X,q)+\tau}{2}\biggr{\\}}.\end{split}$

###### Proof

Let us denote the linear functions in $\theta$ that compose $o_{\theta}(X,y)$
as

$\begin{split}g_{1}(\theta)&=c\left(\theta-d(X,p)-y\right),\\\
g_{2}(\theta)&=c\left(\theta-d(X,q)-\tau+y\right),\\\
g_{3}(\theta)&=c\left(2\theta-d(X,p)-d(X,q)-\tau\right).\end{split}$

By equations (5) and (23), $\Theta(X,y)$ is the minimum value of $\theta$
among values that make each of $g_{1}(\theta)$, $g_{2}(\theta)$, and
$g_{3}(\theta)$ equal to $w(X)$. By solving $g_{1}(\theta)=w(X)$ for $\theta$,
we obtain $\theta=y+\frac{w(X)}{c}+d(X,p)$. Similarly, $g_{2}(\theta)=w(X)$
when $\theta=-y+\tau+\frac{w(X)}{c}+d(X,q)$, and $g_{3}(\theta)=w(X)$ when
$\theta=\frac{w(X)}{2c}+\frac{d(X,p)+d(X,q)+\tau}{2}$. Therefore, the minimum
value among these three values is $\Theta(X,y)$. Additionally, with respect to
$y$, each of the above $\theta$ is linear function. Hence, $\Theta(X,y)$ is a
piecewise linear function in $y$ consisting of at most three segments. ∎

###  3.2 Dominant source sets

For two subsets of $S^{+}$, $X$ and $X^{\prime}$, $X$ is said to dominate
$X^{\prime}$ if $\Theta(X,y)\geq\Theta(X^{\prime},y)$ holds for all $y$. Let
us consider the sufficient condition for a source set $X$ to dominate another
source set $X^{\prime}$. If $d(X,p)=d(X^{\prime},p)$ and
$d(X,q)=d(X^{\prime},q)$, then functions $o_{\theta}(X,y)$ and
$o_{\theta}(X^{\prime},y)$ coincide according to (23). In this case, if
$w(X)\geq w(X^{\prime})$, then $\Theta(X,y)\geq\Theta(X^{\prime},y)$ by
definition (5). Among source sets $X$ such that $d(X,p)$ and $d(X,q)$ remain
the same respectively, one with the maximum $w(X)$ is called a dominant source
set. For two integers $i$ and $j$, let $X_{i,j}$ be the dominant source set
such that $d(X_{i,j},p)=i\tau$ and $d(X_{i,j},q)=j\tau$. Then, $X_{i,j}$ is
formally defined by the following:

$X_{i,j}\coloneqq\\{x\in V(y)\mid d(x,p)\geq i\tau,\ d(x,q)\geq j\tau\\}.$

Note that $X_{i,j}=\emptyset$ may hold when $i$ or/and $j$ are large enough.
Let $\mathcal{X}$ be the family of dominant source sets. Note that both
$d(X,p)$ and $d(X,q)$ are at most $2\tau\sqrt{n}$, we thus have

$\mathcal{X}\coloneqq\\{X_{i,j}\mid 0\leq i\leq 2\sqrt{n},0\leq j\leq
2\sqrt{n}\\}\setminus\\{\emptyset\\},$ (24)

and also

$\Theta^{*}(y)=\max\\{\Theta(X_{i,j},y)\mid X_{i,j}\in\mathcal{X}\\}.$ (25)

Although the number of source sets in $\mathcal{X}$ seems to be $O(n)$ by
definition (24), we can provide a more precise upper bound on the number of
$X_{i,j}$. The following lemma implies that when considering $X_{i,j}$ for a
fixed $i$, only the constant number of $j$ are enough to be considered.

###### Lemma 2

For any source set $X\subseteq{S^{+}}$, it holds in
$\overline{\mathcal{N}}(y)$

$d(X,p)-d(X,q)=\tau i,$

where $i$ is one of integers $-3,-2,-1,0,1,2,3$.

###### Proof

Figure 5: Division of the regions based on the inclusion of node $x$

Let us define subsets of node set in $\overline{\mathcal{N}}(y)$ as follows
(see Fig. 5):

$\begin{split}W_{1}&\coloneqq\\{(i,j)\in V(y)\mid i\leq 0,j\neq 0\\},\\\
W_{2}&\coloneqq\\{(i,j)\in V(y)\mid i\geq 1,j\neq 0\\},\\\
Z_{1}&\coloneqq\\{(i,j)\in V(y)\mid i\leq 0,j=0\\},\\\
Z_{2}&\coloneqq\\{(i,j)\in V(y)\mid i\geq 1,j=0\\}.\\\ \end{split}$

Furthermore, let $d(X,\\{p,q\\})$ denote the minimum value among $d(X,p)$ and
$d(X,q)$, that is,

$d(X,\\{p,q\\})=\min\\{d(X,p),d(X,q)\\}.$ (26)

Here, we know that there exists a node $x\in X$ such that

$d(X,\\{p,q\\})=\min\\{d(x,p),d(x,q)\\}.$ (27)

Let us consider the following four cases regarding the position of such $x$:
(Case 1) $x\in{W_{1}}$, (Case 2) $x\in{Z_{1}}$, (Case 3) $x\in{W_{2}}$, and
(Case 4) $x\in{Z_{2}}$. We give the proof only for Cases 1 and 2 because Cases
3 and 4 shown in a similar way.

(Case 1) When $x\in{W_{1}}$, we have

$d(x,q)=d(x,p)+\tau>d(x,p).$ (28)

Then by equation (26), (27) and (28), we have

$d(x,p)=d(X,\\{p,q\\})\leq d(X,p).$ (29)

Since $d(X,p)\leq d(x,p)$ holds, we have

$d(x,p)=d(X,\\{p,q\\})=d(X,p).$ (30)

Equation (30) implies that

$d(X,p)\leq d(X,q)$ (31)

holds. Furthermore, since $d(X,q)\leq d(x,q)$, we have

$d(X,q)\leq d(x,p)+\tau=d(X,p)+\tau.$ (32)

Considering equations (31) and (32), we obtain

$\begin{split}-\tau\leq d(X,p)-d(X,q)\leq 0.\end{split}$ (33)

(Case 2) When $x\in{Z_{1}}$, because we need to go through $W_{1}$ and $W_{2}$
from $x$ to $q$, we have

$d(x,q)=d(x,p)+3\tau>d(x,p).$ (34)

In a similar way to Case 1, we can conclude that equations (29) and (31) hold.
Since $d(X,q)\leq d(x,q)$, by equations (29) and (34), we have

$d(X,q)\leq d(X,p)+3\tau.$ (35)

Combining equations (31) and (35), we obtain

$\begin{split}-3\tau\leq d(X,p)-d(X,q)\leq 0.\end{split}$ (36)

(Case 3) When $x\in{W_{2}}$, similar to Case 1, we obtain

$\begin{split}0\leq d(X,p)-d(X,q)\leq\tau.\end{split}$ (37)

(Case 4) When $x\in{Z_{2}}$, similar to Case 2, we obtain

$\begin{split}0\leq d(X,p)-d(X,q)\leq 3\tau.\end{split}$ (38)

Summarizing the above cases, for all node sets $X\subseteq S^{+}$, we have

$-3\tau\leq d(X,p)-d(X,q)\leq 3\tau.$ (39)

Recall that $d(X,p)$ and $d(X,q)$ are integer multiples of $\tau$, which
concludes the proof. ∎

Recall that $d(X,p)$ can take at most $O(\sqrt{n})$ values, and for a fixed
value of $d(X,p)$, the number of values which $d(X,q)$ can take is at most
seven by Lemma 2. Therefore, we obtain the following theorem.

###### Theorem 3.2

The number of dominant source sets in $\mathcal{X}$ is $O(\sqrt{n})$.

###  3.3 Algorithms

Our algorithm consists of three steps. First, we compute $w(X_{i,j})$ for all
$X_{i,j}\in\mathcal{X}$. Next, we calculate the segments that compose
$\Theta(X_{i,j},y)$ for all $X_{i,j}\in\mathcal{X}$. Finally, we find the
value of $y$ that minimizes the upper envelope of the function family
$\\{\Theta(X_{i,j},y)\mid X_{i,j}\in\mathcal{X}\\}$.

Let us first describe a method to compute $w(X_{i,j})$ for all
$X_{i,j}\in\mathcal{X}$. Since the number of nodes is $n$, we can compute
$w(X_{i,j})$ in $O(n)$ time for any $i,j$. Moreover, we have
$|\mathcal{X}|=O(\sqrt{n})$ by Theorem 3.2, so all $w(X_{i,j})$ can be
computed in $O(n\sqrt{n})$ time. However, $w(X_{i,j})$ can be calculated more
efficiently by using dynamic programming. See the following lemma.

Figure 6: The regions of $X_{3,1}$ and $X_{4,1}$ consist of nodes within the
gray area respectively. The nodes located on the dashed line in (b) represent
the difference between $X_{3,1}$ and $X_{4,1}$.

###### Lemma 3

Values $w(X_{i,j})$ for all $X_{i,j}\in\mathcal{X}$ can be computed in $O(n)$
time.

###### Proof

Basically we calculate $w(X_{i,j})$ by subtracting the sum of supplies of
$O(\sqrt{n})$ nodes from $w(X_{i-1,j})$ or $w(X_{i,j-1})$ as shown in Fig. 6,
which can be done in $O(\sqrt{n})$ time. We then show the order of calculation
of $w(X_{i,j})$ as follows. First of all, we calculate $w(X_{0,0})$ in $O(n)$
time, and successively updating from $w(X_{0,0})$, we calculate $w(X_{1,0})$
and $w(X_{2,0})$. For $i=0,1,2$, successively updating from $w(X_{i,0})$, we
calculate $w(X_{i,1}),\ldots,w(X_{i,i+3})$. For $i\geq 3$, we update
$w(X_{i-1,i-3})$ to $w(X_{i,i-3})$, and successively updating from
$w(X_{i,i-3})$, calculate $w(X_{i,i-2}),\ldots,w(X_{i,i+3})$. Note that for
each $i$, we consider only at most seven subsets $X_{i,j}$ according to Lemma
2.

By Theorem 3.2, the total number of subsets $X_{i,j}$ considered is
$O(\sqrt{n})$, therefore, all values of $w(X_{i,j})$ can be computed in
$O(n)+O(\sqrt{n})\times O(\sqrt{n})=O(n)$ time. ∎

In the second step, we calculate the segments that compose $\Theta(X_{i,j},y)$
for all $X_{i,j}\in\mathcal{X}$. Recall that $d(X_{i,j},p)=i\tau$ and
$d(X_{i,j},q)=j\tau$. By Theorem 3.1, $\Theta(X_{i,j},y)$ is piecewise linear
function represented by

$\begin{split}\Theta(X_{i,j},y)=\min\biggl{\\{}&y+\frac{w(X_{i,j})}{c}+i\tau,-y+\tau+\frac{w(X_{i,j})}{c}+j\tau,\\\
&\frac{w(X_{i,j})}{2c}+\frac{i+j+1}{2}\tau\biggr{\\}}.\end{split}$ (40)

Using $w(X_{i,j})$ obtained at the first step, $\Theta(X_{i,j},y)$ can be
calculated in $O(1)$ time. According to equation (25), $\Theta^{*}(y)$ is the
upper envelope of the function family $\\{\Theta(X_{i,j},y)\mid
X_{i,j}\in\mathcal{X}\\}$. By equation (40), $\Theta(X_{i,j},y)$ consists of
at most three line segments, hence by Theorem 3.2, $\Theta^{*}(y)$ is
determined by the upper envelope of $O(\sqrt{n})$ line segments, which can be
computed in $O(\sqrt{n}\log n)$ time by Corollary 1. Thus, the value of $y$
that minimizes $\Theta^{*}(y)$ can also be found in $O(\sqrt{n}\log n)$ time
in the third step.

Summarizing the above discussions, the first step requires $O(n)$ time, the
second step requires $O(\sqrt{n})$ time, and the third step requires
$O(\sqrt{n}\log n)$ time. Therefore, we obtain the following theorem.

###### Theorem 3.3

Given a dynamic flow grid network of $N\times N$ nodes (without any sink) with
uniform edge capacity and transit time, the 1-sink location problem on a
particular edge can be solved in $O(n)$ time.

Solving problem SLE for each edge, the 1-sink location problem for all edges
can be solved in $O(n^{2})$ time. In general, the evacuation completion time
can be computed in $\tilde{O}(m^{2}k^{4}+m^{2}nk)$ time using the algorithm
developed by Schlöter et al. [19]. Here, $m$ is the number of edges and $k$ is
the total number of sources and sinks, and the $\tilde{O}$ omits the
logarithmic factor from Big $O$ notation. Adapting this algorithm to our
problem, the evacuation completion time for each node can be computed in
$\tilde{O}(n^{6})$ time since $m=O(n)$ and $k=O(n)$. If we focus on the sink
location on only nodes, by applying this algorithm to each node, the optimal
1-sink location can be computed in $\tilde{O}(n^{7})$ time. Therefore, we
obtain the following corollary.

###### Corollary 2

Given a dynamic flow grid network of $N\times N$ nodes (without any sink) with
uniform edge capacity and transit time, the 1-sink location problem can be
solved in $\tilde{O}(n^{7})$ time.

##  4 Conclusion

The sink location problem is a kind of evacuation problems that seeks to
determine the location of sinks that minimize the time it takes for each
evacuee to arrive at one of the evacuation facilities. In this paper, we
proposed a polynomial-time algorithm for the 1-sink location problem on a
dynamic flow grid network with uniform edge capacity and transit time. This is
the first polynomial-time algorithm for networks that contain a number of
cycles. We remark that our provided approach can be extended to the 1-sink
location problem on $M\times N$ grid networks, which gives a polynomial-time
algorithm even when $M\neq N$.

It would be interesting to develop polynomial-time algorithms for the sink
location problems on more complex networks, especially, a grid networks with
multiple number of capacities and transit times. In addition, from the
viewpoint of real world applications, developing polynomial-time algorithms
for grid networks with holes or the case of multiple sinks is a future
challenge.

## References

  * [1] Benkoczi, R., Bhattacharya, B., Higashikawa, Y., Kameda, T., Katoh, N., Teruyama, J.: Locating evacuation centers optimally in path and cycle networks. In: 21st Symposium on Algorithmic Approaches for Transportation Modelling, Optimization, and Systems (ATMOS 2021). Open Access Series in Informatics (OASIcs), vol. 96, pp. 13:1–13:19 (2021)
  * [2] Bhattacharya, B., Golin, M.J., Higashikawa, Y., Kameda, T., Katoh, N.: Improved algorithms for computing k-sink on dynamic flow path networks. In: Workshop on Algorithms and Data Structures. pp. 133–144. Springer (2017)
  * [3] Busacker, R.G., Gowen, P.J.: A procedure for determining a family of minimum-cost network flow patterns. Tech. rep., RESEARCH ANALYSIS CORP MCLEAN VA (1960)
  * [4] Chen, D., Golin, M.: Sink evacuation on trees with dynamic confluent flows. In: 27th International Symposium on Algorithms and Computation (ISAAC 2016). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik (2016)
  * [5] Disaster Management, C.O.: Basic act on disaster management (2017), retrieved June 30, 2023 from https://www.bousai.go.jp/taisaku/kihonhou/index.html
  * [6] Fleischer, L., Tardos, É.: Efficient continuous-time dynamic network flow algorithms. Operations Research Letters 23(3-5), 71–80 (1998)
  * [7] Ford Jr, L.R., Fulkerson, D.R.: Constructing maximal dynamic flows from static flows. Operations research 6(3), 419–433 (1958)
  * [8] Hart, S., Sharir, M.: Nonlinearity of davenport-schinzel sequences and of generalized path compression schemes. Combinatorica 6(2), 151–177 (1986)
  * [9] Hershberger, J.: Finding the upper envelope of n line segments in ${O}(n\log n)$ time. Information Processing Letters 33(4), 169–174 (1989)
  * [10] Higashikawa, Y., Golin, M.J., Katoh, N.: Minimax regret sink location problem in dynamic tree networks with uniform capacity. In: International Workshop on Algorithms and Computation. pp. 125–137. Springer (2014)
  * [11] Higashikawa, Y., Golin, M.J., Katoh, N.: Multiple sink location problems in dynamic path networks. Theoretical Computer Science 607, 2–15 (2015)
  * [12] Higashikawa, Y., Katoh, N.: A survey on facility location problems in dynamic flow networks. The Review of Socionetwork Strategies 13(2), 163–208 (2019)
  * [13] Hoppe, B., Tardos, É.: The quickest transshipment problem. Mathematics of Operations Research 25(1), 36–62 (2000)
  * [14] Iri, M.: A new method of solving transportation -network problems. Journal of the Operations Research Society of Japan 3(1), 2 (1960)
  * [15] Jewell, W.S.: Optimal flow through networks. In: Operations Research. vol. 6, pp. 633–633 (1958)
  * [16] Kamiyama, N.: Discrete newton methods for the evacuation problem. Theoretical Computer Science 795, 510–519 (2019)
  * [17] Kinjo, K., Matsumoto, M.: An analysis of the impact of individuals and society to evacuation based on the great east japan earthquake survey. Japan Society for Urbanology Annual Report 45, 104–112 (2011 (in Japanese))
  * [18] Mamada, S., Uno, T., Makino, K., Fujishige, S.: An ${O}(n\log_{2}n)$ algorithm for the optimal sink location problem in dynamic tree networks. Discrete Applied Mathematics 154(16), 2387–2401 (2006)
  * [19] Schlöter, M., Skutella, M., Van Tran, K.: A faster algorithm for quickest transshipments via an extended discrete newton method. In: Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA). pp. 90–102. SIAM (2022)