    # Minsum k-Sink Problem on Path Networks

We consider the problem of locating a set of k sinks on a path network with general edge capacities that minimizes the sum of the evacuation times of all evacuees. We first present an O(kn^4n) time algorithm when the edge capacities are non-uniform, where n is the number of vertices. We then present an O(kn^3 n) time algorithm when the edge capacities are uniform. We also present an O(n n) time algorithm for the special case where k=1 and the edge capacities are non-uniform.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Due to many recent disasters such as earthquakes, volcanic eruptions, hurricanes, and nuclear plant accidents, evacuation planning is getting increasing attention. The evacuation -sink problem is an attempt to model evacuation in such emergency situations cheng2013; hamacher2002. In this paper, a -sink means a set of sinks that minimizes the sum of the evacuation time of every evacuee to a sink.

Researchers have worked mainly on two objective functions. One is the evacuation completion time (minmax criterion), and the other is the sum of the evacuation times of all the evacuees (minsum criterion). It is assumed that all evacuees from a vertex evacuate to the same sink.

Mamada et al. mamada2006 solved the minmax 1-sink problem for tree networks in time under the condition that only a vertex can be a sink. When edge capacities are uniform, Higashikawa et al. higashikawa2014b and Bhattacharya and Kameda bhattacharya2015b presented time algorithms with a more relaxed condition that the sink can be on an edge. Chen and Golin chen2016b solved the minmax -sink problem on tree networks in time when the edge capacities are non-uniform. Regarding the minmax -sink on path networks, Higashikawa et al. higashikawa2015a present an algorithm to compute a -sink in time if the edge capacities are uniform. In the general edge capacity case, Arumugam et al. arumugam2014 showed that a -sink can be found in time. Bhattacharya et al. bhattacharya2017a recently improved these results to time in the uniform edge capacity case, and to time in the general case. Table 1 is a list of known algorithms and their time complexities.

The minsum objective function for the sink problems is motivated, among others, by the desire to minimize the transportation cost of evacuation or the total amount of psychological duress suffered by the evacuees. It is more difficult than the minmax variety because the objective cost function is not unimodal along a path, and, to the best of our knowledge, practically nothing is known about this problem on more general networks than path networks. A path network, although simple, can model an airplane aisle, a hall way in a building, a street, a highway, etc., to name a few. For the simplest case of and uniform edge capacities, Higashikawa et al. higashikawa2015a proposed an time algorithm. In Sec. 4 of this paper we present an time algorithm for the case of and general edge capacities. For the case of general and uniform edge capacities, Higashikawa et al. higashikawa2015a showed that a -sink can be found in time bounded by and . Bhattacharya et al. bhattacharya2018a recently showed that a minsum 1-sink in path networks with uniform edge capacities that achieves minmax regret kouvelis1997 can be computed in time.

A somewhat related problem, the quickest transshipment problem, is defined by dynamic flows with a given set of sources and sinks: each source has a fixed amount of supply, and each sink has a fixed demand. The problem is to send exactly the right amount of supply from each source to each sink with the minimum overall time. This problem has been studied for over fifty years since the seminal work by Ford and Fulkerson ford1958. The standard technique to solve this problem is to consider discrete time steps and make a copy of the original network for every time unit from time zero to a time horizon . This process produces a time-expanded network ford1958. Gale gale1959 posed the earliest arrival --flow problem, which is the problem of maximizing the amount of flow that reaches the single sink by any time , . See Wilkinson wilkinson1971, and Minieka minieka1973, Baumann and Skutella baumann2006 for more recent results.

The main contributions of this paper are and time algorithms for computing a minsum -sink, in the non-uniform and uniform edge capacity cases, respectively. These are improvements over our and time algorithms we presented at Iwoca 2018 benkoczi2018a. In achieving these results, we introduce two innovative methods. One is used in Sec. 5.2 to efficiently optimize functions formulated as Dynamic Programming (DP), and the other is used in Sec. 6.2 to compute the cost of moving the evacuees on a subpath to a potential sink. We also present an time algorithm for the special case where and the edge capacities are non-uniform.

This paper is organized as follows. In the next section, we define some terms that are used throughout this paper, and present a few basic facts. Sec. 3 introduces the concepts of a cluster and section (intuitively, a bunched group of moving evacuees), which play a key role in subsequent discussions. In Sec. 4, we design an algorithm which finds a minsum 1-sink. Sec. 5 formulates the framework for solving the minsum -sink problem, utilizing DP. In Sec. 6 we implement a solution method to the DP formulation, and analyze its complexity. Finally, Sec. LABEL:sec:conclusion concludes the paper.

## 2 Preliminaries

### 2.1 Model

Let denote a given path network, where the vertex set consists of , which we assume to be arranged in this order, from left to right horizontally. Vertex has weight (set of positive reals), representing the number of evacuees initially located at , and edge has length or distance and capacity , which is the upper limit on the flow rate through in persons/unit time. By , we mean that point lies on . For we write if lies to the left of . For two points , the subpath between them is denoted by , and (resp. ) denotes its length (resp. the minimum capacity of the edges on ). It takes each evacuee units of time to travel a unit distance on any edge.

Our model assumes that the evacuees at every vertex start evacuation at the same time at the rate limited by the capacity of its outgoing edge. We sometimes use the term “cost” to refer to the aggregate evacuation time of a group of evacuees to a certain destination. A -sink, which means a set of sinks, shares the following property of the median problem kariv1979b.

###### Lemma 1.

higashikawa2015a There is a -sink such that all the sinks are at vertices.

If we plot the arrival flow rate at, or departure flow rate from, a vertex as a function of time, it consists of a sequence of (temporal) clusters. The duration of a cluster is the length of time in which the flow rate corresponding to the cluster is greater than zero. A cluster consists of a sequence of sections, such that any adjacent pair of sections in it have different heights. In other words, a section is a maximal part of a cluster with the same height (= flow rate). A simple cluster consists of just one section. Clearly, in the uniform edge capacity case, all clusters are simple. A time interval of flow rate 0 between adjacent clusters is called a gap. Unless otherwise specified, we assume that evacuees arrive at vertex from vertex . The case where the evacuees move rightward can be treated symmetrically.

The head vertex of a cluster/section is the vertex from which the evacuee corresponding to the start time of the cluster/section originates. The offset of a cluster with respect to vertex is the time until the first evacuee belonging to the cluster arrives at . We say that a cluster/section carries (the evacuees from) vertex , if those evacuees provide flow to the cluster/section.

We define the following cost functions.

 ΦL(x) ≜ cost contribution to x from~{}P[v1,vi]~{}if % vi≺x⪯vi+1, ΦR(x) ≜ cost contribution to x from~{}P[vi+1,vn]~{}% if vi⪯x≺vi+1,

and

 Φ(x)=ΦL(x)+ΦR(x). (1)

A point that minimizes is called a minsum 1-sink.

The total cost is the sum of the costs of all the sections. The cost of a section of height with offset and duration is given by

 λt0+λ22c, (2)

where is the number of evacuees carried by the section higashikawa2015c. The average evacuation time for an evacuee carried by this section is , where represents the average delay before departure from the head vertex of the section, and the aggregate is given by , which yields (2). The first term in (2) represents the time for all evacuees carried by the section to travel from the head vertex of the section to the destination, which is away, and is called the extra cost of the section. The second term in (2) represents the time for all evacuees carried by the section to travel from their origin vertices to the head vertex of the section, and is called the intra cost of the section. To be exact, the ceiling function () must be used to compute costs, because the last group of evacuees to leave the head vertex may not occupy the full capacity of the outgoing edge. but we omit it for simplicity, and adopt (2) as our objective function cheng2013. Or we can consider each molecule of a fluid-like material as an “evacuee.”

A minsum k-sink partitions the path into subpaths, and places a 1-sink on each subpath in such a way that the cost of the max-cost 1-sink is minimized. We shall solve the easier 1-sink problem first and the -sink problem later.

## 3 Cluster/section sequence

In the 1-sink problem, to compute the intra and extra costs at , we obviously need to know the arrival section sequence at . Let (resp. ) denote the arrival section sequence at (resp. departure section sequence from) vertex from Right, both moving left. It is clear that , where denotes the empty sequence. Let us compute successively for . If the height of an arriving section at vertex is higher than , the evacuees carried by that section cannot depart from at the arrival rate. We see that the duration of the section gets stretched in this case, by the ceiling operation mamada2006. From now on, we use the verb ceil to mean performing a ceiling operation. Moreover, when the first evacuee of arrives at , there may be a backlog of delayed evacuees still waiting to depart from . We use those delayed evacuees to fill gaps in . This is only for convenience for analysis purposes, and we actually use the latest arrivals to fill gaps.

###### Proposition 1.

The order among the evacuees within a cluster can be rearranged arbitrarily without affecting the cost of the cluster.

###### Proof.

Consider evacuee (resp. ) who takes (resp. ) time to arrive at a destination. If we interchange the positions of these evacuees, the sum of their contributions to the total cost remains the same at .

If the height of an arriving section at is less than (the capacity of the exit edge from ), then we use the underutilized capacity to accommodate as many of the delayed evacuees as possible, together with the evacuees carried by the section. Based on Proposition 1, we fix the beginning of each arriving section at its original position in the time sequence.

###### Example 1.

Let consist of sections, . Fig. 1 illustrates a situation where the heights of some arriving sections are higher than . Figure 1: (a) αR(vi); (b) Amount equal to the light gray parts fill the gray parts; (c) Result.

Fig. 1(a) shows that the evacuees at leave at the rate of , forming section , and the first evacuee carried by the first section from arrives at before ends. ∎

The above example shows that the following two situations can arise, when is converted into .

#### 3.0.1 Observations

1. A stretched section, due to ceiling, ends in a gap. (Fig. 1(b) shows that the overlapping part of and , as well as the part of that exceeds capacity (shown in light gray) are merged into one section (renamed ), and get stretched in time, partially filling the gap before in Fig. 1(a).)

2. A section may shrink due to the stretched section preceding it, with its start time pushed to a later time. (In Fig. 1(b) and (c), the stretched “swallows” a part of and shrinks. The next section (such as in this example), if any, undergoes no change, since its height is less than .) ∎

## 4 1-sink problem

In this section, we first show how to compute the arrival and departure sequences efficiently. Based on them, we compute the extra and intra costs at each vertex, which provide the costs . Finally, a 1-sink is found by looking for a point that minimizes .

### 4.1 Computing section sequences

From Observations (a) and (b) above, we can easily infer the following lemmas.

###### Proposition 2.

For , the number of sections in is at most one more than that in .

###### Proposition 3.

For each , the heights of the sections in , due to the evacuees on subpath , are non-increasing with time.

We similarly define (resp. ), i.e., the arrival section sequence at (resp. departure section sequence from) from Left, due the evacuees on (resp. ). To compute and , we start from vertex . It is easy to construct , which consists of just a single section of height and duration that starts at time 0 (according to the local time at ), and , which is a gap (=offset) of length , followed by (according to the local time at ). Let us consider and for a general , , where . There are evacuees who depart from for at the rate of . Their departure takes time, and if or

 wi/diτ≤ci−1, (3)

then there is no interaction at between the cluster carrying the evacuees from and the first cluster in . In addition, if then the heights (=flow rates) of some clusters in get reduced and their durations become stretched, to become clusters in . Fig. 1 illustrates this situation. If , on the other hand, the first evacuee in arrives at when there is a backlog of evacuees from , still waiting for departure.

Let denote the total number of evacuees carried by cluster . For the cluster in , we record the value , where is the time difference between the start times of and the next succeeding section , which equals times the distance between the first vertices of and . In what follows, we omit the superscript from these quantities, since it will be obvious from the context. If there is no succeeding section after , we set . For each , we define the critical capacity

 λ(Cm)/tm. (4) Figure 2: The sum of the light gray areas equals the sum of the gray areas.

Fig. 2 is similar to Fig. 1, except that it does not show cluster and contains many symbols for easy reference. Our intention is to separate the ceiling operation from the gap-filling operation for the backlog at . In the figure, the start (resp. end) time of section is indicated by (resp. ). We thus have

 τm=m∑j=1tj+τ1.

The sum of the light gray areas above the capacity equals the sum of the gray areas, that end in the middle of section , which is partially “swallowed up” by the gray area. Let us now find the right end of the gray areas step by step. Since , we merge and , filling the gap between and , because the first evacuee carried by arrives at before the last evacuee carried by leaves .

Now and have been merged into one section. We repeat this process by remembering the ratio , and so forth. In general, we find the smallest111It is not unique without the “smallest” condition. So binary search cannot be used. such that satisfies

 m∑j=1λ(Cj)/m∑j=1tj≤ci−1. (5)

If (5) still holds when in it is replaced by , then the gray area ends at , where (as in Fig. 2), otherwise . In either case, it is easy to determine in (additional) constant time.

Once we determine the end of the gray area as above, we continue with the ceiling operations for subsequent sections in . To facilitate it, we introduce max-heap , in which we place critical capacities of (4), precomputed at the time is constructed. We can thus extract from the largest ratio . As long as contains a ratio larger than , we need to merge the corresponding sections.

###### Lemma 2.

We can compute section sequences in time.

###### Proof.

For we initialize max-heap , and we update it as we compute for . so that we can easily find the largest from in constant time and compare it with . If then we merge and . Since , the weight-time ratio for the combined section satisfies

 λ(Cm+1)/tm+1≤{λ(Cm)+λ(Cm+1)}/(tm+tm+1)≤λ(Cm)/tm.

We thus put in , which may or may not be the largest item in . Each insertion into heap takes time. Then we process the next largest item in and proceed with extending/merging. Proposition 2 implies that only new sections are created in total, as changes from to , each incurring constant processing time. Therefore, the total time for creating all new sections is , taking the insertion/extraction operations into/from into account.

We now go back to deal with the overlapping area (with backlog ) between sections and , shown in Fig. 1(a). Let , be the sections after ceiling by . To find the end time, , of the extended , we first determine, using binary search, the smallest that satisfies

 (t−τ1)ci−1≥w+l∑m=1λ(Cm), (6)

where . Once such a is found then it is straightforward to compute . This gap-filling operation takes time per vertex , and the total for all vertices is .222Sequential search is even faster, taking time in total. But this method is useful in solving the -sink problem.

### 4.2 Computing ΦR(vi) and ΦL(vi)

We discussed above how to convert into , which becomes with an offset. We want to compute cost for , where (resp. ) is the extra cost (resp. intra cost) of , defined after (2).333 can be computed similarly. It is clear that . It is also easy to see that and . Define an array of weights by , which can be precomputed in time for all . By definition, we have for ,

 Ei=Ei+1+diW[i+1]τ. (7)

Thus computing from takes constant time.

###### Lemma 3.

The extra costs, , can be computed in time.

Let us now consider intra cost , which is more difficult to compute. All the sections of the same height belonging to different clusters are said to form a height group, or just a group for short. They are always consecutive. Let consist of a set of groups of sections. The heights of the groups are decreasing from the first group to the last group in . Suppose that a group of height consists of sections , and, for , let denote the sum of the weights of the vertices carried by . Then the sum of the intra costs of the sections in is given by , and the total intra cost is . When sections in are merged due to a lower capacity , clearly changes for some groups .

During the construction of , we need to update the intra costs of the affected section groups. To this end, we add the squared weight of the newly stretched section and the squared weight of the part not swallowed up totally ( in Fig. 2). Then subtract the sum of the squared weights of the totally or partially swallowed up sections. We maintain with the sum of squared weight for each . Since a swallowed up section no longer contributes to the updating time, the total computation cost is . Thus the amortized updating time for the intra cost per section is constant.

Suppose that a group of sections have the same height , and their heights get reduced due to a new, smaller capacity . Let them have weights , so that the sum of their intra costs on departure from , i.e., in , is given by

 12ci−1{λ(S1)2+⋯+λ(Sg)2}=^λG2ci−1, (8)

where

 ^λG=λ(S1)2+⋯+λ(Sg)2. (9)

This implies that we can treat all the sections in as if they were just one virtual section with effective intra cost . This way, we can maintain just the sum of squared weights for . This idea works fine if sections do not merge with each other. Assume now that sections do merge due to a lower capacity encountered. It implies that sections disappear. We subtract the contributions of the disappeared sections from and replace them by the squared weight of the resulting large section. Then the new intra cost is given by . Since a section can disappear no more than once, the total computation cost of updating is . Clearly the intra costs of the groups that are not stretched do not change, so they incur no updating cost. The above argument proves the following lemma.

###### Lemma 4.

The intra costs, , can be computed in time.

Since , Lemmas 2, 3, and 4 directly imply the following corollary.

###### Corollary 1.

We can compute in time. Similarly, we can compute in time.

We then compute in additional time. A vertex that achieves is an optimal 1-sink. By Corollary 1, we have

###### Theorem 1.

A minsum 1-sink in path networks (with general edge capacities) can be found in time.

We formally state our algorithm as Algorithm 1 below.

## 5 k-sink problem

In this section, we solve the -sink problem, which entails applying a number of somewhat sophisticated methods.

### 5.1 DP formulation

We first present a dynamic programming (DP) formulation that follows the template of recursive functions proposed by Hassin and Tamir hassin1991 for the -median problem. Our innovation consists in the manner in which we process the recursive computations efficiently, given that the cost functions for the sink location problem are significantly more difficult to compute than those for the regular median problem, due to the time-dependent nature of the evacuee flow and congestion.

Let , , denote the minsum cost when sinks are placed on subpath . Similarly, define , , as the minsum cost when sinks are placed on subpath , and is the rightmost sink. We can start with , since for .

Let , denote the cost of evacuating from Right to all the evacuees on . Similarly, let , denote the cost of evacuating from Left to all the evacuees on . More generally, let us define and for . By definition, we have

 Fp(j) = minp≤i≤j{Gp(i)+R(i,j)}, (10) Gp(j) = minp≤i≤j{Fp−1(i)+L(i+1,j)}. (11)

To initialize the above recursive computations, we use , which can be computed from given by Corollary 1.444We could also start with , since , which can be obtained from Corollary 1. We can then compute , using (11), and , using (10), and so forth. We can thus compute , provided and are readily made available. Moreover, to obtain a DP algorithm with time complexity sub-quadratic in , we also need to quickly find the index that minimizes the recurrence relations (10) and (11).

Let us address the latter issue first. We refer to the evaluation of (10) and (11) as phase computation. To visualize our approach to computing (10), for a given phase , let us plot points in a 2-dimensional coordinate system for all , for a fixed vertex . See Fig. 3.

If we superimpose the line represented by for a given value in the same coordinate system, it is a line. If we increase from 0, this line eventually touches one of the plotted points. It is clear that the first point it touches gives the optimal value for that minimizes . In Fig. 3, this optimal is given by the point , hence . For convenience, let us refer to point as point .

We now explain that this representation provides us very useful information. To see it, for each point , define the -area that lies above the line and to the left of the vertical line through it as shown as shaded areas in Fig. 3. We say that a point situated in the -area of another point is dominated by , since the cost of point is higher than the cost of . In what follows, we identify a vertex with its index , and may say vertex to refer to . We sometimes say that is dominated by , when is clear from the context. Thus the points at the bottoms of the V-areas are the only non-dominated points. For subpath let , where and are the set of all points at the bottoms of the V-areas. As observed above, we have

###### Proposition 4.

For , holds.

Function can be computed in a similar manner, based on (11). Since , vertex is farther from than it is from , if . We thus have

 R(is,j+1)−R(is,j)≥R(it,j+1)−R(it,j)~{}~{} for~{}s

The upward arrows in Fig. 3 indicate the increase for different ’s. Moreover, if is dominated by , then point will also be dominated by for any . This implies that once it is determined that , then will not belong to for any . We will discuss how to update to , as increases, in the next subsection.

### 5.2 Computing switching points

To compute by Proposition 4, we maintain index set of non-dominated vertices. Let us denote by , , the switching point, namely the leftmost vertex , if any, for which dominates . If such an index does not exist, it means that never dominates and therefore we need not remember . Formally, we define, for ,

 x(is−1,is)=min{j′>is∣Gp(is)+R(is,j′)≤Gp(is−1)+R(is−1,j′)}. (13)

Define a sequence of switching points

 Xp[j]≜{⟨x(i1,i2),…,x(ig(j)−1,ig(j))⟩ if |Ip[j]|≥2Λ if |Ip[j]|=1. (14)

Maintaining together with allows us to update easily to .

###### Lemma 5.

The switching points in satisfy

 x(i1,i2)
###### Proof.

Assume for example that holds, where . Then will never be an optimal vertex, because for large enough which makes dominate , vertex already dominates , since . This implies that should not be in .

Procedure Add-Vertex() updates to for any phase . The while clause is justified by Lemma 5. algocf[htbp]     Let denote the time needed to compute for an arbitrary pair , .

###### Lemma 6.

Running Procedure Add-Vertex for all takes time.

###### Proof.

We need to compute in Line 6 of Add-Vertex(). It is computed once whenever a vertex is entered or removed from an index set, thus times in total.

Let denote the time required to compute for an arbitrary pair . Then we have

.

###### Proof.

To compute , we perform binary search for the pivot (in (13)) in the range . Each probe requires evaluating . By the above lemma, we are naturally interested in making as small as possible, which is the topic of Sec. 6.2.

## 6 Algorithms

### 6.1 Statement

Algorithm 2, given below, computes the minsum -sink. It requires a procedure to compute , which will be discussed in Sec. 6.2. It maintains index list, , and a list of switching points, , in each phase .

###### Lemma 8.

The minsum -sink in path networks with general edge capacities can be computed in plus preprocessing time.

###### Proof.

In each phase (), the most time consuming operation is the updating of and . Algorithm 2 performs iterations in Lines 2 and 4. The for loop (Line 7) invokes Add-Vertex() times. In these invocations, an element is added to at most once, so cannot receive more than a total of elements. Thus is computed times in each repetition of the outer for loop (Line 2), costing time by Lemmas 6 and 7. Similarly Line 9 costs per repetition of the outer for loop. The lemma follows from Lemma 6, since all other steps need less time.

### 6.2 Computing R(i,j) and L(i,j)

Line 6 of Algorithm 2 uses , based on Eq. (10), and Line 12 uses , based on Eq.(11). We only discuss below how to compute , which takes time by definition, since can be computed in similar time, .

#### 6.2.1 Preprocessing

As preprocessing before computing , we construct a binary tree , called the cluster sequence tree, whose leaves are the vertices of , arranged from to . For node555We use the term “node” for to distinguish them from the vertices of . The leaf nodes of are the vertices of . of , let (resp. ) denote the leftmost (resp. rightmost) vertex that belong to subtree . We say that spans subpath . Let (resp.