ON THE CAPACITY PROVISIONING ON DYNAMIC NETWORKS
OLUWASEUN FRANCIS LIJOKA
Bachelor of Science, Federal University of Technology, Akure, Nigeria, 2009
Master of Science, Kwame Nkrumah University of Technology/AIMS-Ghana,
Ghana, 2013
Ph.D. Mathematics, Heriot-Watt University United Kingdom, 2017
A thesis submitted
in partial fulfilment of the requirements for the degree of
DOCTOR OF PHILOSOPHY
in
THEORETICAL AND COMPUTATIONAL SCIENCE
Department of Mathematics and Computer Science
University of Lethbridge
LETHBRIDGE, ALBERTA, CANADA
© OLUWASEUN FRANCIS LIJOKA, 2022
ON THE CAPACITY PROVISIONING ON DYNAMIC NETWORKS
OLUWASEUN FRANCIS LIJOKA
Date of Defence:September 23, 2022
Dr. R. Benkoczi Professor Ph.D.
Thesis Supervisor
Dr. H. Shahadat Professor Ph.D.
Thesis Examination Committee Member
Dr. M. Khan Assistant Professor Ph.D.
Thesis Examination Committee Member
Dr. M. Tata Professor Ph.D.
Internal Examiner
Faculty of Arts and Science
Department of Neuroscience
Dr. M. Golin Professor Ph.D.
External Examiner
Hong Kong University of Science and
Technology, Hong Kong.
Dr. W. Osborn Associate Professor Ph.D.
Chair, Thesis Examination Committee
What we know is a drop, what we don’t know is an ocean
Sir Isaac Newton.
Dedication
I dedicate this project to God Almighty the creator of the universe and to my late mother
Mrs Teniola Grace Lijoka.
iv
Abstract
In this thesis, we consider the development of algorithms suitable for designing evacuation
procedures in sparse or remote communities. The works are extensions of sink location
problems on dynamic networks, which are motivated by real-life disaster events such as
the Tohoku Japanese Tsunami, the Australian wildfire and many more. The available algo-
rithms in this context consider the location of the sinks (safe-havens) with the assumptions
that the evacuation by foot is possible, which is reasonable when immediate evacuation
is needed in urban settings. However, for remote communities, emergency vehicles may
need to be dispatched or situated strategically for an efficient evacuation process. With
the assumption removed, our problems transform to the task of allocating capacities on
the edges of dynamic networks given a budget capacity c. We first of all consider this
problem on a dynamic path network of n vertices with the objective of minimizing the
completion time (minmax criterion) given that the position of the sink is known. This leads
to an O(n logn+ n log(c/ξ)) time, where ξ is a refinement or precision parameter for an
additional binary search in the worst case scenario. Next, we extend the problem to star
topologies. The case where the sink is located at the middle of the star network follows
the same approach for the path network. However, when the sink is located on a leaf node,
the problem becomes more complicated when the number of links (edges) exceeds three.
The second phase of this thesis focuses on allocating capacities on the edges of dynamic
path networks with the objective of minimizing the total evacuation time (minsum criterion)
given the position of the sink and the budget (fixed) capacity. In general, minsum problems
are more difficult than minmax problems in the context of sink location problems. Due to
few combinatorial properties discovered together with the possibility of changing objective
v
function configuration in the course of the optimization process, we consider the develop-
ment of numerical procedure which involves the use of sequential quadratic programming
(SQP). The sequential quadratic programming employed allows the specification of an ar-
bitrary initial capacities and also helps in monitoring the changing configuration of the
objective function. We propose to consider these problems on more complex topolgies
such as trees and general graph in future.
vi
Acknowledgments
I would like to acknowledge the contribution of individual that has made this work a success
in one way or the other. The greatest share of my gratitude goes to my supervisor Professor
Robert Benkoczi for being there always as a constant source of support, idea, courage and
encouragement. I also appreciate him for his patience and for being a good teacher whom
every student would dream to have. I can’t thank him enough for teaching me how to think
outside the box, analyse a problem, and eventually provide a solution to the problem. I
also appreciate his assistance towards my family in general as we are lucky to have Mrs
Benkoczi as our family doctor. I appreciate you a lot!
My next big appreciation goes to my committee members, Dr. Hossain Shahadat and
Dr. Muhammad Khan for their contributions in discussion, positive criticism and ideas. I
thank you both for always being there when I need to talk about my project and also during
presentation. You will always be part of my life. My utmost gratitude also goes to Professor
Wendy Osborn for teaching me advanced data analysis and for supporting all my applica-
tions for scholarships. You will always be part of my story. I will never forget to appreciate
our admin support, Barb Hodgson for her gentleness and help within the department. I
appreciate all the staff of math and computer science, including Professor Daya Gaou, Dr.
Jana Archibald, and Dr. John Sheriff for their support. I will not forget to appreciate Leila
Karimi and Sajad Fathi for organizing time for brainstorming and discussing ideas.
Next, I would like to appreciate my wife Margaret Ayoola for her support and under-
standing during the duration of the program. I really appreciate your drive and ambition to
move the family forward. I cherish you a lot. I am also grateful to Bolakunmi Banjo and
Lara Banjo for their support in taking care of my son Joshua and for always being there
vii
ACKNOWLEDGMENTS
when my family is in need of one help or the other. I am grateful to Olushoga Fashuyi for
supporting me financially and the family of Mrs Toyin Dada and Mrs Soboyejo for granting
us accommodation when we arrived in Lethbridge.
I gratuitously acknowledged the presence of some individuals that add fun to my life
during the program. In particular, my son Oluwasinaayomi Joshua who arrived just at the
beginning of my program. I would like to appreciate David Adebeshin, Fashuyi Ayodeji,
Dr. Felix Ighagbon, and all members of Canadian Red Cross in Lethbridge. I also want to
appreciate my cousin Dr. Akin Ojagbemi and my father inlaw Dr. Mathew Ayoola for their
support and encouragement. Finally, I appreciate the God Almighty who is the source of
all wisdom.
viii
Contents
Dedication iv
Abstract v
Acknowledgments vii
List of Tables xi
List of Figures xii
1 Introduction and Background 1
1.1 Concept of the minmax criterion with a fixed sink . . . . . . . . . . . . . . 3
1.2 Concept of minsum sink location problem with a fixed sink . . . . . . . . . 6
2 Background information and literature review 9
2.1 Overview of the sink location problem . . . . . . . . . . . . . . . . . . . . 9
2.2 Review of data structure for minmax evacuation criterion . . . . . . . . . . 14
2.2.1 The critical and upper envelope tree structure . . . . . . . . . . . . 15
2.3 Review of data structures for minsum evacuation criterion . . . . . . . . . 18
2.3.1 The 1-sink minsum problem . . . . . . . . . . . . . . . . . . . . . 19
2.4 Numerical approach to non-linear constrained problems . . . . . . . . . . . 22
2.4.1 The optimality conditions . . . . . . . . . . . . . . . . . . . . . . 25
2.4.2 Obtaining global optimum solution . . . . . . . . . . . . . . . . . 28
Convex programming . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.3 Unconstrained techniques for solving constrained problems . . . . 30
Sequential unconstrained minimization techniques (SUMT) . . . . 31
The multiplier methods . . . . . . . . . . . . . . . . . . . . . . . . 32
2.5 The sequential quadratic programming (SQP) method . . . . . . . . . . . . 33
2.5.1 SQP notations and model formulation . . . . . . . . . . . . . . . . 34
2.5.2 Quadratic subproblem formulation . . . . . . . . . . . . . . . . . . 36
2.5.3 The SQP algorithm framework . . . . . . . . . . . . . . . . . . . . 38
2.5.4 Newton version of the SQP . . . . . . . . . . . . . . . . . . . . . . 39
2.5.5 Other variants of SQP methods . . . . . . . . . . . . . . . . . . . . 41
2.5.6 Full Hessian methods . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5.7 The reduced Hessian method . . . . . . . . . . . . . . . . . . . . . 44
2.5.8 Merit function for SQP . . . . . . . . . . . . . . . . . . . . . . . . 46
ix
CONTENTS
3 Minmax capacity provisioning on path networks 50
3.1 Model formulation and problem definition . . . . . . . . . . . . . . . . . . 51
3.2 Properties of the optimal solution for an arbitrary path with the sink located
at an endpoint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3 Capacity provisioning for an arbitrary path network with the sink fixed at
an endpoint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.3.1 Implementation details and data structure . . . . . . . . . . . . . . 72
3.4 Capacity provisioning for evacuation of an arbitrary path with one arbitrary
sink fixed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4 Capacity provisioning on star networks 79
4.1 Star network with a sink in the middle . . . . . . . . . . . . . . . . . . . . 79
4.1.1 Algorithm for star networks with a sink in the middle . . . . . . . 81
4.2 Star networks with more edges and a sink located at the centre . . . . . . . 81
4.3 Three edged star network with the sink located on a leaf . . . . . . . . . . . 84
4.3.1 Algorithm strategy . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5 Minsum capacity provisioning on path networks 92
5.1 Preliminaries and model formulation for the minsum capacity provisioning
concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.2 Properties of the optimal solution vector and cost . . . . . . . . . . . . . . 103
5.3 Extension to longer path networks . . . . . . . . . . . . . . . . . . . . . . 105
5.4 Solution by sequential quadratic programming . . . . . . . . . . . . . . . 110
5.4.1 The SQP model and algorithm formulation . . . . . . . . . . . . . 111
5.4.2 The Hessian approximation . . . . . . . . . . . . . . . . . . . . . 115
5.4.3 Numerical experiment . . . . . . . . . . . . . . . . . . . . . . . . 117
6 Conclusion and future work 120
x
List of Tables
2.1 Table of minmax algorithms with their topologies and complexity in litera-
ture. [U] represents uniform edge capacities and [G] general edge capacities. 13
2.2 Table of minsum algorithms with their topologies and complexity in litera-
ture. [U] represents uniform edge capacities and [G] general edge capacities. 13
5.1 Table showing the solution for 2-edged network with fixed budgeted capac-
ity c = 20 and 4 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.2 Table showing the solution for 2-edged network with fixed budgeted capac-
ity c = 10. We reduce the budgeted capacity to let the initial z indicate a
congestion scenario. The optimal solution however indicate no congestion.
The algorithm takes 4 iterations to complete. . . . . . . . . . . . . . . . . . 117
5.3 Table showing the solution for 3-edge network with 5 iterations. . . . . . . 118
5.4 Table showing the solution for 4-edged network with c= 20. The algorithm
takes 8 iterations to complete. . . . . . . . . . . . . . . . . . . . . . . . . 119
5.5 Solution for the case where the sink is in the middle of the path. The so-
lutions are monotone from either sides. The algorithm takes 4 iterations to
complete. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
xi
List of Figures
1.1 Two-edge network with a sink fixed at the right end. . . . . . . . . . . . . . 3
1.2 A path network with a sink located on an edge. . . . . . . . . . . . . . . . 5
1.3 Time-weight graph modeling the objective function for the two-edged path
network of Figure 1.1 illustrating congestion free scenario. . . . . . . . . . 7
1.4 Time-weight graph modeling the objective function for the two-edged path
network of Figure 1.1 illustrating congestion scenario. . . . . . . . . . . . 8
2.1 Three sinks located on a path network. The picture shows how subpaths are
isolated with a sink fixed in each. . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Balanced tree constructed on the vertices. . . . . . . . . . . . . . . . . . . 16
3.1 Simple path network with sink in the middle. . . . . . . . . . . . . . . . . 52
3.2 Illustration of capacities allocation on a two-edged path network with the
sink located at the end vertex on the right. . . . . . . . . . . . . . . . . . . 55
3.3 Four-edged path network with sink at the middle. . . . . . . . . . . . . . . 59
3.4 Binary search procedure for the optimal capacity budget X1. . . . . . . . . 62
3.5 Illustration for the proof of part (ii) from Theorem 3.7. . . . . . . . . . . . 67
3.6 Example to illustrate longer path network with sink at the end. . . . . . . . 67
3.7 Graph illustrating the optimal computation of x2. . . . . . . . . . . . . . . 70
3.8 Explanation of how αi is computed. We let α′ , . . . ,α′1 t be the breakpoints of
yi−1 and p′0, . . . , p
′
t the corresponding linear pieces. Only the line segment
[αi,∞) is stored in the segment tree. . . . . . . . . . . . . . . . . . . . . . 71
3.9 Seven arbitrary y′i functions with their corresponding breakpoints. Only
one-sided line segments [αi,∞] are stored in the segment tree. . . . . . . . . 74
3.10 Stored line segments [αi,∞) for yi functions in Fig. 3.9. The vertical line
shows an arbitrary query point y = α′1 that cuts through only 3 line segments. 75
3.11 Illustration of the application of Algorithm 5 to an arbitrary path network
with a fixed sink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.1 Simple star graph with sink in the middle. . . . . . . . . . . . . . . . . . . 80
4.2 Regular star network with more edges and sink at the centre . . . . . . . . 82
4.3 Three-edged tree network with a sink at its leaf . . . . . . . . . . . . . . . 84
4.4 Graph illustrating the case when x0 ∈ ( c2 ,c) . . . . . . . . . . . . . . . . . 87
4.5 Graph illustrating the case when x ∈ ( c0 3 ,
c
2). Note that the slope of f12 is
not smaller than the larger of x1 and x2 but limited by x0. . . . . . . . . . . 88
4.6 Graph showing the critical weight WX at which θ2 = θ12. . . . . . . . . . . 90
5.1 Illustration of capacities allocation on a path network with two edges. The
sink is assumed to be located at the end vertex on the right. . . . . . . . . . 93
xii
LIST OF FIGURES
5.2 Time-weight graph of Figure 5.1. There is no congestion and evacuees get
to the sink with their starting rates . . . . . . . . . . . . . . . . . . . . . . 94
5.3 Congestion scenario for the path network in Figure 5.1. The intersection of
the two straight lines indicate the time at which θ2(x2,η) = θ1(x1,η) . . . . 96
5.4 The case when η = ηmin which corresponds to maximum value of x if the
optimal solution is in congestion scenario . . . . . . . . . . . . . . . . . . 101
5.5 Path network with three edges and sink located on the right end. . . . . . . 105
5.6 Congestion free scenario where η3,2 ≤ w3 and η2,1 ≤W [2,3] . . . . . . . . 106
5.7 Congestion only at v3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.8 Congestion only at v2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.9 Congestion at both vertices . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.10 Path network with two edges. . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.11 Path network with three edges. Fixed capacity at c = 20. . . . . . . . . . . 118
5.12 Path network with four edges. . . . . . . . . . . . . . . . . . . . . . . . . 118
5.13 Path network with four edges and a sink in the middle. . . . . . . . . . . . 119
xiii
Chapter 1
Introduction and Background
The importance of evacuation planning or strategy is becoming increasing in our world to-
day. Given the havoc and loss caused by recent natural disasters such as tsunami (in Tohoku
city, Japan), earthquakes (in Haiti), flood (in parts of England) and wildfire (in Alberta,
Canada) just to mention a few, governments of developed nations are particularly interested
in efficient measures or strategies that could aid the quick evacuation of people, livestock
and valuables immediately there is a sign or warning that demands such an evacuation.
The extent of destruction and loss of valuable properties and infrastructures may not be
quantifiable after the disaster. Bridges may be swept away, mansions may be demolished,
and train tracks wrecked just to mention a few. Above all, lack of an effective evacuation
procedure may result in a large death toll which can add to the aftermath psychological
trauma of the people and also cause a huge increase in the government disaster recovery
cost. If a natural disaster is unavoidable, as the name indicates, huge cost of recovery and
death toll are two negative impacts every government wants to avoid.
The required evacuation strategy in the event of emergency does not only involve time,
but also the identification of safe locations within a street, region or the whole city at large.
Depending on the scale of a disaster, evacuation strategy may also be planned within a
school, hospital environment, engineering environment and any working environment.
In this thesis, we focus on developing algorithms suitable for the effective evacuation in
remote or sparse communities. The existing algorithms in this context are developed with
the assumption that evacuees can get to the sink by foot, which may not be feasible when
1
1. INTRODUCTION AND BACKGROUND
considering remote communities as the sink may be located more than hundred kilometers
away. Therefore, emergency vehicles may have to be stationed strategically for an effec-
tive evacuation procedure. This leads to the problem of capacity provisioning on dynamic
networks, given a budgeted (fixed) capacity and with the assumption that the sink position
is known. Just as in traditional sink location problems, all algorithms considered in this
work are classified into two optimal criteria. The first with more flurry of work is the min-
imax sink location which focuses on evacuation completion time while the second called
the minsum sink location focuses on the sum of the evacuees’ evacuation times.
The minmax and minsum sink location problems can be seen as the generalization of
the NP−hard unweighted center and median problems respectively with respect to static
networks [24]. In fact, if the edge capacities are large enough, both sink location prob-
lems reduce to the respective static problems. Extensive work up to tree networks have
been considered for the case of the minimax criterion while only path networks have been
investigated in the case of the minsum criterion in literature.
For the rest of this thesis, we shall use the following definitions and parameters. Let
P = (V,E) represent a path graph with a vertex set V = {v1,v2, . . . ,vn+1} and an edge
set E = {ei = (vi,vi+1), i = 1, . . . ,n}. Each vertex vi has a non-negative weight wi which
represents the number of evacuees situated at vi, and each edge ei has a fixed length li but
unknown capacity xi, which represents the upper limit on the number of evacuees that can
be evacuated through the edge ei in a unit time. The parameter c denotes a given budgeted
capacity to be shared optimally between all the edge lengths and τ denotes the transit time
per unit distance.
We let p ∈ P denote any point p that lies on the path graph P, and for p,q ∈ P, p < q
(resp. p > q) is interpreted as point p lies at the left (resp. right) of point q. We define
d(p,q) as the sum of edge lengths between p and q and if p and/or q are on an edge, we
use prorated distance. We also denote the minimum capacity of the edges in a subpath
connecting points p and q by x(p,q). For evacuees to travel from p to q, τd(p,q) time is
2
1.1. CONCEPT OF THE MINMAX CRITERION WITH A FIXED SINK
expended.
1.1 Concept of the minmax criterion with a fixed sink
Let us review the concept of the minmax evacuation criterion, which simply minimizes
the evacuation maximum time or completion time in any topology of interest while also
considering the effect of congestion that may ensue as a result of the edge capacities. The
minmax problem has been studied up to tree topology and the approximation algorithm
approach to a general graph network has been proposed recently in literature [35]. We
review the concept by focusing on path networks with a fixed sink.
Problem 1.1. Given a path network P with a fixed sink (located anywhere on the path either
a vertex or on an edge) and known number of evacuees wi located at each vertex vi. Let us
assume that the edge length li and the edge capacity ci for each edge ei = (vi,vi+1), i ≥ 1
are available. Let also assume that the transit time τ to travel a unit distance is known. The
task is to minimize the maximum time of evacuation on the path network.
The above task is also equivalent to minimizing the evacuation time of the last evacuee
on the network. Let us look at the following path network for example, where evacuees
wi and w2 are located at vertices v1 and v2 respectively. The sink denoted by S is located
(fixed) on the right end.
w1 l w21 l2
v1 c1 v2 c2 S
Figure 1.1: Two-edge network with a sink fixed at the right end.
Now, the time for the last evacuees at v2 to get to the sink is
w
τl 22 + , (1.1)c2
where we refer to the first term as the distance time and the second term the waiting time.
For the evacuees moving from v1, we have to consider two cases:
3
1.1. CONCEPT OF THE MINMAX CRITERION WITH A FIXED SINK
Case 1 (No congestion): This means evacuees from v1 are not impeded at v2 as they journey
to the sink. Mathematically this implies τl1 ≥ w2c . With this assumption, we have the2
evacuation time to be
w
T 1v1 = + τ(l1 + l2), (1.2)min{c1,c2}
where we note that that the waiting time to the sink is dictated by the minimum of the edge
capacities.
Case 2 (Congestion): This is when evacuees from v1 have to stop at v2 due to backlog of
evacuees waiting to evacuate through the edge e2 to the sink. Mathematically this implies
τl < w21 c . In this case the evacuation time is given by2
w
T τl 1
+w2
v2 = 2 + . (1.3)c2
Before we present a more interesting example, we introduce the following array
i
W [i] =W [1, i] = ∑ w j, for i = 1,2, . . . ,n, (1.4)
j=1
so that for i ≤ j we have
W [i, j] =W [ j]−W [i−1]. (1.5)
Similarly we define the prefix sum of the edge lengths
i
L[i] = L[1, i] = ∑ l j, for i = 1,2, . . . ,n, (1.6)
j=1
so that for i ≤ j we have
L[i, j] = L[ j]−L[i−1]. (1.7)
Also we let c(i, j) denote the smallest capacity between subpath P[vi,v j]. Now let us con-
sider the following path network where the sink is fixed in the interior of the path
4
1.1. CONCEPT OF THE MINMAX CRITERION WITH A FIXED SINK
w1 = 15 w2 = 12 w3 = 10 w4 = 9 w5 = 7 w6 = 5 w7 = 11 w8 = 6 w9 = 5 w10 = 7
l1 = 4 l2 = 3 l3 = 5 l4 = 2 l5 = 3 l6 = 2 3 2 l8 = 5 l9 = 3
v1 v2 v3 v4 v5 v6 v7 S v8 v9 v10
c1 = 4 c2 = 5 c3 = 2 c4 = 6 c5 = 4 c6 = 3 c7 = 5 c8 = 3 c9 = 2
Figure 1.2: A path network with a sink located on an edge.
We have the left evacuation (where evacuees move from left to right) time with respect
to a vertex p ∈ [vi,v j] to the sink as
T [i, j]
W [i, p]
L (S, p) = τL[p,S]+ , for S > j (1.8)c(p,S)
and the right evacuation (where evacuees move from right to left) time is given as
T [i, j]
W [p, j]
R (S, p) = τL[S, p]+ , for S < i. (1.9)c(S, p)
As the sink is located on an edge say e j′ = (v j′,v j′+1) (resp. ei′ = (vi′,vi′+1)), a prorated
distance is used to compute the distance between the last vertex v′j and S (resp. the dis-
tance between the last vertex v′i and S). Let us compute the left evacuation time for the
subpath P v v of the network above. This means we need to compute T [2,7][ 2, 7] L (S, p) for
p = 2,3, . . . ,7. We have the following calculations based on equation (1.8)
12 22
T [2,7]L (S,2) = 18+ = 24 T
[2,7]
L (S,3) = 15+ = 262 2
[2,7] 31 38TL (S,4) = 10+ = 20 33 T
[2,7]
. L (S,5) = 8+ = 20.673 3
T [2,7]
43
S 6 5 19 33 T [2,7]
54
L ( , ) = + = . L (S,7) = 3+ = 13.83 5
From the above calculations, the maximum evacuation time for the subpath is 26 given by
T [2,7]L (S,3).
Definition 1.2. The vertex v ∈ v v that maximizes T [i, j]p [ i, j] L (S, p) is called the critical vertex
[5].
For example, the critical vertex for the subpath P[v2,v7] is v3. The example above shows
5
1.2. CONCEPT OF MINSUM SINK LOCATION PROBLEM WITH A FIXED SINK
that once the critical vertex within a subpath is determined, then the maximum evacuation
time for the subpath can be quickly computed. With n being the number of vertices in the
network, Bhattacharya et al. in [5] constructed a data structure called the capacities and
upper envelope (CUE) tree to aid the determination of the critical vertices in O(n log2 n)
time while Higashikawa et al. in [24] employed a dynamic programming data structure to
determine the critical vertex for the case when the edge capacities are uniform.
1.2 Concept of minsum sink location problem with a fixed sink
The minsum evacuation procedure focuses on minimizing the total evacuation time of
all evacuees on a specified network. This criterion is more challenging than the minmax
counterpart and all algorithms on the problem (in literature) are limited to path networks.
Let us review the concept with respect to a fixed sink on a path network. Given a path
network with a fixed sink and known number of evacuees wi located at each vertex vi.
Suppose also that the edge length li and edge capacity ci for each edge ei are known. The
task is to minimize the total evacuation time of all the evacuees to the sink. An important
fundamental property known in literature for this problem is given in Theorem 1.3 below
Theorem 1.3. The minsum sink must be on a vertex [24].
Let us consider the two-edged path network in Figure 1.1. We have two cases
Case 1 (No congestion): We construct the following time-weight graph to illustrate the
model or objective function for the case when τl ≥ w21 c .2
6
1.2. CONCEPT OF MINSUM SINK LOCATION PROBLEM WITH A FIXED SINK
T
f2(w)
Jump
f1(w) w2
τ l c22
0 w2 w1 + w2 W
Figure 1.3: Time-weight graph modeling the objective function for the two-edged path
network of Figure 1.1 illustrating congestion free scenario.
The x−axis is the prefix sum array of the weights while the y−axis is the time to the
sink. For example, the first evacuees from v2 get to the sink in time τl2 while the rest have
to wait w2c unit of time before they proceed to the sink. The graph also shows a gap which2
corresponds to the time of zero flow-rate and no congestion at v2. The areas of the two
trapezoids in the time-weight graph gives the required evacuation total time.
w2 w2
z = w2τl2 +w1τ(l1 + l2)+ 2 + 1 . (1.10)2c2 2min(c1,c2)
Note also that we can derive the model by considering the time for the last evacuees at both
v1 and v2 and computing the integral with respect to the evacuees. Thus we have
w
fv2(w) = τl2 + , w > 0c2
w−w (1.11)
fv1(w) = τ
2
(l1 + l2)+ , w > w ,min 2(c1,c2)
and ∫ w ∫2 w1+w2
z = fv2dw+ fv1dw
0 w2
2 2 (1.12)w w
= τl2w2 + τ(l 2 11 + l2)w1 + + .2c2 2min(c1,c2)
7
1.2. CONCEPT OF MINSUM SINK LOCATION PROBLEM WITH A FIXED SINK
Case 2 Congestion case: Let us assume that τl < w21 c , that is, evacuees from v1 are2
delayed at v2 due to some evacuees waiting to depart v1 to the sink. The following time-
weight graph illustrates the scenario:
T
f1(w)
f2(w)
τ(l1 + l2)
τ l2
w2 η w1 + w2
Figure 1.4: Time-weight graph modeling the objective function for the two-edged path
network of Figure 1.1 illustrating congestion scenario.
Using the diagram above and Equation (1.11), we have the objective function to be
1 1
z = ( fv2(0)+ fv2(η))η+ ( fv1(η)+ fv2(w1 +w2))(w1 +w2 −η), (1.13)2 2
where η is the prefix sum weight which includes the evacuees from v1 that get delayed at
v2. We can easily solved for η in the equation fv1(η) = fv2(η).
8
Chapter 2
Background information and literature
review
2.1 Overview of the sink location problem
The theory of network flow problems originated from a trivial notion of moving from
one place to another together with the objective of reducing some cost function. This theory
however opened a wide door for numerous interesting problems, where algorithm designers
often wish to design procedures reliable to transport supplies from fixed sources to some
target locations or sinks (not necessarily pre-established).
One of the well known applications of this theory is in the quickest transshipment prob-
lems. In this problem, a set of sources with specified amount of supplies and a set of sinks
with pre-determined demands are given. The objective is to send exactly the right amount
of supplies to the sinks in a shortest possible time. Ford and Fulkerson [29] in their quest
to tackle this problem (for one source and sink) introduced the concept of dynamic flow
network which has been exploited to tackle different evacuation problems. A dynamic
network is a graph (path, tree or general graph) which models the movement of supplies
(people, commodities, or other movable valuables) in a network. Each vertex in the graph
holds some amount of supply and each edge is equipped with non-negative capacity, which
limit the flow rate of supply through it. Also, the edges are often equipped with transit time
which represents the time of travel in a unit distance.
Hoppe and Tardos in [10] investigated the problem further for multiple sources and
sinks under the assumption of integral transit time. They presented the first polynomial
9
2.1. OVERVIEW OF THE SINK LOCATION PROBLEM
time algorithm for the problem. However, the complexity of their procedure is inapplicable
in reality and as far as we know, finding a practical polynomial algorithm for the quickest
transshipment problem is still open.
A related application of the dynamic flow theory is the evacuation problem (see survey
paper [25] ). In this problem, discrete supplies are given at some set sources and in contrast
to the quickest transshipment problem, the demands for the set of sinks are not specified.
The objective is to determine appropriate locations of the set of sinks in the network so as
to minimize the travelling time of all evacuees on the network to these sinks. Two optimal
criteria have been categorised in literature for the evacuation problem. The first known
as the minmax criterion focuses on locating sinks on the network so as to minimize the
evacuation completion time. Flurry of works have been published in this category even up
to tree topology. The second category which focuses on the total evacuation time of all
evacuees in the network is known as the minsum criterion. This category is harder than
the former and as a result, few publications limited only to path networks have appeared in
literature.
One important concept which characterized the flows in a network is the confluent or
non-confluent property. Whenever there is a restriction that all flows entering a vertex
must leave along same edge in the network, then we have confluent flows, otherwise non-
confluent flows. Confluent flows are appreciated in evacuation problems for orderliness
and to avoid confusion in the evacuation process. For both minmax and minsum optimal
criteria, it is difficult to construct an optimal confluent flow when the underlying network
is a general graph. In fact, if P ̸= NP, a constant-factor approximate scheme for confluent
flows cannot be achieved in polynomial time [39, 15]. The result of Hoppe and Tardos in
[10] and PTAS result of Belmonte et al in [35] are based on non-confluent flow property,
since there is no restriction on the flows in the network. Examples based on confluent flows
can be seen in (evacuation) sink location problems where the underlying networks are trees
and paths (see [12, 13, 24, 7] for examples.)
10
2.1. OVERVIEW OF THE SINK LOCATION PROBLEM
Under the minmax criterion with confluent flow, Mamada et al in [38] investigated the
location of 1-sink on a dynamic tree network with general edge capacities by employing a
table data structure. In particular, they constructed arriving and sending tables (at a specific
destination vertex) which denote the arrival and departure rate as functions of time. The end
of the arriving table clearly indicates the evacuation completion time. They assumed that the
sink can only be on a vertex and developed an O(n log2 n) time algorithm for the problem.
Also on the tree network, Higashikawa et al in [23] and Bhattacharya and Kameda in [7]
presented an O(n logn) time algorithms under the condition that the sink can be located
on any part of the network. They exploited the unimodality of the objective function and
the concept of tree’s centroid in developing their algorithms. The most recent multiple
sinks location algorithms under the tree topology was by Chen and Golin in [12] and [13].
They presented an O(max{k, logn}·kn log4 n) time algorithms with general edge capacities
while investigating the location of k-sink in the network. Systematically, their algorithm
isolates each subtree with enough sinks in a bottom up fashion until no subtree is left.
Under the non-confluent characterization of flows with minmax criterion, Belmonte et al
in [35] presented an approximation algorithm for locating sinks on general graph. Their
algorithm employed the polynomial-time algorithm of Hoppe and Tardos [10] as subroutine
to develop a fully polynomial time approximation scheme (FPTAS) [17, 32] for k-sink
location on general network. They further showed that the problem is W [1]-hard [17, 32]
when parametrized by k (the number of required sinks in the network). However, it is
known that the complexity of their underlying approach is practically inapplicable.
Considering the path networks, Higashikawa et al in [24] presented an O(kn) time al-
gorithm under the condition that the edge capacities are equal and the flow is confluent.
They formulated a recursive formulation for their approach and employed a 1-sink algo-
rithm to locate a sink in a subpath divided by a (k−1)-dimensional divider. Improving on
this was Bhattacharya et al in [5]. They presented O(min{n log3 n,n logn+ k2 log4 n}) and
O(min{n logn,n+ k2 log2 n}) time algorithms for general and uniform edge capacities re-
11
2.1. OVERVIEW OF THE SINK LOCATION PROBLEM
spectively. Their algorithms are based on the Megiddo’s parametric search [31] and sorted
matrix method of Frederickson and Johnson [18]. They employed a novel tree data structure
called the capacities and upper envelope tree (CUE tree) to facilitate their feasibility test as
well as efficient methods to extract information from the tree data structure. The minmax
k-sink location problem is a generalization of the NP-Hardness of the k-center problems in
general graph [27, 20].
Regarding the minsum sink location problems, only the path algorithms exist so far in
literature. The difficulty is due to the non-unimodal nature of the objective function. Hi-
gashikawa et al in [24] presented an O(n) time algorithm for the case of one sink on a path
√
network. They further showed that k−sink can be found in O(min{kn2,2O( logk log logn)n2}).
Benkoczi et al in [37] presented improved algorithms for the minsum sink location on the
same topology. They showed that 1−sink minsum problem on path with general edge ca-
pacities can be solved in O(n logn) time complexity. They further derived O(kn log3 n) and
O(kn log4 n) time algorithms for the k−sink location on path topology with uniform and
general edge capacities respectively. These algorithms are based on efficient optimized dy-
namic programming formulation and a tree data structure called the cluster sequence tree
which allowed efficient preprocessing of necessary information.
In comparison with the classical median problem, the minsum criterion can be reduced
to the classical unweighted median problem by making the edge capacities as large as possi-
ble. Hence minsum optimal criterion is a generalization of the NP-hardness of the classical
median problem [24, 37]. Tables 2.1 and 2.2 summarize the results known to date for the
minmax and minsum sink location problems.
12
2.1. OVERVIEW OF THE SINK LOCATION PROBLEM
Topology Problem Time complexity
1−sink [U] O(n) [24]
Path k−sink [U] O(kn) [24], O(n+ k2 log2 n) [5], O(n logn) [7]
k−sink [G] O(n logn+ k2 log4 n) [5], O(n log3 n) [5]
1−sink [U] O(n logn) [7, 23]
1−sink [G] O(n log2 n) [38]
Tree
k−sink [U] O(kn2 log4 n) [12], O(max{k, logn}kn log3 n) [13]
k−sink [G] O(kn2 log5 n) [12], O(max{k, logn}kn log4 n) [12]
1−sink [U] O(n) [4]
Cycle
1−sink [G] O(n logn) [4]
Table 2.1: Table of minmax algorithms with their topologies and complexity in literature.
[U] represents uniform edge capacities and [G] general edge capacities.
Topology Problem Time complexity
1−sink [U] O(n) [24]
√ √
k−sink [U] O(n2 ·min{ k logn+ logn,2 logk log logn}) [24], O(kn2 log3 n) [36]
Path
k−sink [U] O(kn log3 n) [37]
k−sink [G] O(kn2 log2 n) [36], O(kn log4 n) [37], O(kn2 log2 n)
Table 2.2: Table of minsum algorithms with their topologies and complexity in literature.
[U] represents uniform edge capacities and [G] general edge capacities.
All the novel algorithms discussed above emerged from the perspective of evacuating
densely populated or urban communities. However, remote areas or sparsely populated
communities, such as in British Columbia, Alberta, California and New South Wales are
not left out from disasters, especially wildfires and floods. These areas are far from urban
settings and evacuation by foot may not be efficient during unprecedented disasters, as
13
2.2. REVIEW OF DATA STRUCTURE FOR MINMAX EVACUATION CRITERION
the safehaven maybe located very far away. Therefore in our research, we shall consider
developing models and algorithms suitable for remote regions, where effective evacuation
may not be possible by foot and emergency vehicles may have to be deployed or situated
strategically to evacuate victims to safe haven (sink) hundreds of kilometers away.
2.2 Review of data structure for minmax evacuation criterion
The study of sink location problems present reliable procedures backed up by rigorous
theoretical analysis for city planners to locate safe-havens in strategic positions and to aid
the efficient evacuation process in the event of unforeseen disaster. In this section, we
review an efficient algorithm proposed in [6] for locating k sinks on a path network, in such
a way that the evacuation completion time is minimized. We acknowledge that Aramgutan
et al. in [3] had presented an O(kn log2 n) algorithm in 2014 while an O(kn) algorithm
was presented in [24] in 2015 for the case where the edge capacities are uniform, using
a special dynamic programming data structure. However, the data structure in [6] proved
more efficient as far as path network is concerned. As for tree networks, Mamada et al
in [30] presented an O(n log2 n) algorithm for 1-sink location with the assumption that the
sink is located on a vertex, using a special table data structure, while an O(n logn) time
algorithm for the case where the sink could be on any part of the tree was presented in [7].
With the important understanding of how crucial the critical vertices are to the compu-
tation of optimal cost (see definition 1.2) they construct an efficient data structure that can
speed up the determination of these vertices, and also to aid the important task called the
feasibility test. Given a time t, a problem instance is t−feasible if there exists a k−sink
such that every evacuees can evacuate to a sink in time not exceeding t. The figure below
shows how 3 sinks are located on a path network with 9 evacuee vertices.
14
2.2. REVIEW OF DATA STRUCTURE FOR MINMAX EVACUATION CRITERION
S1 S2 S3
w1 w2 w4 w6 w7
v1 v2 v3 v4 v5 v6 v7 v8 v9
Figure 2.1: Three sinks located on a path network. The picture shows how subpaths are
isolated with a sink fixed in each.
We generalize that k < n otherwise, we can assign a sink on each vertex to have optimal
cost of zero.
Let us describe the t−feasibility test before we explain the data structure. Starting with
initial vertex va, find vb such that all evacuees on the subpath P[va,vb] can be evacuated
within time t but evacuees on P[va,vb+1] fail the test. This is called the left feasibility test
and θL is used. Next locate a sink s on (vb,vb+1] such that θ
[a,b]
L (v
[a,b]
b,CL )+ xτ = t and
d(vb,s) = x. If s ∈ (vb,vb+1), set index c = b+ 1. If s = vb+1 then set c = b+ 2. Then
locate a vertex vd such that evacuees on P[vc,vd] can be evacuated to s within time t. This
we call the right feasibility test and θR is used. The feasibility test is repeated and at each
round, subpaths are isolated with a fixed sink until the last vertex vn belongs to the last or
kth isolated subpath.
2.2.1 The critical and upper envelope tree structure
To perform the t−feasibility test with respect to k sinks efficiently and to compute the
critical vertices as quickly as possible, they construct a balance tree on the path network
with the vertices being the leaves of the tree. Let v+j be the right end of vertex v j and let ul
and ur be the left and right children of node u of the tree T .
15
2.2. REVIEW OF DATA STRUCTURE FOR MINMAX EVACUATION CRITERION
ρ′
π(ρ′, vj)
π(ρ′, vi)
u
vi vl vu ru vj
Figure 2.2: Balanced tree constructed on the vertices.
During preprocessing, the minimum capacities as well as the cumulative weights asso-
ciated with any given subpath are stored at the nodes using a heap-like data structure bottom
up in O(n) time. We shall explain the left evacuation as the right evacuation procedure is
symmetric. We recall that for any given subpath P[vi,v j] the time to evacuate the evacuees
within the subpath to a sink s beyond v j is given by
{ }
θL = max d(vp,v j)τ+W [vi,vp]/c(vp,v j+1) (2.1)
vp∈[vi,v j]
Let us denote by π(u,u′) the path that connects the nodes u and u′ of the tree T and let ρ′
be the lowest common ancestor of any pair of leaves vi,v j. If we collect all the nodes that
are the right children of the nodes on the path π(vi,ρ′) and the left children of the nodes
on the path π(vi,ρ′) including vi and v j, (but ignoring the ones on the path) then we have a
total of O(logn) such nodes. Let N[vi,v j] be the set of the nodes and P (vi,v j) be the set of
O(logn) subpaths spanned by the nodes in N[vi,v j]. With the existence of N[vi,v j], formula
16
2.2. REVIEW OF DATA STRUCTURE FOR MINMAX EVACUATION CRITERION
(2.1) is broken into
{ }
max d(vp,vru)τ+(W +W [vlu,vp])/min{c(vp,vru+1),c} , (2.2)
vp∈V [vlu ,vru ]
where W =W [vi,vlu−1] and c = c(vru+1,v j+1). Equation (2.2) is further broken down into
two upper envelope functions with respect to each u ∈ N[vi,v j]. The first is a function of W
given by
u { }θL,1(W ) = max d{(vp,vr(u))τ+(W +W [vl(u),vp])/c(vp},vr(u)+1)vp∈V [vlu ,vru ] (2.3)
max θ[l(u),r(u)]
W
= L (v
+
r(u),vp)+ ,vp∈V [vl(u),vr(u)] c(vp,vr(u)+1)
and the second is a function of c given by
{ }
u W [vl(u),vθ τ p
]
L,2(c) = max d(vp,vr(u)) + . (2.4)
vp∈V [vl(u),vr(u)] c
Both formulas are derived by considering the smaller of c(vp,vr(u)+1) and c= c(vr(u)+1,v j).
They are piecewise linear and continuous in their respective variables and therefore can be
stored as sorted points using the idea of line duality and the computation of convex hull
achievable in linear time by Graham scan algorithm.
Lemma 2.1. Given a dynamic path network with n vertices, the capacity and upper enve-
lope tree T can be constructed with associated data in O(n logn) in time and space.
Once a pair of vertices (vi,v j) is given, the values of W and c are determined and the
critical vertices can be determined by binary search over the bending points of θuL,1(W )
and θuL,2(c), which is obtainable in O(logn) time. Thus, the L−critical vertices for P[vi,v j]
can be computed in O(log2 n) because we need to compare O(logn) candidates of N[vi,v j].
With the CUE tree prepared, it takes O(log3 n) time to locate a sink in an isolated sub-
path given a starting vertex va. This is because it takes O(log2 n) to compute the critical
17
2.3. REVIEW OF DATA STRUCTURES FOR MINSUM EVACUATION CRITERION
vertices as well as the left evacuation time and right evacuation time and then repeat the
computation O(logn) time. Finally for general edge capacities, the t−feasibility test to
insert k sinks on a path network can be carried out in O(min{n log2 n,k log3 n}) time. The
O(k log3 n) complexity can be achieved by locating 1 sink in the first partition of the path
network starting from v in O(log31 n) and then repeating the operation in k−1 more times.
The O(n log2 n) complexity is achievable by computing each left and right evacuation time
sequentially which takes O(n) and each computation takes O(log2 n) time.
When the edge capacities are uniform, the critical vertices can be computed by concate-
nating two subpaths in constant time and for each node in T bottom up, which takes O(n)
time to construct the CUE tree. Therefore when the CUE tree is available, the left and right
critical vertices for P[vi,v j], 1 ≤ i ≤ j ≤ n can be computed in O(logn) by exploring the
O(logn) nodes in N[vi,v j]. Finally the t−feasibility test for the case of uniform edge can
be carried out in O(min{n,k logn}) time provided the CUE tree is available.
2.3 Review of data structures for minsum evacuation criterion
We turn attention to review the problem of locating sinks on a path network in such a
way that the evacuation total time is minimized. In particular, we review the 1-sink minsum
work presented in [36] where special innovative method is employed, and for the k > 1 sinks
the innovative method is combined with dynamic programming to solve the problem effi-
ciently. We acknowledge that Higashikawa et al in [24] invented a dynamic programming
approach for the case of uniform edge capacity network. They later converted the prob-
lem to an equivalent problem of searching for k−edge path in a complete weighted directed
√ √
acyclic graph (DAG), an approach which yielded an O(n2 ·min{ k logn+logn,2 logk log logn})
time. When k = 1, their algorithm runs in O(n) time. The minsum sink location problem is
more difficult than the minmax sink location on any topology and its studies so far in liter-
ature is limited to path topologies. One important property of the minsum problem which
also differentiates it from the minmax is that the optimal location of the sink can only be on
18
2.3. REVIEW OF DATA STRUCTURES FOR MINSUM EVACUATION CRITERION
a vertex.
With reference to a particular vertex vi which could also be the location of the sink, if
the arrival flow rate at vi (or departure flow rate from vi) is plotted against time, then we
have a sequence of clusters. A cluster is regarded as the maximal time interval for which
the arrival flow rate at vi is continuously greater than zero. A cluster consists of sequence
of sections which may occur as a result of congestion at the arriving vertex. Sections in a
cluster have different heights, where the heights are given by the exiting edge capacities.
We have a simple cluster when the cluster contains only 1 section, for example, when the
edge capacity is uniform. The head vertex of a cluster is the vertex from which the evacuees
corresponding to the start time of the cluster originates. The head vertex contains the first
set of evacuees who are not impeded on their way to the reference vertex. Thus, the time
for which the first set of evacuees from a cluster reaches a reference vertex vi is called the
offset of the cluster. Finally, the time interval of flow rate zero between adjacent clusters is
termed a gap.
2.3.1 The 1-sink minsum problem
Let us consider the case when k = 1. We define the following time function (or cost)
θL(x) = Cost at x due to the vertices on P[v1,vi] if vi < x ≤ vi+1
(2.5)
θR(x) = Cost at x due to the vertices on P[vi+1,vn] if vi ≤ x < vi+1,
so that
θ(x) = θL(x)+θR(x). (2.6)
A point x ∈ V that minimizes θ(x) is called a minsum 1−sink. Thus, the total cost can be
computed by summing all the cost due to the vertices of all sections. Therefore, we need
a model to compute a cost due to a section. Given a section of height c with offset t0 and
19
2.3. REVIEW OF DATA STRUCTURES FOR MINSUM EVACUATION CRITERION
duration δt , the cost of the section is given by
λ2
λt0 + , (2.7)2c
where λ = cδt is the weight carried by the section. The first term of equation (2.7) is called
the extra cost, which is the time required for the evacuees carried by the section to move
from the vertex head to their destination. The second term is referred to as the intra-cost,
and it is the time required for all evacuees in the section to move from their origin vertices
to the head vertex.
Now, to solve the 1−sink minsum problem, we need an efficient way of detecting the
arrival section sequence denoted by αR(vi) and departure section sequence βR(vi) (both
assumed to be moving to the left). Computing αL(vi) and βL(vi) can be taken symmetrical.
We note that αR(vn) = βR(v1) = Λ, where Λ is an empty sequence. This is because
there is no right arrival at vn and no right departure from v1. We also note that if the height
of a section arriving vi is higher than the exiting edge capacity ci−1 of edge ei−1 then the
evacuees carried by that section cannot depart through vi, hence the duration of such section
gets stretched. There is also a possibility that an arriving section encounter a congestion at
vi due to some backlog of evacuees waiting to depart vi. This prompts the use of gap filling
operation where the delayed evacuees are used to fill the gaps in αR(vi). A section with
higher height than ci−1 is adjusted using the ceiling operation to stretch the duration of the
section. However, if the height of an arriving section is less than ci−1, the underutilized
capacity is used to accommodate as many delayed evacuees as possible. The following
observations are paramount to the algorithm:
(i) A stretched section may end in a gap
(ii) A section may shrink due to the stretched preceding it while its start time pushed to
a later time.
From the above observation, we realise that the number of sections in βR(vi) is at most one
20
2.3. REVIEW OF DATA STRUCTURES FOR MINSUM EVACUATION CRITERION
more than that in αR(vi) for i = n,n− 1, . . . . Moreover, since the evacuees coming from
further vertices would encounter smaller capacities, the heights of the sections in βR(vi)
due to the evacuees on subpath P[vi,vn] are non-increasing with time.
With the above properties at hand, we let λ(C) denote the number of evacuees carried
by a cluster C and for mth cluster in αR(vi), the ratio λ(Cm)/tm is recorded where tm is the
time difference between the start of Cm and the succeeding section Cm+1. We term the ratio
λ(Cm)
t the critical capacity.m
Now, clusters are merged depending on if there is an interaction (congestion) at the
departure vertices, using the filling and ceiling (if incoming height is greater than outgoing
height) operation. To also facilitate the computation of the sequences, max-heap H is
introduced to store the critical capacities. This lead to the following lemma
Lemma 2.2. The section sequences {αR(vi),βR(vi)|i= 1, . . . ,n} can be computed in O(n logn)
time.
For the efficient computation of the extra cost, we first of all create an array of weights
W [i] = ∑nj=i w j. Now using the formula
Ei = Ei+1 +diW [i+1]τ, (2.8)
for i = n−1,n−2, . . . ,1, we can compute the extra cost. In fact Ei−1 can be computed from
Ei in constant time. The weight array W [i] can be precomputed in linear time.
Lemma 2.3. The extra costs {Ei|1 ≤ i ≤ n} can be computed in O(n) time.
To compute the intra cost efficiently, we classify all sections of the same height in
different clusters as sections of the same group or height group. Suppose the arrival section
sequence αR(vi) consists of groups of sections, say G of which the heights of the groups
are non-increasing from the first to the last in G . Let G be a group of G with height h(G)
consisting of sections S1,S2, . . .Sg such that λ(Sq) is the sum of the weights carried by a
21
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
section Sq, 1 ≤ q ≤ g. Easily, we can compute the sum of the intra cost of the group G by
g
I(G) = ∑ λ(S 2q) /2h(G). (2.9)
q=1
Then the total intra cost can be computed using
I(G) = ∑ I(G). (2.10)
G∈G
However, changes occur in G when sections in αR(vi) are merged due to smaller capacity
ci−1 encountered along ei−1 and adjustments have to be done for accurate computation. To
do this, we add the square weight of the new stretched section and the square weight of the
section not totally swallowed, and we subtract the sum of the square weight of the totally or
partially swallowed up sections. This computation takes O(n) time because the swallowed
up section takes no part in the updating time. Based on the explanation above, we have the
following lemma
Lemma 2.4. The intra cost {Ii|1 ≤ i ≤ n} can be computed in O(n logn) time,
and all together we have the following theorem
Theorem 2.5. A minsum 1−sink in path network can be found in O(n logn) time.
2.4 Numerical approach to non-linear constrained problems
The general model for a constrained optimization can be written as follow
min f (x)
x∈Rn
s.t hi(xk) = 0, i ∈ E (2.11)
g j(x)≤ 0, j ∈ I ,
where E and I are sets of indices for equality and inequality constraints respectively. Let
us denote by p the size of E and by m, the size of I . We define the feasible region as
22
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
the set of points that satisfies all the constraints, and we denote by x⋆ the solution of the
problem. The solution vector x⋆ can be local or global. It is a global solution when f (x⋆) is
the lowest function value when compared to the values at all other feasible points. A local
solution vector is a set of points at which the objective function is smaller than the values of
all other feasible nearby points or neighborhood. Many algorithms especially, for nonlinear
problems often seek local solutions because the global solutions are often difficult to locate.
However, when the functions are convex, local solution x⋆ is also global.
Let us denote the feasible set by Ω given by
Ω = {x | hi(x = 0), i ∈ E ;g j(x)≤ 0, j ∈ I}, (2.12)
thus equation (2.11) can be rewritten compactly as
min f (x). (2.13)
x∈Ω
Optimality conditions are of two types. The first is the necessary conditions which must
be satisfied by any solution point (under certain assumptions) and the second is the sufficient
conditions which guarantee that x⋆ is in fact a solution. For unconstrained problems, the
optimality conditions are straight forward. The necessary condition is that the gradient of
the function at the optimal solution must be equal to zero, while the sufficient condition
is that the second derivative at the optimal solution must be positive definite. A solution
vector x⋆ is a local solution of problem (2.13) if x⋆ ∈ Ω and there exists a neighborhood
N of x⋆ such that f (x) ≥ f (x⋆) for all x ∈ N ∩Ω, x ̸= x⋆. The solution x⋆ is said to be
an isolated local solution if x⋆ ∈ Ω and there exists N of x⋆ such that x⋆ is the only local
solution in N ∩Ω.
One of the crucial properties for characterising solutions in constrained optimization
problems is the smoothness of the objective function and the constraints. The presence
of constraints in the problem definition often introduces many kinks and sharp edges and
23
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
that does not mean that the constraint functions that describe the feasible region are non-
smooth [33]. This is because a non-smooth constraint can be broken down into set of
smooth (piecewise) constraints. For example ∥x∥1 = |x1|+ |x2| ≤ 1 can be broken down
into x1+x2 ≤ 1, −x1+x2 ≤ 1, x1−x2 ≤ 1, and −x1−x2 ≤ 1. Also, a non-smooth objective
function can be reformulated as a smooth constrained problem.
One idea of dealing with the inequality constraints in the problem definition (2.11) is
transforming them to a set of equality constraints by adding a set of variables called the
slack variables [2]. Since the inequality constraints in our problem definition is of the
form ”≤”, their values are either negative or zero at a feasible point. Therefore, the slack
variables value must be either positive or zero at optimality. Thus, g j(x)≤ 0 can be written
as g j(x)+ s j = 0, s ≥ 0. The slack variables are then treated as part of the unknowns and
their values are determined as part of the solution set.
Note that when s j for a particular constraint is zero, then such constraint g j is active,
otherwise we have an inactive constraint. Note also that the size of the unknown variables
increases by this idea. To avoid additional constraints s j ≥ 0, we can add the square of the
slack variables to the constraint set to have g(x)+ s2 = 0.
Before we present the necessary and sufficient conditions for constrained problem (2.11),
we give the following definitions.
Definition 2.6. [33, 2, 8] The active set A(x) at any feasible x consists of the equality
constraint indices from E together with the indices of the inequality constraints j for which
g j(x) = 0. That is
A(x) = E ∪{ j ∈ I | g j(x) = 0}. (2.14)
At a feasible point x, the inequality constraints g j is said to be active if g j(x) = 0 and
inactive if otherwise.
Definition 2.7 (LICQ [33]). Given a point x and the active set A(x), the linear indepen-
dent constraint qualification (LICQ) holds if the set of the gradient of active constraints
24
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
{ }
∇g j(x), j ∈ A(x) is linearly independent. If LICQ holds, all active constraints gradients
are nonzero.
Definition 2.8 (Regular points). A point x⋆ satisfying the constraint h(x⋆) = 0 is said to
be a regular point of the feasible set if f (x⋆) is differentiable and gradient vectors of all
constraints at x⋆ are linearly independent.
Definition 2.9. Given a feasible point x and the active constraint set A(x), the set of lin-
earized feasible directions F (x) contains the vectors d that satisfy
dT ∇hi(x) = 0 ∀ i = 1, . . . , pF (x) = (2.15)dT ∇g j(x)≤ 0 ∀ j ∈ A(x)∩ I .
2.4.1 The optimality conditions
We review the necessary and sufficient conditions needed to be satisfied by x⋆ to be a
local solution of (2.11). The first order conditions are connected with the properties of the
gradients of the objective (in relation to unconstrained problems) and constraint functions
(extended to constrained problems). We have conditions in the following theorem.
Theorem 2.10 (First order necessary conditions [33, 2, 8]). Let x⋆ be a regular point of the
feasible set that is a local minimum for problem (2.11). Then there exists Lagrange multi-
pliers µ⋆ (a p-vector) and λ⋆ (an m-vector) such that the Lagrangian function is stationary
with respect to xk, k = 1, . . . ,n, µi i = 1, . . . , p, λ j j = 1, . . . ,m and s j, j = 1, . . . ,m at point
x⋆.
(i) The Lagrangian function for the problem is given by
p m
L(x,µ,λ,s) = f (x)+∑ µihi(x)+ ∑ λ j(g j(x)+ s2j)
i=1 j=1 (2.16)
= f (x)+µT h(x)+λT (g(x)+ s2)
25
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
(ii) Gradient conditions
∂L ∂ f p ∂h mi
= + µ⋆ +
∂x ∂x ∑ i ∂x ∑
⋆∂gµ ji = 0, k = 1, . . . ,n
k k i=1 k j=1 ∂xk
∂L
= 0 =⇒ hi(x⋆) = 0, i = 1, . . . , p (2.17)∂µi
∂L
= 0 =⇒ (g (x⋆)+ s2) = 0, j = 1, . . . ,m
∂λ i jj
(iii) Feasiblity check for inequality constraints. That is s2j ≥ 0 or g j ≤ 0 for j = 1, . . . ,m.
(iv) Switching conditions or the complementary slackness conditions:
∂L
= 0 =⇒ 2λ js j = 0, j = 1, . . . ,m. (2.18)∂s j
(v) Non-negativity of the Lagrange multiplier for the inequality constraints: λ⋆j ≥ 0 j =
1, . . . ,m.
(vi) Regularity check: that is the matrix
G = (∇h1,∇h2, . . . ,∇hp,∇g j1,∇g j2, . . . ,∇gqx) (2.19)
is linearly independent. The index qx is the size of the set of active constraints as-
sociated with the inequality constraints. If this condition holds, then the Lagrange
multipliers for the constraints are unique.
The switching conditions in the above theorem also tells us that if the inequality con-
straint g j(x) ≤ 0 is active at the candidate minimum point x⋆ then the corresponding La-
grange multiplier must be non-negative [2].
The first order conditions help us to understand the relationship between the derivatives
of the objective function and the active constraints at the optimal solution x⋆. When these
conditions are satisfied, a move along any vector d from F (x) either makes dT ∇ f (x⋆)> 0
26
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
or dT ∇ f (x⋆) = 0. When we have dT ∇ f (x⋆) = 0, the second order derivative informa-
tion helps us to understand whether a move along this direction increases or decreases the
objective function f .
To establish the second order conditions, we first impose a stronger smoothness assump-
tion on f , h and g. We require the functions to be twice continuously differentiable. Now,
given F (x) and some µ⋆ and λ⋆ that satisfy the KKT conditions, we define the critical cone
as
C (x⋆,µ⋆,λ⋆) = {d∈F (x) | ∇[hi g Tj] d= 0, i= 1, . . . , p and all j ∈A(x⋆)∩ I with λ⋆j > 0}.
(2.20)
Theorem 2.11 (Second order necessary conditions [2]). Suppose that x⋆ is a local solution
of (2.11) and that the LICQ holds. Let (λ⋆,µ⋆) be the Lagrangian multiplier vector for
which the KKT conditions hold. Then
dT ∇2L(x⋆,µ⋆,λ⋆)d ≥ 0, ∀ d ∈ C (x⋆,µ⋆,λ⋆). (2.21)
Theorem 2.12. (Second order sufficient condition [2]) Suppose that for some feasible point
x⋆ ∈ Rn there are Lagrange multiplier vectors λ⋆ and µ⋆ such that the KKT conditions hold.
Suppose also that dT ∇2L(x⋆,µ⋆,λ⋆)d > 0 for all d ∈ C (x⋆,µ⋆,λ⋆), d ̸= 0. Then x⋆ is a
strict local solution of problem (2.11).
The second order conditions can be reformulated in a slightly weaker version by using
two-sided projection of the Hessian of the Lagrangian HL(x⋆,µ⋆,λ⋆) onto the subspaces
that are related to C (x⋆,µ⋆,λ⋆). When the strict complementary conditions and LICQ hold,
we have
C (x⋆,µ⋆,λ⋆) = Null[∇[hx⋆, g(x⋆)]T ] j∈A(x⋆) = Null(G(x⋆)) (2.22)
where G is defined in (2.19). Let Z be a full column rank whose column span the space
27
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
C (x⋆,µ⋆,λ⋆), that is
⋆
C (x⋆,µ⋆,λ⋆) = {Zu | u ∈ R|A(x )|}, (2.23)
then the necessary condition can be written as
uT ZT HL(x⋆,µ⋆,λ⋆)Zu ≥ 0 ∀ u, (2.24)
or ZT HL(x⋆,µ⋆,λ⋆)Z is semi positive definite. Similarly, the second order sufficient the-
orem becomes ZT HL(x⋆,µ⋆,λ⋆)Z is positive definite. The matrix Z can be computed nu-
merically using the QR-factorization of the active constraints gradient.
2.4.2 Obtaining global optimum solution
Global solution to a constrained non-linear optimization problem can be obtained in
two ways [2, 33]:
(i) Exhaustive search method, which involves calculating all the local minimum points
and computing the cost function at those points. The point with the least value of cost
function is then taken as the global minimum. This approach can be employed if the
cost function is closed and continuous on a bounded feasible set
(ii) Convexity: the KKT conditions are sufficient to compute the global minimum points
if we can show that the optimization problem is convex.
We focus on the second method in this subsection.
Definition 2.13. A set S is said to be convex if two points p1 and p2 are any points in S
and the line segment joining the two points is entirely in S [2, 33].
Mathematically for two points x(1) and x(2) the line segment joining them is given by
x = αx(2)+(1−α)x(1), 0 ≤ α ≤ 1. (2.25)
28
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
If the whole line segment (2.25) is in S , then S is convex.
We can determine if a given function f (x) is convex in relation to a convex set.
Definition 2.14. A function f (x) is called convex on the convex set S if the graph of the
function lie below the line joining any two points on the curve f (x). For n−dimensional
space, f (x) defined on a convex set S is said to be convex if the following inequality is
satisfied [2]
f (α(x(2)+(1−x(1))))≤ α f (x(2))+(1−α) f (x(1)), 0 ≤ α ≤ 1. (2.26)
To check or determine the convexity of a function, equation (2.26) is not often used as
it may require an infinite number of pairs of points. The following theorem is handy to use
Theorem 2.15. A n−dimensional variable function f (x1,x2, . . . ,xn) defined on a convex
set S is convex if and only if the Hessian matrix of the function is positive semi-definite
or positive definite at all points in the set. The function f is strictly convex if the Hessian
matrix is positive definite for all points. The Hessian matrix of f is given by [2, 33]
 ∂2 f ∂2 f 2 ∂ 2 ∂x ∂x · · ·
∂ f
 x1 1 2 ∂x1∂xnH f =  ... ... · · · ...  (2.27)
∂2 f ∂2 f · · · ∂
2 f
∂x2n ∂xn∂x2 ∂x2n
The above theorem is necessary and sufficient condition to show if a function is convex.
Convex programming
For the constrained optimization problem, if function g j(x) is convex, then the singleton
set g j(x)≤ e j (e j a constant) is convex [2]. Now, if all g j for j = 1, . . . ,m are convex, then
the set {g j(x)≤ e j, j = 1, . . . ,m} is also convex, since the intersection of convex functions
is convex.
29
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
Theorem 2.16. Let Ss constitute the points that satisfy the constraint functions. That is
}
Ss = {x | hi(x) = 0, i = 1, . . . , p g j(x) = 0, j = 1, . . . ,m , (2.28)
then Ss is convex if functions set g are convex and functions h are linear [2, 33].
Feasible set Ss is also convex if the set is defined by linear equality or inequality func-
tions. Now if the cost function is convex over the feasible set, we have a convex program-
ming problem. The following theorem gives us the information needed for the local and
global solution for convex programming problems
Theorem 2.17. If f (x⋆) is a local minimum for a convex function f (x) defined on a convex
feasible set S , then f (x⋆) is also a global minimum and x⋆ is the global minimum vector
[2, 33].
2.4.3 Unconstrained techniques for solving constrained problems
The study of unconstrained non-linear optimization problems has opened a bank of
techniques that can be exploited to solve constrained optimization problems. Many of the
methods include the Newton method, which offers a second-order approximation and the
quasi-Newton methods, which combine the attractive features of the steepest descent and
Newton method. To employ the unconstrained methods to solve constrained problems,
some sort of transformations need to be done. The objective function of the constrained
problem is often wedded together with its associated constraint functions to form a compos-
ite function. Penalty parameters are often incorporated in the composite function formed,
and are included to penalise the new objective function if the constraints are violated. The
values of the penalty parameters are directly proportional to the extent of the violation. That
is, the bigger the violation, the larger the penalty value [2, 33].
When the composite objective function is constructed successfully with some set of
penalty parameters, it can be minimized using any of the unconstrained methods. At every
30
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
iteration, the penalty parameters are adjusted inline with some conditions and the composite
objective function is again redefined and minimized. The process continues until desired
solutions are obtained.
We review the transformation method under two categories [2, 33]:
(i) Sequential unconstrained minimization technique (SUMT)
(ii) Multiplier method.
Both transformation categories transform the problem to the following form [2]
Φ(x;r) = f (x)+P(h(x),g(x),r), (2.29)
where r denotes a vector of penalty parameters and P is a real-valued function that imposes
a penalty on the cost function with the help of r. We discuss the two categories in turn.
Sequential unconstrained minimization techniques (SUMT)
Under this categories are the penalty function method and barrier function method. The
penalty function method constructs the real-valued function P in (2.29) in such a way that a
positive value of P is added to f (x) whenever there is a violation of constraint. An example
is the quadratic loss function given as [33]
p m
P(h(x),g(x),r) = r{∑[h 2i(x)] + ∑ [g+i (x)]2}
i=1 j=1 (2.30)
g+i (x) = max(0,gi(x)),
where r is a scalar penalty parameter. It is easy to observe that the function P is always
positive in value. We also observe that g+i (x) = 0 if the inequality constraint is active
(gi(x) = 0) or inactive (gi(x)< 0) but positive if the inequality constraint is violated.
The advantages of the method include the ability to search through infeasible region and
the possibility of assigning an arbitrary initial solutions for both the solutions and penalty
31
2.4. NUMERICAL APPROACH TO NON-LINEAR CONSTRAINED PROBLEMS
parameter. The drawbacks of the method include the fact that the iterative process may stop
prematurely and the final solution may not be feasible if this happens.
The barrier function methods designed for inequality-constrained problems include [33,
2]
(i) The inverse barrier function
1 m −1
P(g(x),r) = ∑ (2.31)r j=1 gi(x)
(ii) The Log barrier function
1 m
P(g(x),r) =
r ∑ log(−gi(x)). (2.32)j=1
Both functions above give an infinite value whenever any of the inequality function is active.
Therefore it is not advisable to start the method with some feasible solutions due to large
barrier set around the feasible set. We can observe that as r → ∞, x(r)→ x⋆. Some of the
weaknesses of the barrier methods include instability (ill-behavior) near the boundary of the
feasible set, unclear method of choosing or selecting the iterates of the penalty parameters
and the ill-conditioning of the Hessian of the unconstrained function as r → ∞
The multiplier methods
The methods under this category provide ways of circumventing the difficulties of the
methods under the SUMT category. To make sure that the composite function Φ has a good
conditioning, the transformation is carried out in such a way that the penalty parameters
need not go to infinity. The methods also have the advantage of faster convergence rate
compared to the SUMT methods with an arbitrary starting solution. The multiplier methods
are also referred to as the augmented Lagrangian methods in literature. Below is an example
of the method, whose construction introduces a penalty and multiplier for each constraint
32
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
[2]
1 p m
P(h(x),g(x),r,θ) = ∑ r′(h +θ′ 1i i i)2 + r ((g +2 2 ∑ j j +θ j) ), (2.33)i=1 j=1
where θ j > 0, r > 0, and θ′j i, r′i > 0 are parameters associated with the jth inequality and
ith equality constraints respectively. Also, g( +i+θi) = max(0,gi+θi). It is easy to observe
that if θ′i = θ j = 0 and r′i = r′j = r, then we have the quadratic loss function (2.30) which
requires r → ∞ in order to enforce convergence.
For a special case of equality-constrained problem, the augmented or composite func-
tion can be defined as
p 1
ΘE(x,h(x),r) = f (x)+∑(µih (x)+ rh2i (x)), (2.34)
i=1 2
i
where r > 0 is a penalty parameter and µi is the Lagrange multiplier for the ith equality
constraint. In general for equality-inequality constrained problem, the multiplier method
can be written as [33]
λ g 1 2 λm jj j(x)+ rg (x) if g j + ≥ 0
Φ(x,g(x),r) = ΦE(x,h(x),r)+ ∑  2 j r (2.35)j=1 − 1 λ2 λ, if g + j2r j j r < 0,
where λ j ≥ 0 represents the Lagrange multiplier for the jth inequality constraint.
2.5 The sequential quadratic programming (SQP) method
The sequential quadratic programming (SQP) is arguably the most successful method
for solving nonlinear constrained optimization problems [8, 2, 33]. It is regarded as a con-
ceptual method from which different algorithms have emerged. The SQP was seemed to
have been invented in 1963 by R.B Wilson in his PhD thesis titled ”A simplicial algorithm
for concave programming” [41], where he called the method Newton-SQP algorithm. How-
ever, the method became popular in the late 70s through the works of Garcia Palomares and
33
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
Mangasarian [34] and Han [21], where methods for unconstrained optimization problems
were embedded in the SQP concept.
The SQP method basically models the nonlinear constrained problem for a given iterate
xk, k ∈ N0 by a quadratic programming (QP) subproblem, and then uses the solution to the
subproblem to construct a better iterate xk+1. Iterating this process leads to a sequence of
iterates that is hoped will converge to a local solution x⋆ as k → ∞. The method shares
some properties of Newton-type methods such as rapid convergence, especially when the
iterates are close to the solution. It may behave erratically when the iterates are far from the
solution, however some techniques have been invented to control this phenomenon [8].
The SQP method allows infeasible initial starting solution which is a big advantage
when solving a nonlinear programming problem with nonlinear constraints. The method
also depend largely on the existence of rapid and accurate algorithms and techniques for
solving quadratic programming problems. Fortunately, there are many good procedures for
solving quadratic programming problems.
2.5.1 SQP notations and model formulation
We briefly discuss here the required notations, process, and subquadratic model for the
SQP method. We also list some assumptions needed for the method, which include the
blanket assumption that all functions in the problem definition are three times continuously
differentiable. We denote by ∇ f (x) to be the gradient of a scalar function f , i.e., ∇ f (x) =
[ ∂ f , · · · , ∂ f ]T∂x ∂x . We also employ the same notation ∇ to denote the Jacobian of a vector-1 n
valued function, for example, ∇h(x) = (∇h1(x),∇h2(x), . . . ,∇hp(x)). We let H denote the
Hessian of a scalar-valued function, which is a symmetric matrix with entries of the form
2
H f (x) = ∂ f (x)i, j ∂xi∂x .j
Paramount to the SQP method is the Lagrangian function
L(x,µ,λ) = f (x)+µT h(x)+λT g(x), (2.36)
34
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
where µ ∈ Rp and λ ∈ Rm+ are the Lagrangian multiplier vectors for the equality and in-
equality constraints respectively. We define the index set of active constraints by
I(x) = {i : gi(x) = 0}, (2.37)
and we denote by G(x) ∈ Rn×(p+qx) the matrix given by
G(x) = (∇h1(x), . . . ,∇hp(x),∇gi1(x), . . . ,∇giq (x)), (2.38)x
where qx = |I(x)|. We denote by x⋆ ∈ Rn the local minimum point of the problem and
we write HL(x⋆,µ⋆,λ⋆), where µ⋆ ∈ Rp, λ⋆ ∈ Rm are the optimal multiplier vectors. We
assume the following conditions pertinent to the optimal solution (see [8, 2, 33] ).
(i) The first order necessary optimality condition: That is, there exist Lagrangian multi-
pliers µ⋆ ∈ Rp and λ⋆ ∈ Rm+ such that
∇L(x⋆,µ⋆,λ⋆) = ∇ f (x⋆)+∇h(x⋆)µ⋆+∇g(x⋆)λ⋆ = 0. (2.39)
This condition is also known as the KKT conditions. If a feasible point x satisfies the
KKT conditions, we call the point a critical point of the problem. Critical point may
not be a local minimum point.
(ii) The columns of G(x⋆) are linearly independent.
(iii) The strict complementary slackness holds at x⋆, that is, gi(x⋆)λ⋆i = 0, for j = 1, . . . ,m.
But if gi(x⋆) = 0, λ⋆i > 0.
(iv) The positive definite condition on HL⋆. That is, the Hessian of the Lagrangian with
respect to x is positive definite on the null space of G(x)T . That is dT HL⋆d > 0 ∀
d ̸= 0 such that G(x⋆)T d = 0.
35
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
The last three conditions are referred to as the second-order sufficient optimality conditions
which ensure that x⋆ is an isolated local solution of the problem and that the multipliers are
unique.
2.5.2 Quadratic subproblem formulation
The quadratic subproblem is constructed in a way that it reflects local attributes of the
original problem. First, we express f as its local quadratic approximation
1
f (x)≃ f (xk)+∇ f (xk)d + d Tx x BKdx, (2.40)2
where dx = x− xk and the direct or natural definition of matrix Bk is the Hessian of f at
iterate xk. The constraints are also expressed by their local linear approximations
g(x)≃ g(xk)+∇g(xk)dx
(2.41)
h(x)≃ h(xK)+∇h(xk)dx.
This leads to the following quadratic programming subproblem
1
min ∇ f (x )T Tk dx + dx B2 k
dx
dx
s.t ∇h(xk)dx +h(xk) = 0 (2.42)
∇g(xk)dx +g(xk)≤ 0.
The above formulation is reasonable and useful if the constraints are linear. The presence
of nonlinear constraints in many constrained optimization problems requires the use of a re-
lated approach that involves the Lagrangian as the objective function. The related approach
is given below:
1
min ∇L(x T Tk,µk,λk) dx + dx Bkd2 xdx
s.t ∇h(xk)dx +h(xk) = 0 (2.43)
∇g(xk)dx +g(xk)≤ 0.
36
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
The subproblems Eq.(2.42) and (2.43) are equivalent in some instances. First, if there is no
inequality constraints. This is because the term ∇h(xk)T dx is a constant and the functional
in (2.43) reduces to that in (2.42). Secondly, if the multiplier λk in (2.43) satisfies λik = 0,
for i∈ Iin(xk). This is because only the active inequality constraints are important and since
terms ∇h(xk)T dx and ∇g(xk)dx, i ∈ Iac(xk) are constants, we have equivalency.
Another idea of formulating the subproblem suggests the addition of slack variables
to the inequality constraints before constructing the subproblem. This idea leads to the
following
1
min ∇ f (xk)T dx + d Tx B2 k
dx
dx,ds
s.t ∇h(xk)dx +h(x (2.44)k) = 0
g(xk)+ sk +∇g(x )Tk dx +ds = 0,
where s = t2. The formulation is equivalent to the one below where the Lagrangian is used
as the objective functional:
1
min ∇L(x T Tk,µk,λk) dx + dx Bkdx
dx,ds 2
s.t ∇h(xk)dx +h(x ) = 0 (2.45)k
g(xk)+ sk +∇g(xk)T dx +ds = 0.
The quadratic subproblem is solved at each iteration and the solution together with the
multipliers is updated using:
xk+1 = xk +αdx
µk+1 = µk +αdµ (2.46)
λk+1 = λk +αdλ,
where α is the steplength and
dµ = µ−µk
(2.47)
dλ = λ−dλ.
37
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
2.5.3 The SQP algorithm framework
Consistency of the quadratic subproblem is very crucial to all SQP algorithm variants.
This implies that the associated system of constraints must have a non-empty feasible set.
This requirement is guaranteed when xk is close to x⋆ supported also by condition imposed
on matrix G(x). However, some difficult problems may fail consistency condition even
when xk is in the neighborhood of x⋆.
Another issue we need to consider is the convergence of the iterates. Iterates may
converge locally or globally. Local convergence is based on the assumptions that the initial
solution is relatively close to a local solution x⋆ and the initial Hessian is close to HL⋆.
Apart from the non-singular requirement of HL⋆, the following are the conditions imposed
on HL⋆ and Bk and are important for local convergence (see [33, 2, 8]):
(i) Uniform positive definiteness of Bk on the null space of ∇h(x )Tk is required. That is,
there exists β1 > 0 such that ∀k, dT Bkd ≥ β1∥d∥2 ∀d such that ∇h(xk)T d = 0.
(ii) The sequence of matrix {Bk} must be uniformly bounded. That is, there exists β2 > 0,
such that ∥βk∥ ≤ β2
(iii) The inverses of matrices Bk must be uniformly bounded. That is there exists β3 > 0
such that ∀k, B−1k exists and ∥B
−1
k ∥ ≤ β3.
For global convergence, we assume that that initial iterate may be remote and a func-
tion that measures the progress of convergence is often incorporated. This function is called
merit function denoted by Φ. It reduction at each step indicates a progress towards a solu-
tion. We present the SQP skeleton algorithm below [8, 2]:
38
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
Algorithm 1: SQP algorithm framework.
Result: Output the approximate solution vector for a given nonlinear constrained
problem
(1) Initialize (x0,µ0,λ0) and B0 ;
(2) Choose a merit function Φ and set k = 0;
(3) Form the QP and solve for (dx,dµ,dλ) ;
(4) Choose steplength α such that Φ(xk +αdx)< Φ(xk);
(5) Set
xk+1 = xk +αdx
µk+1 = µk +αdµ
λk+1 = λk +αdλ
;
(6) Stop if converged ;
(7) Otherwise compute Bk+1, set k = k+1 and return to 3.
Since the full constraint nonlinear problem can be converted to equality-constrained
problem by adding slack variables, we shall focus on the following equality-constrained
quadratic programming for easy explanation of different variants of SQP method
1
min ∇ f (xk)T dx + d T2 x
Bkdx
dx (2.48)
s.t ∇h(xk)dx +h(xk) = 0.
2.5.4 Newton version of the SQP
The Newton-type SQP provides an algortithm with a unit steplength [8, 2]. The lo-
cal convergence proof of SQP is also based on the Newton version. From the first order
39
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
necessary conditions apply to (2.48), we have
∇h(xk)µk+1 +Bkdx =−∇ f (xk)
(2.49)
∇h(xk)T dx =−h(xk),
where µk+1 is a vector of multipliers at k+1. Setting µk+1 = µk +dµ, we have
Bkdx +∇h(xk)dµ =−∇L(xk,µk)
(2.50)
∇h(xk)T dx =−h(xk).
Now, setting Bk = HL(xk,µk) and applying Newton’s method together with the first order
necessary conditions, we have the right-hand side vector as
 ∇L(x,µ)ψ x µ ( , ) =  = 0, (2.51)
h(x)
and the Jacobian matrix at the solution is
 
 HL⋆ ∇h(x⋆)J (x⋆,µ⋆ ) =  (2.52)
∇h(x⋆)T 0.
The matrix J (x⋆,µ⋆) is assumed to be non-singular. The Newton SQP scheme now requires
to solve
J (xk,µk)dx,µ =−ψ(xk,µk), (2.53)
where dx,µ = (dx,dµ) and the scheme is updated with a unit steplength. The initial iterate
for the multiplier can be computed using (see [8])
u0 =−[∇h(x0)T ∇h(x −10)] ∇h(x )T0 ∇ f (x0). (2.54)
40
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
Formula (2.54) follows the first order necessary condition and also referred to as the least
square estimate [2, 8].
The advantages of the Newton-SQP scheme include the benefit of a unit steplength and
rapid convergence with no need of line searches. It has been proven that if the initial iterate
(x0,µ0) is sufficiently close to (x⋆,µ⋆), the scheme converges quadratically. However, it
is practically difficult to choose an initial iterate that will be closer to the true local solu-
tion. Moreover, when the initial iterate is far from the true solution, HL(xk,µk) cannot be
assumed to be positive definite on the required subspace and this may not guarantee the
existence of local solution to (2.48).
2.5.5 Other variants of SQP methods
The updating scheme for the matrix Bk play an important role in SQP algorithms. As we
have noticed earlier, the Newton version uses a direct Hessian approximation as the matrix
Bk, which of course has some practical advantages. Other variants of SQP method emerge
from finding good approximation schemes for Bk. These schemes, however emerge from
the in-depth study of unconstrained optimization problems. The variants are divided into
two categories. The first called the full Hessian approximations approximate the matrix Bk
to be close to HL as possible while the second referred to as reduced Hessian approxima-
tions compute matrices that approximate the Hessian on the null space of the Jacobian of
the constraints [8]. One important property of the sequence {Bk} necessary for the approx-
imation of HL⋆ is the bounded deterioration property defined below:
Definition 2.18. [8, 9] A sequence of matrix approximation, Bk, for the SQP method is said
to have the bounded deterioration property if there exist constants α1 and α2 independent
of k such that
∥Bk+1 −HL⋆∥ ≤ (1+α α )∥B −HL⋆1 k k ∥+α2αk, (2.55)
where
αk = max{∥xk −xk∥,∥xk −x⋆∥,∥µ ⋆k+1 −µ ∥,∥µ −µ⋆k ∥}.
41
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
We discuss the full Hessian approach first, followed by the reduced Hessian method.
2.5.6 Full Hessian methods
Under the full Hessian approximation methods, we have (see [8])
(i) The rank-two Powell-Symmetric-Broyden (PSB) [21]. This method approximate
Bk using
1 − T − T − (y−BkS)
T S
Bk+1 = Bk + T (y BkS)S +S(y BkS) T 2 (SS
T ), (2.56)
(S S) (S S)
where
S = xk+1 −xk (2.57)
and
y = ∇L(xk+1,µk+1)−∇L(xk,µk+1). (2.58)
The sequence of matrices {Bk}, k > 0 may not be positive definite. It has been verified
that if x0 is close to x⋆, and B0 is close to HL⋆, then {Bk} satisfies the bounded
deterioration property and the iterates converges superlinear to a local solution. If
the initial solution is chosen remotely, the scheme may suffer the same drawback of
the Newton SQP, due to the positive indefiniteness of the sequence {Bk}. Moreover,
global convergence scheme maybe difficult to establish.
(ii) The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method. The BFGS scheme fo-
cuses on generating sequence of Bk approximations that are positive definite. The
formula is given by (see [8, 2, 33])
yyT B (SST− k )BB kk+1 = Bk + T T , (2.59)S y S BkS
and
ST y > 0. (2.60)
42
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
The formula possess the bounded deterioration property and if Bk is positive definite
and yT S > 0, then Bk+1 is positive definite. The main drawback is that condition
(2.60) may not be satisfied. The advantages of the scheme are attractive. First, initial
matrix B0 can be chosen to be an identity matrix. Secondly, if the problem is convex,
then it is easy to satisfy the positive definite requirement of HL⋆. Also, due to the
positive definite of Bk, the existence of solution to (2.48) is guaranteed even when
the initial iterates are chosen remotely. Lastly, BFGS scheme is also backed with
rigorous analysis in the realm of unconstrained optimization problems.
(iii) The Powell-SQP method. This scheme modifies the formula (2.59) by replacing y
by
ŷ = θy+(1−θ)BkS, for θ ∈ (0,1]. (2.61)
The modification still preserves the positive definiteness of the Hessian approxima-
tion and the condition (2.60) is always satisfied. However, the secant condition satis-
fied by (2.59) no longer apply. Although the local convergence proof of the method
is yet to be established, the method has performed experimentally very well.
(iv) The augmented Lagrangian method. This method transforms the objective func-
tion of the problem to
η
f (x) = f (x)+ ∥h(x)∥2A , η > 0, (2.62)2
whose Lagrangian function is given by
η
LA(x,µ) = L(x,µ)+ ∥h(x)∥2, (2.63)2
and Hessian at (x⋆,µ⋆) is
HL (x,µ) = HL⋆+η∇h(x⋆A )∇h(x⋆)T . (2.64)
43
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
BFGS method is then applied to the transformed objective function. The drawback
of the problem is the difficulty in choosing η. In some problems, large value of η
may be required to ensure the positive definiteness of HL⋆, which also requires the a
priori knowledge of x⋆. Numerical instability may ensue if an arbitrary large η value
is used.
(v) The SALSA SQP method. To circumvent the difficulties of the argumented La-
grangian SQP idea, the SALSA SQP method was invented and formulated in such a
way that independent of x⋆, a precise estimate of η can be chosen. Let
yA = ∇LA(xk+1,µk+1)−∇LA(xk,µk+1), (2.65)
then we can write
yA = y+η∇h(xk+1)h(xk+1), (2.66)
where y is given in (2.58). Positive definiteness of the augmented BFGS scheme is
preserved if the condition (2.60) is satisfied. The attractive feature also include the
idea that a minimum value of ηk independent on x⋆ can be given to ensure yTAS > 0
for values of x ⋆ ⋆k, µk+1 close to (x ,µ ). The drawback include the possibility that
HL⋆ could fail positive definite test despite the possibility of positive definite of the
sequence {Bk}. The bounded deterioration property may not be satisfied. Lastly, the
local convergence theory is not yet established.
2.5.7 The reduced Hessian method
The reduced Hessian approximation of sequence of matrices {Bk} is based on the con-
dition that HL⋆ is required to be positive definite on the null space of the gradient matrix
G(x)T defined in (2.38) [8, 11]. Thus the methods approximate only the portion of the
Hessian matrix pertinent to the subspace G(x)T . Attractive advantages of the method are
the size of the Jacobian matrix is reduced to n−m and standard positive definite updating
44
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
scheme can be employed. We shall discuss the method in [8, 11].
Let ∇h(xk) estimated at xk be a full rank matrix. Let us denote by Zk and Yk the ma-
trices whose columns form bases for the null space of ∇h(x)T and range space of ∇h(xk)
respectively. Matrices Zk and Yk can be computed by a QR factorization of ∇h(xk). We
also assume that the columns of Zk are orthogonal. The definition of a reduced Hessian is
given below:
Definition 2.19. Let (xk,µk) be a given iterate and ∇h(xk) has a full rank. The ma-
trix ZTk HL(xk,µk)Zk is called a reduced Hessian for the Lagrangian function estimated
at (xk,µk).
The term of the reduced Hessian depends on the choice of the basis for the null space
of the matrix ∇h(xk)T and hence not unique. The reduced Hessian at (x⋆,µ⋆) is positive
definite and if (xk,µk) is chosen close to (x⋆,µ⋆) the reduced Hessian at (xk,µk) will also be
positive definite.
Let us decompose dx as follow:
dx = Zk pz +Yk py, (2.67)
then the constraint in (2.48) becomes
∇h(xk)TYk py =−h(xk), (2.68)
and we can solve for py using
py =−[∇h(xk)TYk]−1h(xk). (2.69)
The QP subproblem in (2.48) reduces to
1
min pTz Z
T T T
p 2 k
BkZk pz +(∇ f (xk + py Bk))Zk pz (2.70)
z
45
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
and pz satisfies
pz =−ZTk BkZk[(∇ f (xk)+ pTy Bk)Zk]. (2.71)
The reduced matrix can be updated directly instead of updating Bk and recomputing ZTk BkZk.
Let Rk be a given (n−m)× (n−m) positive definite matrix at xk and µk, we can compute
the next solution by first computing py using (2.69) and setting
p =−R−1z k Z
T
k (∇ f (xk)+ p
T
y Bk). (2.72)
With pz and py computed, the next iterate is computed by setting
xk+1 = xk +dx, (2.73)
where dx is given by (2.67). To continue the iteration, a new multiplier µk+1 is needed and
this can be computed using the least square formula (2.54). Also to update Rk, the BFGS
scheme (2.59) can be employed by replacing Bk with Rk, and setting
S = ZTk (xk+1 −xk) = ZTk Zk pz
(2.74)
y = ZTk [∇L(xk +Zk pz,µk)−∇L(xk,µk)].
2.5.8 Merit function for SQP
Merit functions denoted Φ are often incorporated in SQP algorithms to monitor global
convergence (see [8]). They are usually constructed in a way that the solutions of the
nonlinear constrained problem are the unconstrained minimizers of the merit functions and
satisfy
Φ(xk +αdx)< Φ(xk), (2.75)
where α is an appropriate steplength. Merit functions in constrained problems have pa-
rameters that control the step reduction of f and the feasibility of the solution xk at each
iteration.
46
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
To ensure a steady decrease of a merit function, the following conditions are imposed []
(i) The initial vector x0 and all subsequent iterates xk, k ∈N are located in some compact
set C ∈ Rn
(ii) The columns of ∇h(x) are linearly independent for all x ∈ C .
We present the two commonly used merit functions below.
(i) The augmented Lagrangian merit function
The augmented merit function of the form [40, 8]
η
ΦF(x;η) = f (x)+h(x)T µ(x)+ ∥h(x)∥2, η > 0 (2.76)2
have been employed in literature. The multiplier µ(x) can be estimated using the least
square equation (2.54), and therefore µ(x⋆) = µ⋆. Based on the above conditions, ΦF
and µ are differentiable and ΦF is bounded below on C for sufficiently large η. It has
been proven that x⋆ ∈ C is a strong local minimum of ΦF if and only if x⋆ is a strong
local minimum of the original nonlinear problem, provided that the above conditions
are satisfied. However, if x is not a critical point of the original problem, then dx is a
descent direction for ΦF .
The augmented merit function can be extended directly to handle inequality con-
straints in two ways. First by the active set strategy where some set of inequality
constraints are treated as equality constraints. Let
−λ
I ik = {i : gi(xk)≥ }, (2.77)η
where at xk, Ik contains all unsatisfied constraints but no satisfied ones. The merit
47
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
function can then be defined as
1 1 1
ΦFI(x,µ,λ;η)= f (x)+h(x)T µ(x)+ η∥h(x)∥2+ ∑ (gi(x)λi(x)+ ηgi(x)2)+ (λ(x))2.2 i∈I 2 2η ∑ ik i∈/Ik
(2.78)
The second extension uses the idea of slack variable formulation of the quadratic
subproblem, where squared slack variables are added to g(x)
min f (xk)
x,t
s.t h(x) = 0 (2.79)
g(x)+ t2 = 0.
If we let z = t2j j , for j = 1, . . . ,m, we have an adjusted equality constraint vector
  h(x)h(x,z) =  (2.80)
g(x)+ z.
The merit function is then constructed as
T ηΦ f z(x,z) = f (x)+h(x,z) µ(x,z)+ ∥h(x,z)∥2, (2.81)2
where µ(x,z) can be computed using the least square estimate (2.54) and the slack
variables step dz can be computed after each iteration using
dz =−(∇g(x Tk) dx +g(xk)+ zk). (2.82)
Note that equation (2.82) is derived from the linearization of g(x)+ z = 0.
(ii) The ℓ1 merit function [22, 8] is another applicable merit function for the SQP al-
gorithm. When the problem is defined with only equality constraints, we have it to
48
2.5. THE SEQUENTIAL QUADRATIC PROGRAMMING (SQP) METHOD
be
Φ1(x;ρ) = f (x)+ρ∥h(x)∥1, ρ > 0 (2.83)
while for the full constrained nonlinear problem we have it to be
Φ1(x) = f (x)+ρ[∥h(x)∥ +∥g+1 (x)∥1], (2.84)
where 
g+(x) =0 if gi(x)≤ 0 (2.85)gi(x) if gi(x)> 0
The ℓ1 merit function is a form of penalty function just like the augmented La-
grangian. That is, there exists ρ⋆ such that for all ρ ≥ ρ⋆, if x⋆ is an unconstrained
minimum of Φ1 then x⋆ is a solution of the nonlinear problem. Note that Eq. (2.83) is
not differentiable on the feasible set and the term ρ∥h(x)∥1 cannot be squared because
the exact penalty property will be lost. The merit function is constructed in such a
way that as ρ → ∞, the minimizer of Φ1 converges to the solution of the nonlinear
constrained problem.
49
Chapter 3
Minmax capacity provisioning on path
networks
In this chapter, we develop and analyse a novel algorithm and data-structure suitable for
evacuation procedure in remote and sparsely populated communities. When considering
remotely or sparsely populated areas, such as the sparse areas of British Columbia, Cali-
fornia and New South Wales, evacuation process may not be feasible by foot as assumed
in many (traditional) sink location algorithms (see Chapter 1). As a matter of urgency,
emergency vehicles or aircraft may have to be dispatched (or already situated) in order
to evacuate the distressed population. In this case also, congestion may not be the primary
concern of the evacuation planners and hence, there is need for a new evacuation procedure.
To explore this problem, we introduce the concept of capacity provisioning for evac-
uation on dynamic networks as an extension of the traditional sinks location problems,
suitable for the effective evacuation planning of the aforementioned communities. Unlike
the usual sink location problem (irrespective of the criterion) which always incorporate road
capacities in the model development, this new concept involves the determination of the op-
timal capacities distribution from a fixed budgeted capacity under the assumption that the
location of the sink is known. In other words, the budgeted capacity is viewed as a resource
to be allocated optimally to the edges of the network such that the evacuation time of the
evacuees to the sinks is minimized. We formerly present the problem definition and model
in the next section, specifically for a path topology.
50
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
3.1 Model formulation and problem definition
Using the parameters defined in Chapter 1, We define our problem as follows:
Problem 3.1. Given a dynamic path network N = (P,c,τ,w), a fixed sink s ∈ P, and a total
capacity budget c. Let z(x) denote the evacuation time which represents the time for the
last evacuee to get to sink s given the capacity vector x = (x1,x2, . . . ,xn) where xi represents
the capacity assigned to edge ei on the path. Find the optimal capacity vector x⋆ such that
z(x⋆) is minimum and ∑ni=1 xi = c.
Just like the traditional sink location problem on a path, we observe that two inde-
pendent flows emerge on either side of the sink (see Chapter 1). If we denote by zL(x)
(resp. zR(x)) the evacuation time of evacuees on the left (resp. right) side of the sink, then
z(x) = max{zL(x),zRx}. We shall focus on function zL(x) and treat zR(x) symmetrically.
For simplicity, we assume a path with n vertices with the sink s connected to vertex vn via
edge en. Thus, all evacuee vertices are on the left side of the sink. From Higashikawa et
al. [24] and Bhattacharya et al. [5], we learn that the objective function or cost can be
expressed in terms of n evacuation time functions θi(x)
zL(x) = max {θi(xi)}, (3.1)
i∈[1,n]
where
∑i
θ j=1
w j
i(xi) = + τd(vi,s) (3.2)x(vi,s)
and x(vi,s) denotes the minimum edge capacity between a vertex vi and the sink s, that is
x(vi,s) =mini≤ j≤n x j. The first term in Eq. (3.2) computes the waiting time (at vi) of the last
evacuee from v1 if the congestion occurs at vi. We next illustrate several simple instances
of the capacity provisioning problem.
Illustration 3.1.0.1 (Sink in the middle with two evacuee nodes). Let us consider the simple
path network shown in Fig. 3.1, where the sink is located at the middle of the path. We
denote the budgeted capacity by c and the number of evacuees at v1 and v2 by w1 and w2
51
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
respectively. Our objective is to find optimal values of edge capacities x1 and x2 such that
the completion time of all evacuees is minimized. Note also that x1 + x2 = c and both edge
capacities are non-negative. For simplicity, let x2 = x so that x1 = c− x.
w1 c− x x w2
v1 l1 s l2 v2
Figure 3.1: Simple path network with sink in the middle.
We understand that the evacuation times (or costs of evacuation) for evacuees on both
sides of the sink are given by the following evacuation functions (see [24, 5]).
w
θ x 11( ) = + τlc− x 1 (3.3)
w
θ 22(x) = + τlx 2
.
Intuitively, the first term of the functions represents the time the last evacuee needs to wait
at the node before it can proceed on the edge. This time is determined by the edge capacity.
The second term represents the travel time through the edge. We notice that the optimal
solution x∗ has the following property: θ (x∗1 ) = θ2(x∗). If this is not the case, assume
θ1(x∗) > θ2(x∗), thus z(x∗) = θ1(x∗). Then there exists ε > 0 such that the solution x∗−
ε, which is obtained by transferring ε capacity from edge (s,v2) to edge (s,v1), is still
determined by θ (x∗1 − ε). But then θ2(x∗− ε) ≤ θ1(x∗− ε) = z(x∗− ε) < θ1(x∗) = z(x∗)
which contradicts the assumption that x∗ is optimal. Based on this property, we have the
following quadratic equation
x2(l2 − l1)τ+ x(w1 +w2 + τl1c− τl2c)−w2c = 0. (3.4)
If l1 = l2, then we obtain x2 = x = w2c w1cw1+w and x = c− x = .2 1 w1+w2
Without loss of generality, suppose that l1 > l2. Then, the optimal solution is given by
one of the roots of Eq. (3.4).
52
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
We can show that only one of the roots of the equation is feasible for our problem. The
roots of Eq. (3.4) are
√( )
w 2
x 1
+w2 c 1 w1 +w2 4w2c
a,b = + ± + c −2τ(l1 − l2) 2 2√ τ(l1 − l2) τ(l1 − l2)
w1 +w2 c 1 (w1 +w2)2
= + ± c2 2(w1 −w2)c+
2τ(l1 − l2) 2 2 τ2(l1 − l2)2
+ .
τ(l1 − l2)
We analyse the first root x = xa as follow :
Case 1: When w1 > w2, we have √
w1 +w c 1 (w +w )2 2(w −w )cx 2 c2 1 2 1 2a = + + + +2τ(l1 − l2) 2 2√( τ2(l1 − l2)2) τ(l1 − l2)
w1 +w2 c 1 (w1 −w2) 2x (3.5)a > + + c+2τ(l1 − l2) 2 2 τ(l1 − l2)
w
= c 1+ .
τ(l1 − l2)
Case 2: When w1 < w2, we have √
w1 +w2 c 1 (w1 +w2)2xa = + + √c2 + 2 2 −
2(w2 −w1)c
2τ(l1 − l2) 2 2 ( τ (l1 − l2) ) τ(l1 − l2)
w1 +w2 c 1 − (w2 −w )
2
x c 1 (3.6)a > + +2τ(l1 − l2) 2 2 τ(l1 − l2)
w
= c 1+ .
τ(l1 − l2)
The case when w1 = w2 is obvious. Thus xa is more than the total allocated budget c and is
infeasible. The same idea can be employed to show that 0 < xb < c. We present the proof
53
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
for only the case w1 > w2, as the case w2 > w1 follows the same idea. We have
√
w1 +w2 c 2xb = + −
1 2 (w1 +wc 2
) 2(w1 −w2)c
+ +
2τ(l1 − l2) 2 2√ τ2(l1 − l2)2 τ(l1 − l2)
w1 +w2 c − 1 2 (wc 1 +w2)
2 2(w1 +w2)c
> + + + (3.7)
2τ(l1 − l2) 2 2√( τ2(l1 − l2)2) τ(l1 − l2)
w 21 +w2 c − 1 (w1 +w2)= + c+ = 0.
2τ(l1 − l2) 2 2 τ(l1 − l2)
Next we show that xb < c √
w
x 1
+w2 c 1 2 (w1 +w2)
2 2(w1 −w2)c
b = + − c + +2τ(l1 − l2) 2 2√ τ2(l1 − l2)2 τ(l1 − l2)
w1 +w2 c − 1 2 (w1 +w2)
2
− 2(w1 +w2)c< + c + (3.8)
2τ(l1 − l2) 2 2√( τ2(l1 − l2)2) τ(l1 − l2)
w1 +w2 c − 1 w1 +w
2
2
= + − c = c.
2τ(l1 − l2) 2 2 τ(l1 − l2)
Algorithm 2: Capacity allocation algorithm for dynamic path network with two
evacuee nodes and a sink at the middle.
Result: Output optimal cost and vector of edge capacities for the network in
Fig. 3.1
Inputs: w1, w2, c, τ, l1, l2;
if l1 = l2 then
Assign x⋆ w2c2 = w +w ;1 2
else
Assign x⋆2 = xb, where xb is the root of Eq. (3.4) that lies in (0,c) ;
end
Compute z(x⋆ ⋆2) = θ2(x2) ;
Output the optimal cost z(x⋆2) and vector (c− x⋆2,x⋆2)
Illustration 3.1.0.2 (Path with two edges and sink at an endpoint). Let us now consider
54
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
another simple path network (shown in Fig. 3.2)
w1 l1 w2 l2
v1 c− x v2 x s
Figure 3.2: Illustration of capacities allocation on a two-edged path network with the sink
located at the end vertex on the right.
where the sink is denoted by s on the right. The case where the sink is located on the left
can be treated symmetrically. We again denote x2 by x and x1 = c− x. From the figure, the
time taken for the first evacuee at v2 to reach the sink is τl2. Since the last evacuee at v2
must wait up to w2x , then the evacuation time for the last evacuees at v2 is
w2
x + τl2. At the
same time, evacuees from v1 start moving towards the sink. We have to consider two cases:
Case 1: (Non-congestion scenario) Evacuees from v1 arrive at v2 and v2 is empty (its
evacuees have already left v2). For this case, we can write the completion time as
w
z x 1( ) = θ1(x) = + τ(l + lmin{c− x x} 1 2
). (3.9)
,
Case 2: (Congestion scenario) Evacuees from v1 get to v2 but are being delayed due
to backlog of evacuees waiting to evacuate through edge e2.The completion time for this
scenario is (see [5, 24])
w
z x θ x 1
+w2
( ) = 2( ) = + τl2. (3.10)x
According to authors of [5] and [24], the two cases can be merged into a single relation:
z(x) = max{θ1(x),θ2(x)}. (3.11)
We observe the following property (see Theorem 3.2) of the optimal solution vector
called the monotonicity property. The property states that the allocation of capacities must
be non-decreasing towards the sink. This helps us further simplify the expression for the
evacuation time. In addition, we now know that x⋆ ≥ c2 . We shall discuss the proof of the
55
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
property in Section 3.2.
Theorem 3.2. Given a path P with n+ 1 vertices where v1, . . .vn are occupied by evac-
uees w1, . . . ,wn and the sink is located at v ⋆n+1. Let x = (x⋆, . . . ,x⋆1 n) be the optimal edge-
capacities vector, then x⋆ ⋆i ≥ x j if i > j.
Based on Theorem 3.2, we rewrite both time functions as follows:
w
θ x 11( ) = + τ(l1 + l2)c− x (3.12)
w1 +wθ 22(x) = + τl2.x
To obtain the optimal solution for z(x) in Eq. (3.11) with θ1 and θ2 defined in Eq. (3.12),
we first examine the nature of both evacuation time functions θ1 and θ2. We observe
that evacuation function θ1(x) (resp. θ2(x)) is increasing (resp. decreasing) with x and
limx→c θ1(x) = ∞. Two situations are possible. The two evacuation functions may either
intersect (θ1( c2)≤ θ (
c
2 2)) or not (θ1(
c
2)> θ2(
c
2)) for x∈ [
c
2 ,c]. We explain the consequences
of the two situations below:
Case 1 : If θ1 and θ2 do not intersect for x ∈ [ c2 ,c], then, z(x) = θ1(x) and x =
c
2 is the
optimal solution with the optimal cost θ1( c2).
Case 2 : If θ1 and θ2 intersect, then the intersection point determines the optimal solution.
These two cases are determined by the values of θ1 and θ2 at point c2 . To identify these
cases, we solve the following equation.
(c) (− c) 2wθ1 θ2 = τl − 21 = 0, (3.13)2 2 c
and consequently we have
2w
c 2= . (3.14)
τl1
We call the value of the capacity in Eq. (3.14) critical capacity and we denote it by ccritical .
The above shows that if we have any budgeted capacity c greater than ccritical , allocating
56
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
equal capacities is enough to attain optimality. However, if the two functions intersect at a
point in (c/2,c), then this point is the value of x that minimizes the cost function z(x). To
compute this point, we require θ1(x)− θ2(x) = 0, which leads to the following quadratic
equation
τl x21 − (2w1 +w2 + τl1c)x+(w1 +w2)c = 0. (3.15)
The solutions to Eq. (3.15) are
√
2w
x 1
+w2 + τl1c ± (w2 − τl1c)
2 +4w1(w1 +w2)
a,b = 2τl1 4τ2l2
. (3.16)
1
We now show in the lemma below that only one of the solutions in Eq. (3.16) lies in the
open interval (c/2,c) and the other can be discarded.
Lemma 3.3. Given a path network with 3 vertices {v1,v2,v3} where the sink is located on
v3. If the total capacity c is less than the critical value ccritical , then the optimal solution is
one of the zeros of Eq. (3.15) presented in Eq. (3.16).
Proof. Based on the condition τl < 2w21 c from Eq. (3.14), we first show that the solution
given by √
2w1 +w2 + τl1c (w2 − τl1c)2 +4w1(wx 1 +w2)a = + . (3.17)2τl1 4τ2l21
is more than the budgeted capacity c to be shared. We analyse xa as follow:
√ [ ] ( )
w1 w2 c 1 w
2 2
x 2a = + + √+ −
c w1 w1w2
+ +
τl1 2τl1 2 [4 τl1] 1( ) τl1 (τl1)2
≥ w1 c c 1 c
2
− w
2
1 cw1
+ + √+ [c + ] +τl1 4 2 4 2 τl1 2τl1
3c w1 c2 w1 w1 c (3.18)
= + + + +
4 τl1 √(16 τl1 )τl1 2
3c w1 c w
2
1
= + + +
4 τl1 4 τl1
2w1
= c+ ,
τl1
57
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
Thus xa is more than our budgeted capacity and not feasible.
With c < 2w2τl , we analyse xb as follow:1 √ [ ] ( )
w 21 w2 c 1 w2 c w
2
1 wx 1
w2
b = + + − − + +τl1 2τl1 2 √4 [τl1 1 ] τ(l1 ) (τl1)2
w1 w2 c − 1 w
2 2
2 2w2 w1 w1w2
> + + − + + (3.19)
τl1 2τl1 2 √(4 τl1 τl1) τl1 (τl1)2
w w 21 2 c w2 w1 c
= + + − + = .
τl1 2τl1 2 2τl1 τl1 2
Next,(with c < 2w2τl ) we show that xb < c as follow:1 √ [ ] ( )
w1 w2 c − 1 w2 − c
2 w 21 w1wx 2b = + + √ + +τl1 2τl1 2 [4 τl1] 1 ( )τl1 (τl1)2w 21 c c 2− 1 c − c w1 w1c< + +
τl1 4 2√ + +4 (2 )1 τl1 2τl1
w1 3c − c
2 w 21 w1c
= + √ + + (3.20)τl1 4 16 (τl1) 2τl1
w 3c c2 w 21 − 1 − w1c< + +
τl1 4 √(16 τl1) 2τl1
w1 3c w1 c
2
= + − − = c.
τl1 4 τl1 4
Since x cb ∈ (2 ,c), x = xb is the optimal solution and the optimal solution vector is x
⋆ =
(c− xb,xb).
Definition 3.4. Whenever the optimal capacities for two consecutive edges are not obtain-
able by allocating equal capacity, we call this situation a touching scenario, and the solution
method a touching solution. The solution xb in Eq. (3.19) is an example of a touching solu-
tion.
We can summarize our findings in the following algorithm for the case of a path network
with two evacuee nodes and the sink node located at one of the endpoints of the path.
58
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
Algorithm 3: Capacity allocation algorithm for dynamic path network with two
evacuee nodes and one sink at an endpoint.
Result: Output optimal cost and vector of edge capacities for a two-edged network
Inputs: w1, w2, c, τ, l1, l2;
Compute the critical capacity c 2w2critical = τl ;1
if If c ≥ ccritical then
Assign equal capacity x⋆ = c− x⋆ = c2 ;
else
Let x⋆ = xb, where xb is computed in (3.19);
end
Compute z(x⋆) = θ ⋆1(x ) ;
Output the cost z and vector (c− x⋆,x⋆)
Illustration 3.1.0.3 (Path network with four evacuee nodes, two nodes on each side of the
sink). Let us consider a more interesting case where we have a path network with four
edges with the sink located in the middle of the path (see Fig. 3.3).
w1 l1 w2 l2 l3 w4 l4 w5
v1 x1 v2 x2 s x3 v4 x4 v5
Figure 3.3: Four-edged path network with sink at the middle.
Based on Theorem 3.2, we require a non-decreasing allocation of capacities towards the
sink from both sides. In other words, we require x2 ≥ x1 and x3 ≥ x4. From the path graph,
we have the following evacuation functions
z (xTL ) = z1,2(xT ) = max{θ1(x1),θ2(x2))}
(3.21)
zR(xT) = z T3,4(x ) = max{θ3(x3),θ4(x4))},
59
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
where
w
θ1(x
1
1) = + τ(l1 + l2)x1
w1 +wθ 22(x2) = + τlx 22 (3.22)
w +w
θ 3 43(x3) = + τlx 33
w
θ4(x
4
4) = + τ(l + l ).x 3 44
Before we present an algorithm for this problem, we prove the following lemma which
gives us an insight into how the optimal solutions can be obtained.
Lemma 3.5. At the optimal solution x⋆,
z ⋆1,2(x ) = z ⋆3,4(x ). (3.23)
Proof. Suppose the above is not true and without loss of generality, let us assume z (x⋆1,2 )>
z ⋆3,4(x ). Then there exists ε > 0 so that we decrease x3 and x4 by ε each and increase x1
and x2 by ε each to obtain a new vector x satisfying the conditions z1,2(x) ≥ z3,4(x) and
z ∗ ⋆1,2(x)< z1,2(x ). This contradicts the fact that x is the optimal solution.
Based on the lemma above, we sketch an algorithm to solve this more general problem.
We observe that, each sub-path on either side of the sink reduces to a problem which we
know how to solve if we are given the total budget capacity for the sub-path. Let X1 be
the total capacity allocated to one of the sub-paths, which we call the left sub-path. We
extend our notation to denote the optimal solution cost for the left sub-path if the total
capacity allocated to the left sub-path is X1 by zL(X1). Similarly, we have zR(c−X1) as the
optimal cost for the right sub-path. In addition, the actual solutions x1, . . . ,x4 are parametric
solutions with X1 as the parameter. Algorithm 4 illustrates the steps we can take to solve
this problem instance.
60
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
Algorithm 4: Capacity allocation algorithm for the instance in Fig. 3.3
Result: Output optimal cost zL(X∗1 ) and vector x
∗
Inputs: w1, w2, w3, w4 c, τ, l1, l2, l3, l4;
Compute the critical capacities c = 2w21,2 τl and c
2w4
1 3,4
= τl ;4
Solve the left sub-path problem with budget X1 = c1,2 and the right with budget
X2 = c−X1;
if If zL(X1)≥ zR(X2) then
X⋆ ≥ c and x⋆ = x⋆1 1,2 1 2 =
X1
2 is the parametric solution;
zL(X1)← 2w1X + τ(l1 + l2) ;1
else
X⋆1 < c1,2 and zL(X1)← f2(xb), where x∗2 = xb (given by (3.19)) is a parametric
solution that depends on X1 and x∗1 = X1 − x∗2;
end
Solve the left sub-path problem with cost X1 = c− c3,4 and the right with cost
X2 = c3,4 ;
Repeat the test for the right sub-path;
Obtain function zR(X1) and the corresponding parametric solutions x∗3 and x
∗
4 as
functions of X1;
Solve equation zL(X1) = zR(X1) numerically, and obtain X∗1 ;
Output the optimal cost zL(X⋆1 ) and the capacity vector (x
⋆,x⋆1 2,x
⋆ ⋆
3,x4) by
evaluating the parametric solutions at X∗1 ;
Before we present the general properties of the optimal solution for a given arbitrary
path network, let us briefly discuss the running time of Algorithm 4. We note that, in the
first part of the algorithm, we obtain a parametric optimal solution x∗ which depends on
the unknown capacity budget X1 assigned to the sub-path left of the sink. The time com-
plexity for this computation is O(1). The more challenging phase of the algorithm is to
determine the value for the budget X1 so that zL(X1) = zR(X1). The expressions zL(X1) and
61
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
Θ
zL zR
solution returned
upper bound
X− X∗ X+
0 c
δ
Figure 3.4: Binary search procedure for the optimal capacity budget X1.
zR(X1) are linear functions of the reciprocal of the parametric solution. In the worst case,
the parametric solutions depend on budget X1 as illustrated by Eq. 3.16. It is impractical to
attempt to find an algebraic expression for X1 that solves equation zL(X1) = zR(X1). This
equation can be transformed into a complicated polynomial equation of degree larger than
four in X1 and we know that there is no algebraic expression for the roots of a general
polynomial equation of degree higher than four [1]1. We opt to find the value for X1 ap-
proximately using a simple binary search procedure. We find next an upper bound on the
number of steps needed by this binary search procedure to obtain a solution with evacuation
time Z ≤ (1+ ε)Z∗ for any ε > 0, where Z∗ is the evacuation time of the optimal solution.
Let δ be the difference between the upper and lower bounds on the optimal X∗1 that is main-
tained by the binary search procedure. Since the initial binary search range for X1 is interval
(0,c) where c is the value of the total capacity budget, the number of binary search itera-
tions performed is O(log cδ). Let X
−
1 and X
+
1 be the lower and upper bounds respectively
1This doesn’t mean that an algebraic solution for this particular equation does not exist.
62
3.1. MODEL FORMULATION AND PROBLEM DEFINITION
on the optimal X∗1 maintained by the binary search (see Fig. 3.4). We know that zL(X1) is
a decreasing function of X1 and zR(X1) is increasing. During any iteration of the binary
search procedure, value min{z (X−L 1 ),zR(X
+
1 )} can be returned as an approximate solution
and max{zR(X−1 ),zL(X
+
1 )} is a lower bound on the optimal solution. The binary search can
stop when min{zL(X−),z (X+1 R 1 )} ≤ (1+ ε)max{z
−
R(X1 ),zL(X
+
1 )}. It can be verified that
this condition is satisfied if,
zR(X−1 +δ)≤ (1+ ε) · z (X
−
R 1 ) or (3.24)
z − −L(X1 )≤ (1+ ε) · zL(X1 +δ). (3.25)
We consider first function zL(X1) and condition (3.25) to identify a suitable constraint on
δ. To simplify the calculations, we note from Algorithm 3 that there are two cases for
calculating the capacity assignment of the left sub-path. In both cases, the capacity assigned
to an edge such as (v2,s) in Fig. 3.2 for the left sub-path grows linearly with the budget
X1 (we note that the unknown parameter X1 corresponds to budget c in Eq. (3.16)). The
evacuation time zL(X1) grows asymptotically as a fractional expression in X1. To simplify
the calculations, we take z aL(X1) = gX1+h for some constants a,g,h with a,g > 0. If we
substitute z −L(X1 ) in (3.25), we obtain
a
− ≤
a
(1+ ε) − , . . . (3.26)gX1 +h g(X1 +δ)+h
gX−1 +h≤ hδ ε = εX−1 + ε . (3.27)g g
′
Similarly, we can write z aR(X1) = h′−g′X which we substitute in (3.24).1
a′ ′≤ a(1+ ε) , . . .
h′−g′X+ h′−g′(X+1 1 −δ
h′−g′X+1 h′−g′(X
−−δ) h′−g′X−
δ ≤ ε ′ = ε
1
′ = ε
1 + εδ.
g g g′
63
3.2. PROPERTIES OF THE OPTIMAL SOLUTION FOR AN ARBITRARY PATH
WITH THE SINK LOCATED AT AN ENDPOINT
We consider the tighter inequality that implies the inequality above,
h′−g′X− ′
δ ≤ ε 1 h′ = ε ′ − εX
−
1 . (3.28)g g
We notice that the upper bounds on δ from Eq.(3.27) and (3.28) can be arbitrarily small
′
depending on the value of X−1 . However, the sum of the upper bounds, ε(
h h
g + g′ ), does not
′
depend on X−1 . This means that, when δ ≤
ε h h
2(g′ + g), at least one of eq. (3.24) and (3.25) is
satisfied and we obtain a (1+ ε) approximate solution. We conclude that, the running time
of Algorithm 4 is O(log cε).
3.2 Properties of the optimal solution for an arbitrary path with the
sink located at an endpoint
In this section, we study and analyse the properties of the optimal solution vector x⋆
exploitable to develop an efficient algorithm given an arbitrary path with the sink located
at the endpoint. Note that x⋆ is a vector of capacities which, when allocated to the edges
of the path network, minimizes the completion time and also satisfies the capacity budget
condition x⋆ · 1 = c, where 1 is the vector with all elements equal to 1. We recall the
definition of evacuation functions θi(xi) for the case when the sink is located on point vn+1,
the right endpoint of the network
W [1, i]
θi(xi) = + τL[i], (3.29)x(vi,s)
where W [1, i] = ∑ij=1 w j and L[i] = ∑
n
j=i l j. We also let L[1, i] = ∑
i
j=1 l j.
The first property of the optimal solution is the monotonicity property given in Theo-
rem 3.2. We have exploited and employed this property in Illustration 3.1.0.2 and 3.1.0.3.
We present the proof below.
Proof of Theorem 3.2. Let us assume that the optimal solution, x⋆ is not monotonic in-
64
3.2. PROPERTIES OF THE OPTIMAL SOLUTION FOR AN ARBITRARY PATH
WITH THE SINK LOCATED AT AN ENDPOINT
creasing towards the sink at vn+1. Then there exist xi > xi+1 in the ordering of x⋆. Let x′ be
another vector with the same elements as x⋆, except that x′ ′i+1 = xi and xi = xi+1. We notice
that θ ′j(x ) = θ j(x⋆) for all j ≤ i and j ≥ i+2. We observe that
W [1, i+1]
θi+1(x⋆) = + τL[i+1], (3.30)min{xi+1,x(vi+2,s)}
and
′ W [1, i+1]θi+1(x ) = + τL[i+1]min{x′i+1,x(vi+2,s)} (3.31)
W [1, i+1]
= + τL[i+1].
min{xi,x(vi+2,s)}
Since xi > xi+1, then min{xi,x(vi+2,s)} ≥ min{xi+1,x(vi+2,s)} and θ ′i+1(x ) ≤ θi+1(x⋆).
Therefore, it is either the cost is the same and x′ is OPT or cost decreases which is a
contradiction that x⋆ is OPT.
Based on Theorem 3.2, we have the following definition which redefines the evacuation
time for any index i in terms of the discovered minimum capacity and therefore definition
in Eq.(3.2) changes.
Definition 3.6. Given an instance of a capacity provisioning problem on a path network
with a sink located at the right end, the associated evacuation time for an index i is given by
W [1, i]
θi(xi) = + τL[i]. (3.32)xi
The second theorem gives us an adjustment to the monotonicity property and shows
the possibility of assigning equal capacity to adjacent edges. The theorem also gives us a
candidate of the optimal cost amongst many of the time functions.
Theorem 3.7. Let P be a path network with n+ 1 vertices, with the sink located at the
(n+1)th vertex. Let h⋆ be the largest index defining the optimal cost z⋆, i.e.,
⋆ W [1,h
⋆]
z = θh⋆ = + τL[h⋆].xh⋆
65
3.2. PROPERTIES OF THE OPTIMAL SOLUTION FOR AN ARBITRARY PATH
WITH THE SINK LOCATED AT AN ENDPOINT
Then the optimal capacity vector has the following properties
(i) x⋆ = x⋆ ⋆h⋆ h⋆+1 = . . .= xn
(ii) Let i1, i2, i3 . . . , ik be the indices that define z⋆, enumerated in increasing order, i.e. i1 <
i2 < i3 . . . < ik and θi1 = θi2 = θi3 . . . = θik = θh⋆ . Then, i1 = 1. Moreover, if i j+1 >
i j +1 for some 1 ≤ j < k, then xi j = xi j+1 = . . .= xi j+1−1 for all j ∈ {1, . . .k−1}.
This theorem states two important properties of the optimal edge capacity vector. First,
the evacuation time is always determined by the evacuation function θ1. Second, if other
evacuation functions are equal in value with θ1, namely θi2, . . . ,θik , then these indices de-
termine the only increasing sequence of elements in the optimal capacity vector.
Proof. Suppose there exists an index j ≥ h⋆ such that x j+1 > x j. By definition, θ j > θ j+1.
Therefore, we can write
W [1, j] W [1, j+1]
+ τL[ j]> + τL[ j+1]
x j x j+1 (3.33)
W [1, j] W [1, j+1]
+ τl > .
x jj x j+1
This means that we can decrease the capacity x j+1 by an appropriate value ε > 0 in such
a way that inequality θ j+1 < θ j is maintained. We then distribute the extra capacity ε to
all the other edges on the path, namely to x1, . . . ,x j and to x j+2, . . .xn, in such a way as to
maintain the monotonicity property of the vector of capacities. As a consequence, all of the
evacuation functions θ1, . . .θ j and θ ⋆j+2 . . .θn decrease. Therefore z decreases which is a
contradiction.
To prove the second part of the theorem, we use the diagram in Fig. 3.5.
66
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
Distribute the  capacity
Distribute the  ca-
x pacitya xa+1
wa wa+1 w wa+2 ij+1−1
i a l a+ 1 l a+ 2j a a+1 ij+1 − 1 ij+1
Decrease by  capacity
Figure 3.5: Illustration for the proof of part (ii) from Theorem 3.7.
Suppose xa < xa+1 for some i j ≤ a < i j+1−1. Then we can take an appropriate amount
ε > 0 of capacity from xa+1,xa+2, . . . ,xi j+1−1 and distribute it uniformly to all of the other
edges, namely x1, . . .xa and xi+1, . . .xn without violating the monotonicity property. As a
consequence, θp for all p ∈ {a+1, . . . , i j+1−1} increase while the other values of θp for all
p∈/ {a+1, . . . , i j+1−1} decrease. If ε is appropriately chosen, z⋆ decreases, a contradiction.
Moreover, if θ1 < θi1 , then we can take an ε > 0 capacity from x1 and add it uniformly to all
other capacities so that θ1 ≤ θp and θp decreases for all p ≥ 2. Consequently z⋆ decreases,
which is also a contradiction. Therefore i1 = 1.
3.3 Capacity provisioning for an arbitrary path network with the sink
fixed at an endpoint
We first consider the path illustrated in Fig. 3.6, where we have three vertices with
evacuees and the sink is located at the right endpoint.
w1 l1 w2 l2 w3 l3
v x1 1 v x2 2 v x3 3 s
Figure 3.6: Example to illustrate longer path network with sink at the end.
67
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
Just as before, we need to minimize the following objective function
z(x) = max{θ1,θ2,θ3}, (3.34)
where
w
θ 11 = + τL[1]x1
W [1,2]
θ2 = + τL[2]x2
W [1,3]
θ3 = + τlx 3
.
3
Based on Theorem 3.2, we require x3 ≥ x2 ≥ x1.
We can formulate our problem as a convex, non-linear program, with variables Z,x1,x2,x3
and use a general purpose solver to find the optimal solution:
minZ, subject to:
Z ≥ w1 + τL[1]
x1
≥ W [1,2]Z + τL[2]
x2
≥ W [1,3]Z + τl
x 3
.
3
x1 + x2 + x3 = c
0 ≤ x1 ≤ x2 ≤ x3.
However, we can avoid solving this non-linear program by exploiting Theorem 3.7. We
propose a combinatorial algorithm for an arbitrary path with the sink fixed at an endpoint,
that will output an optimal parametric capacity assignment to the edges of the path. More
precisely, we obtain the solutions x∗2(x1), . . . ,x
∗
n(x1) as function of the parameter x1, the
capacity assigned to the first edge. We then determine the optimal value for x1 by solving
68
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
numerically the budget capacity equation,
x1 + x2(x1)+ . . .+ xn(x1) = c. (3.35)
Moreover, we also obtain specific conditions for Eq. (3.35), which allow us to determine x1
to the desired precision, by binary search.
According to Theorem 3.7, θ1 must be a candidate that determines the optimal cost. If
θ2 is also a candidate then we can write, after some basic algebraic manipulations,
⇒ 1 wθ θ 1 1 τl11 = 2 = = + . (3.36)x2 w1 +w2 x1 w1 +w2
We can substitute 1x = yi so that Eq. (3.36) becomesi
w
y 1
τl1
2 = y1 + . (3.37)w1 +w2 w1 +w2
Since Theorem 3.2 requires x1 ≤ x2, this implies y2 ≤ y1. Thus, we have
w1 τl τly 11 + ≤ y =⇒ y ≥ 1 . (3.38)w1 +w2 w
1 1
1 +w2 w2
The above implies x1 ≤ w2τl and as long as this is true, we can always use Eq. (3.36) to solve1
for x2. Therefore we have following:y1 if y < τl11 w
y 22 = (3.39)w1
w +w y1 +
τl1 τl1
1 2 w1+w
if y ≥
2 1 w2
We illustrate the parametric solution y2 in Fig. 3.7
69
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
yj
y1
y2
α = τ l1 y2 1w2
Figure 3.7: Graph illustrating the optimal computation of x2.
For example, let us recall the two nodes network discussed in Illustration 3.1.0.2. If y1
is in the region below α2, then we just need to solve for y1 in 2y = c. Otherwise, we need1
to solve for y1 in w1+w2 1 1 1y1(2w1+w2)+τl =1 c , since x1+x = c .2
We notice that y2 is a concave piece-wise linear function of y1. In general, we can
compute yi for all i ≥ 2 as a function of y1 using Theorems 3.7 and 3.2 as follows. We
consider two possible cases:
Case 1: If θi is dominated, that is θi < θ1, then xi = xi−1, according to Theorem 3.7.
Case 2: If θi = θ1, then θi is a candidate of the optimal cost. From the relation θi = θ1 and
with xi substituted by 1y , we obtaini
′ w1 L[1, i−1]yi = y + τ . (3.40)W 1 i 1[ , ] W [1, i]
70
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
Putting the two cases together, we have the following:
Lemma 3.8. Capacity for each edge ei for i ≥ 2 can be computed using
 w1 y + τL[1,i−1]W [1,i] 1 W [1,i] , if y 1
≥ αi
yi = (3.41)yi−1, otherwise,
where αi is the value of y1 for which yi = yi−1.
yj
p′ yi−1t
p′j
y′i
p′0
α′j+1
α′ . . . α′ α . . . α′ yi 11 j t
Figure 3.8: Explanation of how α is computed. We let α′ , . . . ,α′i 1 t be the breakpoints of yi−1
and p′0, . . . , p
′
t the corresponding linear pieces. Only the line segment [αi,∞) is stored in the
segment tree.
To construct an efficient data structure for the solution of our problem, we exploit the
property of the linear function derived in Case 2 above (see Eq. (3.40)). Let α′1, . . . ,α
′
t be
the breakpoints of yi−1 and let p′0, . . . , p
′
t be the corresponding linear pieces (see Fig. 3.8).
We denote by αi the value of y1 for which y′i = yi−1 (the intersection point between y′i and
yi−1). First, we argue that this intersection point is unique, because the slope of the linear
function y′, w1i W [1,i] is smaller than the slope of the last linear piece of y ,
w
− 1i 1 W [1,i−1] . This
property is crucial for the computation of all intersection points as well as the development
of our data structure as we shall see later in Subsection 3.3.1
71
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
Once we have the relation (3.41) for all yi with 2 ≤ i ≤ n, all we need is to numerically
solve the following equation in y1:
1 1 1
+ + . . .+ = c. (3.42)
y1 y2 yn
We now only have to discuss the procedure to compute αi and how we store the yi functions
so that Eq. (3.42) can be solved efficiently.
3.3.1 Implementation details and data structure
Suppose the piece-wise linear function yi−1 is known and its linear pieces are available
in order. We will compute yi from yi−1 using Eq. (3.41). To find the intersection point αi,
we first test the intersection with linear piece p′t (see Fig. 3.8). There are two cases. (a)
Function y′i intersects linear piece p
′
t at a point smaller than breakpoint α′t , and (b) y′i and
p′t intersect at a point greater than or equal to breakpoint α′t . Since the slope of function
y′i is smaller that that of p
′
t , case (a) means that y
′ < p′i t (y
′
i dominates yi−1) and thus p
′
t is
not part of function yi. We repeat the test with linear piece p′t−1, and so on, until we find
the linear piece p′ whose intersection with y′j i falls in case (b). That last intersection point
is breakpoint αi we needed to compute. We remark that computing breakpoint αi takes
amortized O(1) time since any linear piece whose intersection falls in case (a) is tested at
most once for all i ∈ {1, . . . ,n}.
The procedure illustrated above will compute the parametric solution yi for all i ∈
{2, . . . ,n} in total time O(n), however, we must be careful how we store this parametric
solution. The trivial approach of storing yi functions separately has O(n2) time and space
complexity. We propose to store the parametric solutions in a segment tree [28] (see Illus-
tration 3.3.1.1). The data structure has O(n) space complexity and can be easily constructed
in time O(n logn), as described next.
1) Compute all of the breakpoints αi, i = 2, . . .n associated with functions yi. We note
that these values may not be obtained in order. This step takes O(n) time as previously
72
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
explained.
2) To construct the segment tree, we sort the elements αi, say in increasing order. For each
element αi sorted, we make sure to retain the corresponding value of the index i. This
step takes O(n logn) time.
3) We construct the segment tree over the αi sorted elements, in O(n) time.
4) For each i ∈ {2, . . . ,n}, we insert interval [αi,∞) in the segment tree, in total O(n logn)
time (see Fig. 3.8 for an example).
We can now explain how we can use this segment tree to evaluate the sum 1 + 1y1 y + . . .+2
1
y at some fixed value y
′
1 = α , in O(logn+n) = O(n) time (see Eq. (3.42)). For simplicity,n
consider the example in Fig. 3.9 and Fig. 3.10.
In this example, the query value α′ is between breakpoints α4 and α3 and the vertical
query line through α′ intersects the intervals in the segment tree corresponding to indices
2, 4, and 7. This means that y2 and y3 are evaluated using linear function y′2; y4, y5, y6 are
evaluated using linear function y′4; and y7 is evaluated using linear function y
′
7.
In general, we describe the following procedure to evaluate the sum of reciprocals of
the parametric solutions yi:
a) Given query value α′, identify two consecutive breakpoints αi < α′ <= α j in the seg-
ment tree. Total time, O(logn).
b) Traverse the segment tree from the leaf node of α j to the root and collect the indices of
all intersected intervals. Total time: O(logn+n) = O(n). Let i1 < i2 < .. . < ik be these
indices in order.
c) We note that the list of indices of intersected intervals may not be obtained in the desired
order from the segment tree. We can sort these indices using counting sort in total time
O(n) [14].
73
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
d) If there are two or more indices: for any two consecutive indices i j and i j+1, evalu-
ate yi j , . . . ,yi j+1−1 using linear function y
′
i . If there are one or more indices: evaluatej
y1, . . .yi1−1 using y1. If i1 does not exist: then evaluate all yi parametric solutions using
y1.
Illustration 3.3.1.1. We show below (see Fig. 3.9) an example for arbitrary seven y′i func-
tions together with their breakpoints.
yj
y1 y2
y3
y4
y5
y
α 63
α6
α5
α2 α4 y7
α7
y1
Figure 3.9: Seven arbitrary y′i functions with their corresponding breakpoints. Only one-
sided line segments [αi,∞] are stored in the segment tree.
74
3.3. CAPACITY PROVISIONING FOR AN ARBITRARY PATH NETWORK WITH
THE SINK FIXED AT AN ENDPOINT
0 α7 α2 α4 α3 α5 α6 y′2(y1)
y′3(y1)
y′4(y1)
y′5(y1)
y′6(y1)
y′7(y1)
y1 = α
′
Figure 3.10: Stored line segments [αi,∞) for yi functions in Fig. 3.9. The vertical line
shows an arbitrary query point y = α′1 that cuts through only 3 line segments.
Fig. 3.10 shows the line segments stored in the segment tree together with a vertical line
that shows how an arbitrary query point y ′1 = α goes through just 3 line segments. From
Fig. 3.10, we notice that segment y′2 and y
′
3 share the same segment of y
′
2 and segments
y′ , y′4 5 and y
′
6 share the same line segment of y
′
4 based on the query point y1 = α
′. For this
example, we resolve Eq. (3.42) in the binary search as follow:
1 2 3 1
′ + ′ ′ + ′ ′ + ′ ′ ≤ (≥)c (3.43)α y2(α ) y4(α ) y7(α )
Using our segment tree, we can now proceed with binary search to identify a single
interval for y1 that contains the optimal value y∗1 in such a way that all yi functions are linear
in this interval. Thus, Eq. (3.42) consists of a sum of O(n) reciprocals of linear functions.
With the segment tree, this binary search procedure takes O(n logn) time.
We present the algorithm in Algorithm 5. In the algorithm, we let W be the sorted
vector of weights on the vertices and L denotes the vector of edge lengths.
75
3.4. CAPACITY PROVISIONING FOR EVACUATION OF AN ARBITRARY PATH
WITH ONE ARBITRARY SINK FIXED
Algorithm 5: Capacity allocation algorithm for an arbitrary path network with a
sink at its end
Result: Output the optimal cost and vector of edge capacities for an arbitrary path
network with the sink at its end
Inputs: W, L, c, τ;
Compute all breakpoints αi of y′i, for i = 2, . . .n and store them in a stack ;
Sort the n−1 breakpoints and construct a segment tree on them ;
Store only the needed n−1 intervals [αi,∞) at the nodes of the tree in the usual
way;
Execute Binary Search to find the interval I∗ for y∗1 and obtain the reciprocal
functions for yi ;
Solve Eq. 3.42 numerically by binary search over interval I∗ and obtain y∗1;
Output the solutions 1y (y∗i 1) ;
Output the cost z = θ1(x⋆1)
We now summarize the time complexity of algorithm 5 in the following theorem:
Theorem 3.9. The optimal capacities vector x⋆ for a path network of n+1 vertices with a
sink at its end can be obtained in O(n logn+n log(c/ξ)) time, where ξ is the precision for
parameter y∗1.
Proof. The first term in the complexity of the algorithm follows from the earlier discussion.
For the second term, we argue that each binary search iteration is resolved by evaluating
1
y′ +
1
y (y′) + · · ·+
1 ′ ∗
2 yn(y′)
≤ (≥)c for some value y in interval I . Alternatively, a numerical
scheme (e.g., Newton method) or binary search procedure can be employed to determine
y⋆ to some chosen precision ξ. Note that the range for I∗1 is of order of c, hence we can run
the search operation at a cost of O(log(c/ξ)) rounds. Since it takes O(n) to evaluate all the
linear functions yi, it follows that the whole operation takes O(n log(c/ξ)).
76
3.4. CAPACITY PROVISIONING FOR EVACUATION OF AN ARBITRARY PATH
WITH ONE ARBITRARY SINK FIXED
3.4 Capacity provisioning for evacuation of an arbitrary path with one
arbitrary sink fixed
We can extend Algorithm 5 to solve the general case. First, we compute, independently
of each other, the segment trees for each sub-path on either side of the sink. Let AL and AR
be the list of α values for the left and right sub-paths respectively. Let θL1(x1) and θ
R
1 (u1) be
the corresponding candidates of the costs for both subnetworks according to Theorem 3.7,
so that zL(X⋆) = θL ⋆1(x1) and zR(U
⋆) = θR ⋆1 (u1), where X
⋆ (resp. U⋆) is the sum of optimal
capacities for the left (resp. right) subpath. If we choose a query value x1 = 1y for the left1
subpath, we can easily solve for u1 using θL( 1 ) = θR( 11 y 1 u ). That is, for a fixed sink s, we1 1
have
1 wLy1 τ
= 1 LR + R (L [1]−L
R[1]), (3.44)
u1 w1 w1
where wL1 (resp. w
R
1 ) is the weight on the first vertex v
L
1 (resp. v
R
1 ) on the left (resp. right)
subpath, LL[1] = d(vL R R1 ,s) and L [1] = d(v1 ,s). Now, we can execute Algorithm 5 on the left
subpath with query 1y , compute the total capacity U of the right subpath using Eq. (3.41)1
and u1 given by Eq. (3.44). The addition of assigned capacities for both left and right is
then compared with c and the binary search procedure continues based on the comparison.
This procedure takes O(n logn+n log(c/ξ)) time.
77
3.4. CAPACITY PROVISIONING FOR EVACUATION OF AN ARBITRARY PATH
WITH ONE ARBITRARY SINK FIXED
θ1(
1 ) θ ( 1
y n−1 )1 u1
w x = 11 1 w2 w xn−1 =
1
y1 n−1 u1 wn
v1 l v2 ( ) S ( ) vn−1 ln−1 vn1
θ 11 y = θ
1
n−1
1 u1
X
X + U ≤ (≥)c U
Figure 3.11: Illustration of the application of Algorithm 5 to an arbitrary path network with
a fixed sink.
78
Chapter 4
Capacity provisioning on star networks
The vital idea we learnt from allocating capacity on path networks is that, if we can find the
optimal allocation for the furthest edge from the sink, we can find the optimal allocation for
all other edges in the network using the parametric equation given in (3.40). In this chapter,
we shall look at how we can extend this idea to star networks where the sink is located in
the middle of the star. We shall also investigate a complicated case where the sink is located
on a leaf node. For the latter case, we shall exploit the flow theorem idea of Klinz employed
in [10] to tackle the complicated case. We restrict our studies for the latter case to a star
network with three links (or edges), as more properties are required to tackle the general
case. The next section presents the case where the sink is located in the centre of the star
graph.
4.1 Star network with a sink in the middle
We investigate the problem of optimal allocation of capacities in star topology (see
Figure 4.1), where the sink is located at the centre. This can be linked in reality to a group
of scattered communities linked together by a star topology with individual edges linked to
a safe-haven at the centre
79
4.1. STAR NETWORK WITH A SINK IN THE MIDDLE
v v2 v13
v l1 x14 x4
l s4 xn
ln
v vn5
v6
Figure 4.1: Simple star graph with sink in the middle.
We observe that the monotonicity property in 3.2 does not apply here as each vertex has
a direct link to the sink. It is however easy to construct the model for the time functions
associated with the edges of the network. That is, the cost functions for all 1 ≤ i ≤ n is
given by
w
θ x ii( i) = + τli. (4.1)xi
Unlike in the path network studied earlier where the time function associated with the fur-
thest edge always dominate other cost functions at optimality, this case is different as given
below in Theorem 4.1
Theorem 4.1. When the capacity allocation is optimal, θ1(x1) = θ2(x2) = · · ·= θn(xn).
We present the proof below.
Proof. Suppose there is an optimal allocation for which θ j(x j)< θi for some indices i and j.
Then we can take an appropriately small amount ε > 0 from the capacity x j and distribute it
equally to the other edges of the star graph, thus decreasing θi(xi) and all θ functions other
than θ j. Therefore max1≤i≤n θi(xi) decreases and so the evacuation time also decreases,
which is a contradiction.
80
4.2. STAR NETWORKS WITH MORE EDGES AND A SINK LOCATED AT THE
CENTRE
4.1.1 Algorithm for star networks with a sink in the middle
Without loss of generality, we consider the leaf nodes of the star network in Figure 4.1
labelled in non-increasing order of their edge length, so l1 ≥ l2 ≥ . . .≥ ln. From the relation
θi(xi) = θ1(x1), we have the following linear function:
w
y 1
l − l
i = y1 + τ
1 i
, ∀1 ≤ i ≤ n, (4.2)
wi wi
where we denote yi = 1x for all 1 ≤ i ≤ n. Notice that the second term in Eq.(4.2) isi
nonnegative and that is the reason for the edge length ordering. These equations describe
an optimal parametric solution to the capacity allocation problem when the sink is located
at the central vertex. We can obtain the optimal value for y1 using Eq. (3.42) and the
technique discussed in Subsection 3.3.1. We let L be the sorted edge length according to
the labelling described above, and W be the corresponding weight vector. We present the
algorithm below:
Algorithm 6: Capacity allocation algorithm for star networks (Figure 4.1) with a
sink at its centre.
Result: Output the optimal cost and vector of edge capacities for the star network
4.1
Inputs: W, L, c, τ;
Employ the technique discussed in Subsection 3.3.1 ;
Output the solutions 1yi(y∗) ;1
Output the cost z = θ1(x⋆1)
4.2 Star networks with more edges and a sink located at the centre
We extend our study to a more complex star network which have more edges in its links
but the sink is still fixed at the centre (see Figure 4.2)
81
4.2. STAR NETWORKS WITH MORE EDGES AND A SINK LOCATED AT THE
CENTRE
vm1
vm−1 11
xm−1 l
m
1 x
m
1 x
1 vn
n
1
lm−11 l
1
n
S
x21 x
2
n
v2
v2 l2 l2
n
1 1 n
xm−1n
x1 m−11 ln
lm xm vm−1l1 n n n1
v11
vmn
Figure 4.2: Regular star network with more edges and sink at the centre
To solve the above problem, we realise that the balancing algorithm discussed in Illus-
tration 3.1.0.3 must come into play. Let us assume that v11 is the furthest node to the sink
amongst all the leaf nodes, then at optimal solution, we require θ1 θ j1 = 1 for j = 2, . . . ,m.
We also understand that θ11 = θ
j
n for j = 2, . . . ,m at optimal solution. Therefore, we can
write
θ1 = θ21 1 = · · ·= θm1 = θ1n = θ2n = · · ·= θmn . (4.3)
More to the properties of the optimal solution, and precisely for the interior edges, we
understand that the monotonicity property (see 3.2) must be observed toward the sink. The
above equation helps us to employ the same idea of restricting the problem to finding the
optimal allocation to the furthest edge provided we can compute all other edge capacities
using this optimal value. Let y j 1i = j and from Eq.(4.3), we can write y
1
n,y
2
n, . . . ,y
m
x n
as
i
82
4.2. STAR NETWORKS WITH MORE EDGES AND A SINK LOCATED AT THE
CENTRE
functions of y1 2 m1,y1, . . . ,y1 . That is,
y1n = g
1
1(y1)
y2n = g2(y
2
1)
(4.4)
... · · · ...
ymn = gm(y
m
1 ).
With θ11 being the function associated with the furthest edge to the sink, we can compute
x j j1 =
1
j and xn = 1j , j = 2,3, . . . ,m as functions of y1y y 1. That is1 n
y21 = h (x
1
1 1)
y31 = h2(x
1
1)
(4.5)
.
.. · · · ...
ym = h (y11 m−1 1)
and
y2 1n = q1(x1)
y3 1n = q2(x1)
(4.6)
.
.. · · · ...
ymn = q
1
m−1(y1).
Note also that the interior edge capacities y12, . . . ,y
1 1
n−1 can be written as functions of y1,
y22, . . . ,y
2 2 m m
n−1 as functions of y1 and y2 , . . . ,yn−1 as functions of y
m
1 . That is,
y12 = r
1
1(y
1
1) y
2 2 1 m m m
2 = r1(y1) · · ·y2 = r1 (y1 )
y1 = r1(y1 2 2 23 2 1) y3 = r2(y1) · · · ym = rm(ym3 2 1 )
(4.7)
... · · · ...
y1 1 2 2 2 m m mn−1 = rn−2(y1) yn−1 = rn−2(y1) · · · yn−1 = rn−2(y1 ).
83
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
Thus all unknown edge capacities can be written as functions of the furthest edge capacity
y11 and we can write
y1 +h (y11 1 1)+ · · ·+hm−1(y1 1 11)+q1(y1)+ · · ·+qm−1(y1)
+ r1(y11 1)+ · · ·+ r1 (y1)+ r2n−1 1 1(h1(y11))+ · · ·+ r2n−2(h1(y11))
+ · · ·+ r2n−2(h (y11 1))+ · · ·+ rm1 (hn−1(y11))+ · · ·+ rm 1 m 1n−2(hn−1(x1))+ · · ·+ rn−2(hm−1(x1)) = c,
(4.8)
and the binary search of the previous chapter (see 3.3.1) can employed.
4.3 Three edged star network with the sink located on a leaf
We extend our study of capacity provisioning on star networks, where in this case the
sink is located on a leaf node (see Figure 4.3)
S
l0 x0
w0
v0
l1
x1
w l21 x2
v1 w2
v2
Figure 4.3: Three-edged tree network with a sink at its leaf
This case is more difficult compared to the cases of the star networks considered earlier
as it was easy to compute all other edge capacities in term of the furthest edge capacity, and
the binary search technique discussed in Subsection 3.3.1 can be employed to search for
the optimal value of the furthest edge capacity that solves the allocation problem and such
that the addition of the capacity is equal to the budgeted capacity c.
Before we proceed, we present the following theorem called the Klinz theorem em-
84
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
ployed in [10] for the quickest transshipment problem. This is because the evacuation
problem is a special case of the quickest transshipment problem [10, 25]. The theorem
gives us an idea of computing the optimal cost of flow:
Theorem 4.2 (Klinz theorem [10, 26]). Let S be the set of terminals (set of supply sources
and sink) and let A be a subset of S. Let o(A) denote the maximum amount of flow that the
sources (of supply) in A can send to the sink in S \A in time θ without considering other
terminals and let v(A) be the total supply in A. Then an evacuation problem is feasible if
and only if o(A)≥ v(A).
Therefore given a dynamic tree network with evacuees on vertices, the optimal cost of
flow is determined by the largest evacuation time when all the subsets that constitute the
powerset of evacuees vertices are considered. From the tree network above, the powerset
consists of 7 elements given by {v0,v1,v2,{v1,v2},{v1,v0},{v2,v0},{v1,v2,v0}}. How-
ever, some of them can be discarded as they are dominated by obvious ones. We note
that the flow from {v1,v0},{v2,v0} and v0 are dominated by the flow from v1,v2 and
{v1,v2,v0} = v012 respectively. Therefore, we need to consider the flow dictated by the
following subsets {v1,v2,v12,v012} and the one with the maximum evacuation time is taken
as the optimal cost. In addition to the property of the optimal solution and cost, the follow-
ing lemma emphasizes a monotone allocation from the edges associated with the leaves to
the sink:
Lemma 4.3. Capacity allocation is monotone non-decreasing from the leave edges to the
sink.
Proof. The flow from v1 and v2 to the sink s can be seen as two path networks and therefore
the same argument of Theorem 3.2 of a path network applies.
Let us denote by f1 (resp. f2) the amount of flow from v1 (resp. v2) that the network can
allow to the sink within time t. According to Klinz theorem, we let θ1 (resp. θ2) to denote
the evacuation time that matches up with the flow from v1 (resp. v2) and the corresponding
85
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
unknown edge capacity is denoted by x1 (resp. x2). Let us assume that v1 is the furthest
to the sink. The flow dictated by {v1,v2} denoted by f1,2 starts with capacity x2 because
the first evacuee arrives from v2 and at τ(l0 + l1), we have evacuees arriving from v1 and
thus more flow. However, the rate of this flow is limited by the capacity x0 of the link e0
or x1 + x2. Using Klinz theorem, we look at the weights w1 and w2 and check the time for
which the flow f12 equals w1 +w2. We denote this evacuation time by θ12.
4.3.1 Algorithm strategy
We note by Lemma 4.3 that c3 ≤ x0 < c. We fix a value for x0 in its range so that x1+x2 =
c− x0. Then we look for the best way to determine x1 and x2 such that max{θ1,θ2,θ12} is
minimized. We then use this approach as a subroutine to determine x0 in such a way that
the optimal solution or cost is attained.
Observation. If x0 is large such that x0 > x1 + x2 then x c0 > 2 . Then the best assignment
is the one in which θ0 = θ1 = θ2 = θ12. Any solution that deviates from this equal cost
will make one of the evacuation time higher and such evacuation time can be reduced by
adjusting the capacities appropriately.
Note that it is possible to have x0 > c2 because of the possibility of a large weight at v0
and it is reasonable to assign a large capacity to the link e0 so that many evacuees can be
evacuated and also reduce congestion at v0.
Case 1: Based on the observation if x0 is fixed in [ c2 ,c), then x1 + x2 is constant. Let θ
be the time for which all evacuation time θ2 = θ1 = θ12. Suppose we have θ1 and θ2 such
that θ does not exist, we can consider the addition of x1 and x2 and construct the solution
that gives θ.
86
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
W
w1 + w
f12
2
f
w 22
f1
w1
τ(l0 + l2) τ(l0 + l1) θ2 = θ1 = θ12 = θ
t
Figure 4.4: Graph illustrating the case when x c0 ∈ (2 ,c)
Let X = c−x0, where x1+x2 = X , then from the above graph (from flow f1 and f2) we
have
w1
= x
θ− τ l l 1( 0 + 1) (4.9)
w2
= x .
θ− τ l 2( 0 + l2)
Then adding the two equations above gives
w1 w2
+ = X . (4.10)
θ− τ(l0 + l1) θ− τ(l0 + l2)
From the above we can solve for θ in term of X to have θ = g(X) = g(c− x0). We only
need to compare g(c− x0) with θ0 which gives us the idea of binary search to find x0 such
that
w1 +w2 +wθ 30 = + τl0 = g(c− x0). (4.11)x0
With x0 fixed, we can compare θ0 with g(c−x0) and if θ0 > g(c−x0) then we can increase
x0 otherwise we decrease x0. If x0 =X and θ0 < g(c−x0), then we have the optimal solution
87
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
already because we cannot increase or decrease x0 anymore.
Case 2: Now suppose x ∈ ( c , c0 3 2) and we still maintain x1+x2 =X . We note that X ∈ (
c , 2c2 3 )
and θ1 dominates θ2 and θ12. This is because if θ2 or θ12 dominates, x1 can be reduce by
ε > 0 capacity which increases x2. We can continue to decrease x1 until θ1 dominates. We
must also remember that the value assigned to x1 and x2 are limited by x0. Let us examine
the dynamic of the flow f12 in relation to the flow f2 using the graph below
W
w1 + w
f12
2
f2
w2 f1
WX
τ(l0 + l2) τ(l0 + l
t
1)
Figure 4.5: Graph illustrating the case when x0 ∈ ( c3 ,
c
2). Note that the slope of f12 is not
smaller than the larger of x1 and x2 but limited by x0.
The slope of f12 is not smaller than the larger of the slope of f1 and f2 but not greater
than x0. If we increase x1 and decrease x2 then θ1 decreases while θ2 increases and also θ12
increases because the parallel line f2 decreases. Now if x2 is increased to its upper bound
c− x0 then θ12 > θ2, however, x2 must not be greater than x0. Geometrically we derive θ12
as follow: From the line f2 and f1 we have
w1
= x =⇒ w = x τ(l − l ). (4.12)
τ(l1 − l
2 1 2 1 2
2)
88
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
Considering the starting point of line f12 and its slope we have
w2 +w1 −w1
= x =⇒ w = θ
θ − τ l l 1 2 12
x1 − x1τ(l0 + l1). (4.13)
12 ( 0 + 1)
Adding Eq. (4.12) and Eq.(4.13) and solving for θ12 gives
w1 +w2 − x2τ(l1 − l2)θ12 = + τ(l0 + l1). (4.14)x1
When x2 is closer to zero, θ2 dominates θ12 and therefore there is a point for x2 where
θ2 = θ12. Therefore if we can determine the point x2 such that θ12 = θ2, then we can
compare with θ1. If θ1 is larger then we can decrease x2 and from there we know that
θ2 > θ12 and we can always compare θ1 and θ2. On the other hand if θ1 < θ12 = θ2 at the
critical point x2, then we can increase x2 and we know that θ12 > θ2. Therefore we can
compare θ1 and θ12 to search for the optimal solution by binary search by comparing θ1
and θ12. When the best x1 and x2 are computed after the binary search, the maximum cost
dictated by θ1 and θ2 or θ2 and θ12 is compared with θ0 using binary search.
Now let us look at the following graph to analyse the weight wX such that θ2 = θ12
given x0 fixed.
89
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
W
Slope x0
w + w θ1,21 2
θ2
w2
WX
τ(l0 + l2) τ(l0 + l )
t
1
Figure 4.6: Graph showing the critical weight WX at which θ2 = θ12.
Using the slope of f2 together with the similar triangle in the graph we have
w2 wX
= . (4.15)
θ2 − τ(l0 + l2) τ(l1 − l2)
Using the line f12 we have
w1 +w2 −wX
= x . (4.16)
θ2 − τ
0
(l0 + l1)
With Eq. (4.15) and Eq. (4.16) we have a quadratic equation in θ2 when wX is eliminated:
x θ20 2 −θ2(x0τ(l0 + l2)+ x0τ(l0 + l1)+(w1 +w2))
(4.17)
+ x0τ2(l0 + l1)(l0 + l2)+(w1 +w2)τ(l0 + l2)+w1τ(l1 − l2) = 0.
Once θ2 is computed, x2 can be computed and we know that θ2 = θ12 at the computed x2.
Then we can compare with θ1 to know which direction to go. Thus the quadratic equation
gives us the critical value of θ2 (or x2) before the comparison (by binary search) procedure.
The case of an extended star network with the sink located at a leaf node is more com-
90
4.3. THREE EDGED STAR NETWORK WITH THE SINK LOCATED ON A LEAF
plicated as there are many flow of the form fi j whose influence in the set of optimal cost
need to be tracked or considered.
91
Chapter 5
Minsum capacity provisioning on path
networks
In this chapter, we turn attention to capacity provisioning problem on path networks, with
the objective of minimizing the total evacuation time of all evacuees in the topology. We
call this problem the minsum capacity provisioning problem in relation to the traditional
minsum sink location problem [24, 37]. Recall that the minsum sink location is more
difficult in any topology of interest compared to the minmax problem (see Section 1.2),
and in fact the study of the problem is limited to path topology in literature (see the survey
paper [25]). However, the minsum problems are more attractive when considering the
psychological duress that may be suffered by the evacuees on their way to safe-haven.
As we have learnt in Chapter 3, the traditional idea of sink location algorithms is based
on the possible evacuation by foot, which is effective in when considering dense commu-
nities. As for remote or sparse communities, evacuation procedure surely require some
adjustments to the traditional model. Just as in in Chapter 3, we shall view the budgeted
capacity as resource to be allocated optimally to the edges of the network in a way that the
evacuation total time is minimized, and we require the sum of the allocated capacities to
be no more than the budgeted capacity. The next section present the former definition of
the problem, as well as the model formulation. We shall spend more time in the section to
properly investigate and establish the problem with some examples.
92
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
5.1 Preliminaries and model formulation for the minsum capacity pro-
visioning concept
Let P = (V,E) represents a path graph with a vertex set V = {v1,v2, . . . ,vn+1} and an
edge set E = {ei = (vi,vi+1), i = 1, . . . ,n} . Each vertex vi has a non-negative weight wi
which represents the number of evacuees situated at vi, and each edge ei has a fixed length
li but unknown capacity xi, which represents the upper limit on the number of evacuees
that can be evacuated through the edge ei in a unit time. The parameter c denotes a given
budgeted capacity to be shared optimally between all the edge lengths and τ denotes the
transit time per unit distance. We define the minsum capacity provisioning formerly as
follow:
Problem 5.1. Given a dynamic path network N = (P,c,τ,w), a fixed sink S ∈ P, and a total
capacity budget c. Let z(x) denote the sum of the evacuation times of all evacuees on P
given the capacity vector x = (x1,x2, . . . ,xn) where xi represents the capacity assigned to
edge ei on the path. Find the optimal capacity vector x⋆ such that z(x⋆) is minimum and
∑ni=1 xi = c.
Let us investigate the problem on the following simple two-edged network (Figure 5.1),
where the sink is fixed at the right end of the network and evacuees from v1 and v2 have to
move at the same time to the sink. We let x2 = x so that x1 = c− x.
w1 l1 w2 l2
v1 x1 v2 x2 S
Figure 5.1: Illustration of capacities allocation on a path network with two edges. The sink
is assumed to be located at the end vertex on the right.
There are two cases to consider:
[Case 1 (No Congestion):] Let us assume that there is no congestion. That is, the first
set of evacuees from v1 gets to v2 and moves forward to the sink without any delay at v2.
93
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
Mathematically, this implies that τl1 ≥ w2x or x ≥
w2
τl . We can derive the objective function1
by adding the areas of the two trapeziums in Figure 5.2. Thus we have
w2 w2
z(x) = 1 + 2 +w τ(l + l )+ τl w . (5.1)
2(c− x) 2x 1 1 2 2 2
Figure 5.2: Time-weight graph of Figure 5.1. There is no congestion and evacuees get to
the sink with their starting rates
From the objective function, if x < c2 then the cost becomes
w2 +w2
z(x) = 1 2 + τl
2x 2
w2 + τ(l1 + l2)w1, (5.2)
which we can easily improve by allocating equal capacity to both edges. Thus trivially, the
94
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
cost
w2 +w2
z(x) = 1 2 + τl2w2 + τ(l1 + lc 2
)w1 (5.3)
is better (lesser) than the cost in (5.2) since 2x < c. Can we say we have reached the optimal
cost and solution with the above? This will be clarified with the analysis below.
First, we observe that the objective function can be re-written as
z(x) = w1τ(l1 + l2)+ τl2w2 + min g(x), (5.4)
x∈[ c2 ,c)
where
w2 w2
g(x) = 1 + 2 . (5.5)
2(c− x) 2x
Since the first two terms of (5.4) are constant, only the function g(x) is needed to be anal-
ysed. We also observe that the function g(x) is convex for x ∈ [c/2,c] and hence must
have one local minimum (correspondingly global minimum in the specified interval). To
check this, we use the second derivative test for the roots of g′(x) = 0 as follow. Computing
g′(x) = 0 leads to the following quadratic equation
x2(w2 2 21 −w2)+w2cx−w22c2 = 0, (5.6)
whose roots are √
−2w2 4 2 2 22c 4w2c +4(w1 −w2)w2c2x 2a,b = ± . (5.7)2(w2 2 2 2 21 −w2) 4(w1 −w2)
Further simplification yields
[ −w2c w ]x x 2ca = , b = . (5.8)w1 −w2 w1 +w2
Now using the second derivative test to examine the nature of each root derived, we first
95
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
discover that
′′ w
2 w2
g (x) = 1 23 + 3 (5.9)(c− x) x
is positive for any value of x in [c/2,c). Now for the first root, we have
′′ −(w 4g x 1 −w2)( a) = 3 , (5.10)c w1w2
which is obviously negative and hence xa corresponds to a local maximum of g(x). The
second test gives
′′ (w1 +w
4
g x 2
)
( b) = 3 , (5.11)c w1w2
which is non-negative and therefore xb corresponds to a local minimum. Furthermore, since
g(x) is convex, xb also corresponds to its global minimum point (and global minimum of z)
in the interval [c/2,c]. { }
Therefore for the problem at hand if c ≥ 2w2τl , the optimal solution is max
c w2c
1 2
, w1+w .2
[Case 2:] Let us now consider the case when c < 2w2τl . In this case, the optimal solution1
may be in congestion case or out of congestion as we shall soon find out. First let us draw
the graph of the scenario.
T
θ1(x,w)
θ2(x,w)
τ(l2 + l1)
τ l2
0
w2 η2,1 w1 + w2 W
Figure 5.3: Congestion scenario for the path network in Figure 5.1. The intersection of the
two straight lines indicate the time at which θ2(x2,η) = θ1(x1,η)
96
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
At x = c/2, the two lines θ1 and θ2 are parallel since they both have the same slope
and the cost can be computed using the dominant of the two functions. However, the cost
can be improved by increasing x to a point where the two lines intersect as shown in Fig-
ure 5.3. Moreover, there is a possibility that the optimal solution vector is in non-congestion
configuration, which is case 1.From the graph we have the following time functions which
represent the time for the last evacuees at v2 and v1 to get to the sink:
w
θ2(x,w) = τl2 + , w ≥ 0x
w−w (5.12)
θ x w 21( , ) = τ(l1 + l2)+ , w > w ,min{c− x,x} 2
however because of monotonicity property (which will be discussed later in Section 5.2) of
the optimal solution and the geometric picture above, we define θ1 as follow:
w−w
θ 21(x,w) = τ(l1 + l2)+ , w > w2. (5.13)c− x
To find the value of η as a function of x, we need to minimize the following time
Θ(x,η) = max{θ1(x,η),θ2(x,η)}, (5.14)
that is, we need to solve for η in
η−w η
τ l l 2( 1 + 2)+ − τlc− x 2
− = 0. (5.15)
x
Solving the above gives
w2x− τl1(c− x)xη = . (5.16)
2x− c
Now equation (5.16) also verify that the two lines are parallel at x = c2 as η does not exist at
this point. Using Figure 5.3, we derive the objective function geometrically by adding the
97
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
areas of the two trapeziums featured in the graph. We have
1 1
z(x) = (θ2(x,0)+θ2(x,η))η+ (θ2 2 1
(x,η)+θ1(x,w1 +w2))(w1 +w2 −η), (5.17)
where we have put x1 = c− x and x = x. Expanding (5.17) with η defined in (5.16) gives( 2 2 )
1 2τl2 +
w2x−τl1(c−x)x
( (2x−c)x
(w2x− τl1(c− x)x)
z(x) = +
2 2x− c )
1 w2x− τl (c− x)x w
2( τl 12 + + τ)(l1 + l 12)+ × (5.18)(2x− c)x c− x
w
w w 2
x− τl1(c− x)x
1 + 2 − .2x− c
Now we compute z′(x) = 0 to have
z′(x) =−2τ2l2x4 +6cτ2l2x3 +(−7c2τ2l2 +2cτw l +4w2 −2w21 1 1 2 1 1 2)x2
+(4c3τ2l2 −4c21 τw2l1 −4cw21 +4cw22)x− (5.19)
− c4τ2l2 3 2 2 2 21 +2c τw2l1 + c w1 −2c w2 = 0.
Compactly,
Ax4 +Bx3 +Cx2 +Dx+E = 0, (5.20)
where A,B,C,D,E are the coefficients in the polynomial (5.19).
The roots of (5.20) that lies in ( c2 ,c) can only be found numerically as the closed form
analytical method does not exist. However, we note that if z(x) is convex in the interval
(c/2,x), then one of the roots of z′(x) (which corresponds to the local minimum) lies in the
interval.
Lemma 5.2. Given a path network with two edges and a sink located at its end. The
minsum capacity objective functions associated with the network for both congestion and
non-congestion cases are convex.
98
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
Proof. Note that monotonicity implies x2 ≥ x1. The Hessian matrix associated with the
non-congestion case is  2HNon = w13 0x1  (5.21)2
0 w2
x3
.
2
The determinant of k× k, k = 1,2 submatrix of HNon are non-negative, hence H is positive
definite. For the congestion case, we have the determinant of its 1×1 matrix to be
l2τ2x3x21 1 2 −2l1τw2x3x 2 3 2 2 2 2 21 2 −w1x1 +3w1x1x2 −3w1x1x2 +w1x3 +w22 2x31 . (5.22)
x31(x2 − x1)3
The denominator is non-negative since x2 ≥ x1. The numerator simplifies to
x3 2 2 31(w2 − l1τx2) +w1(x2 − x1) (5.23)
and also non-negative. The determinant of its 2by2 submatrix gives
w2 21(w2 − τl1x1)
x3
(5.24)
1(x2 − x1)3
which is non-negative. Therefore the Hessian matrices associated with the two-edged net-
work are positive definite.
Let us further examine the problem with some examples.
Example 5.3. Let the budgeted capacity c be equal to 10, weights w1 and w2 be equal to
10 and 18 respectively. Let the transit time per unit length τ be equal to 1, and edge lengths
l1 and l2 be 3 and 4 respectively. Clearly, c < 2w2τl . If we directly assume that the optimal1
solution would be in congestion case, then the objective function computed from (5.17)
99
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
gives
1 x(−12+3x)2 4x(−12+3x)
z =
2( +(2x−10)2 2x−10 (5.25)1 − x(−12+3x))( −12+3x 10 )+ 28 11+ + .
2 2x−10 2x−10 10− x
Now z′(x) = 0 gives
−18x4 +540x3 −5468x2 +23360x−36800 = 0 (5.26)
The zero of (5.26) that lies in (5,10) is a single value x = 6.2825. The second derivative test
shows that this value corresponds to a local minimum. The cost at this point is z(6.2825) =
181.3187. If the optimal solution is in a congestion scenario, then the value of η at this point
must be between w2 = 18 and w1 +w2 = 28. The value of η at this point gives 16.7760
which contradicts our analysis. In fact if the optimal solution lies in a congestion case as we
have assumed, then this point must lie between xmin and xmax which correspond to when η
is equal to ηmax = w1 +w2 and ηmin = w2 respectively. To further verify this, we compute
xmin as follow:
w2x− τl1(c− x)xw1 +w2 = 2x− c
18x−3(10− x)x
28 = (5.27)
2x−10
0 = 3x2 −68x+280.
Solution to the quadratic equation above is x = 17.25 or x = 5.40. Since the first solution is
above the budgeted capacity, we can discard it, hence xmin = 5.40. To compute the cost at
xmin, we need the following objective function
η2
z(xmin,η ) = τl η + maxmax 2 max (5.28)2xmin
100
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
which is derived when the two lines intersect at w1 + w2. Thus the cost at xmin gives
z(xmin) = 184.5925. In similar way, we can compute xmax as follow:
w2x−3(c− x)xw2 = 2x− c
18x−30x+3x2
18 = (5.29)
2x−10
0 = 3x2 −48x+180.
Solving the quadratic equation above gives x = 10 or x = 6. We discard the first as it is
equal to the budgeted capacity and xmax = 6. We can derive the corresponding objective
function from Figure 5.4 by adding the areas of the two trapeziums in the figure. Thus we
have
( ) ( )
1 η 1 η w
z(xmax,ηmin) = 2τl
min
2 + w2 + 2τl2 + τl
min 1
1 + + w1. (5.30)2 xmax 2 xmax c− xmax
T
θ1(x,w)
θ2(x,w)
τ(l1 + l2)
τ l2
w2 = η w1 + w2 W
Figure 5.4: The case when η = ηmin which corresponds to maximum value of x if the
optimal solution is in congestion scenario
101
5.1. PRELIMINARIES AND MODEL FORMULATION FOR THE MINSUM
CAPACITY PROVISIONING CONCEPT
Now at xmax = 6, z(xmax,w2) = 181.50. We discover that the cost derived when x =
6.2825 is better than the cost when x is equal to xmin and xmax. However, this value is not
between xmin and xmax. Note that we have congestion when x ∈ (5,6]. Can we find a value
of x in this interval that gives a better cost? This means we need to minimize
( )( )
1 12−3x 10 12x−3x2
z(x) = ( 4+ +7+ 28−2 10−2x )(10− x ) 10−2x (5.31)
1 12−3x 12x−3x2
+ 4+4+
2 10−2x 10−2x
with x ∈ [5,6]. Unfortunately the equation Z′(x) = 0 has no zero in [5,6]. Now as we
increase x beyond 6, we move away from congestion and the line θ1(x) goes beyond the
line θ2(x). The objective function in this case is (the no congestion case)
1 1
z(x) = (τl2 + f2(x,w2))w2 + (τ(l1 + l2)+ f1(x,w1 +w2))(w1 +w2 −w2 2 2
)
162 50 (5.32)
= 142+ + .
x 10− x
Minimizing the above function with x ∈ [6,10] gives x = 6.4287 with the cost 181.20. Since
this solution gives a better cost compared to others derived previously, we can write (the
optimal solution vector) x⋆ = [3.57,6.43].
The above example shows that the initial situation or configuration may be a congestion
scenario but the optimal solution may not. We can easily verify the no congestion scenario
by checking the inequality τl w21 > x⋆ . The next example shows that the optimal solution may
remain in congestion configuration.
Example 5.4. Let c= 20, τ= 1, w1 = 100, w2 = 125, l1 = 10 and l2 = 8. Then the objective
function gives
−25 4x
3 −708x2 +19215x−125100
z(x) = . (5.33)
4 (x−10)(−20+ x)
102
5.2. PROPERTIES OF THE OPTIMAL SOLUTION VECTOR AND COST
Computing z′(x) = 0 leads to
−200x4 +12000x3 −221250x2 +1650000x−4500000 = 0. (5.34)
The zero of the polynomial (5.34) that lies in (10,20) is derived uniquely (numerically)
as x = 11.7843 with the cost z(x = 11.784) = 4076.5527. The value of η at this point
is 141.4764. Since w2 < η(x) < w1 +w2, the solution is in congestion case. To further
verify that the above solution is the best, we compute xmin = 10.7884 (corresponding to
ηmax = 225) and xmax = 12.5 (corresponding to ηmin = 125 ). Now z(xmax) = 4091.6666
and z(xmin) = 4146.2801. Since x = 11.784 gives the best cost, we write x⋆ = 11.784. We
do not need to search for solution in [12.5,20] because the corresponding ”no congestion”
objective function has no zero in the specified interval.
We present the algorithm for the two-edged minsum capacity provisioning below. Note
that from the analysis above, the candidate solution for x2 are x = chalf 2 , x
w2c
b = w1+w ,2
x ax = w2m τl and the root xnum of the order 4 polynomial (5.19) that lies in (c/2,x1 max).
Algorithm 7: Capacity allocation algorithm for two-edged network with sink lo-
cated at its end.
Result: Output optimal cost and vector of edge capacities for the network in
Fig. 5.1
Inputs: w1, w2, c, τ, l1, l2;
Compute the candidate vector for x2 that is, x2 = {xhalf,xb,xmax,xnum} ;
Compute the corresponding x1, that is x1 = c ·1− x2 ;
Compute the corresponding values of objective functions z(x) ;
Output the vector x = [x1,x2] that gives the minimum cost.
5.2 Properties of the optimal solution vector and cost
From the classical sink location problem we understand that the sink must be on a vertex
(see [24] and Section 1.2), we repeat the property here as a lemma:
103
5.2. PROPERTIES OF THE OPTIMAL SOLUTION VECTOR AND COST
Lemma 5.5. The optimal location of the sink s must be a vertex.
Now for capacity allocation problem, we have the following properties:
Theorem 5.6. Given a path P with n+ 1 vertices where v1, . . .vn are occupied by evac-
uees w , . . . ,w and the sink is located at v . Let x⋆ = (x⋆1 n n+1 1, . . . ,x
⋆
n) be the optimal edge-
capacities vector and z(x⋆) denotes the total evacuation time, then x⋆ ≥ x⋆i j if i > j.
Proof. Let us first assume that the flow is congestion free. WLG let us consider flow from
indexes vi and vi+1 We have the following evacuation times at vi and vi+1 respectively:
w
θ ii = + τd(v ,s)min{xi,xi+1}
i
(5.35)
w
θ i+1i+1 = + τd(vi+1,s),xi+1
and individual contributions to z(x⋆) gives
i w
2
z = i + τd(v ,s)w
2min{xi+1,xi}
i i
2 (5.36)
zi+1
w
= i+1 + τd(v
2x i+1
,s)wi+1.
i+1
Since we can reduce zi+1 ( or overall reduce z(x⋆)) by making x ⋆ ⋆i+1 ≥ xi then xi+1 > xi in
x⋆.
Now suppose there is a congestion at vi+1, where some evacuees from vi encounter a
delay at vi+1 due to backlog of evacuees waiting to depart vi+1, then there exist some supply
ηi+1 that constitutes part of supply from vi that merge with some supply at vi+1. Therefore
we have the following time functions:
w
θ ii = + τd(vmin{x x } i
,s)
i+1, i (5.37)
η
θ i+1i+1 = + τd(vi+1,s),xi+1
104
5.3. EXTENSION TO LONGER PATH NETWORKS
and contribution to z(x⋆) gives
2
zi,i+1
η
= (i+1 + τd(v2x i+1,s)ηi+1i+1 ) (5.38)
1 ηi+1 w
+ τl ii +2τd(vi+1,s)+ + (wi+1 +wi −η2 x i+1
)
i+1 min{xi+1,xi}
Obviously the above contributions can be reduced by making x ⋆i+1 greater than xi in x .
Theorem 5.7. The minsum objective function z for the capacity provisioning problem is
convex.
Proof. Since the components of z can be broken down to set of convex functions, therefore
their addition is convex.
5.3 Extension to longer path networks
We have seen in the above section that the objective function for the minsum capacity
provisioning is convex. However, it configuration could change during the optimization
process. For simplicity, we shall focus on the objective function that defines the flow to the
right zLx and assume that the objective function for flow to the left is zero, i.e., zR(x) = 0.
The flow to the right can be analysed symmetrically. Thus, if the sink is located on any
interior vertex of the path, then
z(x) = zL(x)+ zR(x). (5.39)
We also define for j > i, W [i, j] = ∑ jk=i wk, L[i
j
, j] = ∑k=i lk and for j ≥ 1, W [ j
j
] = ∑k=1 wk
and L[ j] = ∑ jk=1 l j. Let us consider the three edged network shown below:
w1 x1 w2 x2 w3 x3
v1 l1 v2 l2 v3 l3 S
Figure 5.5: Path network with three edges and sink located on the right end.
Four possible initial configurations for z can be derived for this network:
105
5.3. EXTENSION TO LONGER PATH NETWORKS
Case 1: When we have congestion free scenario, that is, η1,2 ≤ w3 and η2,1 ≤ w3 +w2,
then the objective function is
1 1
z(x) = (τl3 +θ3(x3,w3))w3 + (τL[2,3]+θ (x ,W [2,3]))w2 2 2 2 2 (5.40)
1
+ (τL[1,3]+θ1(x1,W [1,3]))w1.2
T
θ1
τL[1, 3]
θ2
τL[2, 3]
θ3
τ l3
w3 W [2, 3] W [1, 3] W
Figure 5.6: Congestion free scenario where η3,2 ≤ w3 and η2,1 ≤W [2,3]
Case 2: When we have Congestion only at v3 but not at v2. That is, w3 < η3,2 ≤ W [2,3]
and η2,1 ≤W [2,3].
106
5.3. EXTENSION TO LONGER PATH NETWORKS
T
θ1
τL[1, 3]
θ2
θ3
τL[2, 3]
τ l3
w3 η3,2 W [2, 3] W [1, 3] W
Figure 5.7: Congestion only at v3
Note that if η3,2 <W [2,3], then we have
1 1
z(x) = (τl3 +θ3(x3,η3,2))η3,2 + (θ2(x2,η3,2)+θ2(x2,W [2,3]))(W [2,3]−η2 2 3,2
)
1
+ (τL[1,3]+θ1(x1,W [1,3]))w1.2
(5.41)
When η3,2 =W [2,3] and η2,1 ≤W [2,3], then
1 1
z(x) = (τl3 +θ3(x3,η3,2))η3,2 + (τL[1,3])+θ1(x1,W [1,3]))w1, (5.42)2 2
which is covered by (5.41).
Case 3: When we have congestion only at v2 and not at v3. That is, W [2,3]<η2,1 ≤W [1,3]
and η3,2 ≤ w3. If η2,1 <W [1,3], then
107
5.3. EXTENSION TO LONGER PATH NETWORKS
T
f1
τL[1, 3]
f2
τL[2, 3]
f3
τ l3
w3 W [2, 3] η2,1 W [1, 3] W
Figure 5.8: Congestion only at v2
1 1
z(x) = (τl3 +θ3(x3,w3))w3 + (τL[2,3]+θ2(x2,η2,1))(η2,1 −w )2 2 3 (5.43)
1
+ (θ2(x2,η2,1)+θ1(x1,W [1,3]))(W [1,3]−η2,1).2
When η2,1 =W [1,3], then
1 1
z(x) = (τl3 +θ3(x3,w3))w3 + (τL[2,3]+θ2(x2,η2 2 2,1
))(η2,1 −w3), (5.44)
which is covered by (5.43).
Case 4: If there is congestion at v3 and v2. That is, w3 < η3,2 ≤ W [2,3] and W [2,3] <
η2,1 ≤W [1,3]. Then we have
108
5.3. EXTENSION TO LONGER PATH NETWORKS
T
θ1
θ2
τL[3]
θ3
τL[2, 3]
τl3
w3 η3,2W [2, 3] η2,1 W [3] W
Figure 5.9: Congestion at both vertices
1 1
z(x) = (τl3 +θ3(x3,η3,2))η3,2 + (θ3(x3,η3,2)+θ2(x2,η2,1))(η2,1 −η3,2)2 2
1
+ (θ
2 2
(x2,η2,1)+θ1(x1,W [3]))(W [3]−η2,1).
(5.45)
The above objective function also covers the cases η3,2 =W [2,3], η2,1 =W [3].
Therefore in general, given any path network of n+1 vertices where the sink is located
on (n+ 1)th vertex, the minsum objective function for the capacity provisioning can be
109
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
written as:
1
z(x,β) = (θn(x2 n
,0)+θn(xn,βn,n−1))βn,n−1
1
+ (θn−1(xn−1,β2 n,n−1
)+θn−1(xn−1,βn−1,n−2))(βn−1,n−2 −βn,n−1)
+ · · · 1 1+ (θ2(x2,β3,2)+θ2(x2,β2,1))(β3,2 −β2,1)+ (θ1(x1,β2,1)+θ1(x1,W [n]))(W [n]−β2,1),2 2
(5.46)
where w if τln−1 ≥ wnn
β xn,n−1 = n (5.47)ηn,n−1, otherwise,
and for 2 ≤ j < n, we have
W [ j,n] if τl
w j
 j−1
≥
β x jj, j−1 = (5.48)η j, j−1, otherwise.
Note again that η j, j−1 is computed by solving for the weight in θ j(x j,w)−θ j−1(x j−1,w) =
0, for 2 ≤ j ≤ n. With the objective function defined, we need an algorithm that can ac-
commodate the changing nature of the function given an initial configuration. We also
need the algorithm to seamlessly accommodate the associated constraints, which include
the positive capacity allocation, monotone capacity allocation, and the sum of the allocated
capacity must be equal to the given budget c. We discuss our approach in the next section.
5.4 Solution by sequential quadratic programming
The sequential quadratic programming SQP is an iterative numerical method often em-
ployed to solve difficult optimization problems with constraints [19, 16]. The method incor-
porates the idea of Lagrange multiplier method to incorporate the constraints and then uses
the idea of quadratic programming to make the numerical method faster (see Chapter 2).
One important advantage of the SQP is that the initial solution vector can be made arbitrary.
110
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
However, because we understand that the optimal solution vector must be monotone, we
incorporate a monotone initial solution vector into the SQP methods.
For our problem, we focus on the flow from left to right zL(x) and assume zR(x) = 0.
We redefine the problem as follow:
min z(x)
S.t. x1 + x2 + · · ·+ xn = c (5.49)
xi − xi+1 ≤ 0, ∀1 ≤ i ≤ n−1. (5.50)
5.4.1 The SQP model and algorithm formulation
We employed the Newton SQP discussed in Chapter 2 . The Lagrangian objective
function which incorporate the constraints with the objective function is given by
Cineq
L = z(x)+λ0(x1 + x + · · ·+ x − c)+ ∑ λ (x − x 22 n j i i+1 + s j), (5.51)
j=1
where s j 1 ≤ j ≤ Cineq, (Cineq = 1+ (|E| − 2)) are slack variables added to convert the
inequality constraints to equality constraints, λ0 is the Lagrange multiplier for the equality
constraint in Eq.(5.49) and λ j, 1 ≤ j ≤Cineq are the Lagrange multipliers for the converted
inequality constraints.
Let ∆x= (x−xk,s− sk)′ and let g be the vector of constraints, then the Taylor expan-
sion of the Lagrangian function at the optimal solution gives
1
L(x⋆,λ⋆)+∇L(x⋆,λ⋆)∆x+ ∆x′B(x⋆,λ⋆)∆x, (5.52)
2
where we have included the Lagrangian multiplier λ0 in λ and B denotes the approximation
of the Hessian of the Lagrangian L . Now, linearizing the constraint vector at iterate xk gives
g(xk)+∇g(x ′k)∆x+ sk +∆ = 0. (5.53)
111
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
We now have the quadratic subproblem (see Chapter 2)
1
min ∇L(xk,λk)∆x+ ∆x′B(x2 k
,λk)∆x (5.54)
∆x
S.t. ∇g(x )′k ∆x+∇g(xk) = 0. (5.55)
The iterative method to be derive follows the same approach as the Newton’s method
as we shall see. Recall that the first order condition condition stipulates that ∇L(x,λ) = 0
(see Subsection 2.5.1). Differentiating the objective function in Eq.(5.54) with respect to
∆x gives
∇z+∇g(x ⋆ ′k)λ +∆x Bk = 0. (5.56)
The above equation can be rearranged to have
λ⋆
′
∇g(xk)+Bk∆x=−∇z. (5.57)
Setting λ⋆ = λk+1 = λk +∆λ and substituting in Eq.(5.57) give
−∇z−λ′k∇g(xk) =∆λ′∇g(xk)+Bk∆x
(5.58)
−∇L(xk,∆ ′k) =∆λ ∇g(xk)+Bk∆x.
Using Eq.(5.58) and the linearized constraint vector Eq.(5.53), we have the system
−∇L(xk,∆k) =∆λ′∇g(xk)+Bk∆x
(5.59)
−g(xk) = ∇g(x ′k)∆x,
which can be arranged in matrix form as

B(xk,λk) JTg (xk)


  
∆x ∇z(x ′k)+λk∇g(xk)=−  (5.60)
Jg(xk) 0 ∆λ g(xk)
112
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
and compactly as 
J 

L ∆x=−∇x,λL(x,λ), (5.61)
∆λ
where Jg denotes the Jacobian of the constraint vector g and ∇x,λ denotes the Jacobian of
L with respect to x and λ. The matrix on the right is the KKT or Jacobian matrix for the
problem. The one-step iterates are then computed using:
xk+1 = xk +∆x
(5.62)
λk+1 = λk +∆λ,
where we have merged the solutions for the slack variables to xk+1. That is xk+1 = [xk+1,sk+1].
We choose the initial Lagrangian multiplier vector λ0 to follow the first order optimality
conditions or Karush-Kuhn-Tucker (KKT) conditions (see Eq.(2.54)):
λ0 = [Jg(x0,s )J′0 g(x ,s )]
−1
0 0 (J ′g(x0,s0)Jz(x0,s0)), (5.63)
and the initial slack variable vector to be a set of ones. More importantly, we define our
merit function to be the augmented function
γ
φ(x,γ) = z(x)+g(x)λ′+ ∥g(x)∥22 (5.64)2
(where γ is a positive real number) to track the convergence of the objective function. Note
that as the configuration of z changes, the merit function also changes too. The algorithm
follows a Newton optimization idea.
113
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
Algorithm 8: Capacity allocation algorithm for an arbitrary path network in min-
sum criterion
Result: Output the optimal cost and vector of edge capacities for the minsum
problem
Inputs: W, L, c, τ;
Allocate equal capacities to determine the initial z configuration;
Construct the associated Lagrangian function L(x,s,λ) ;
Compute the Jacobian matrix JL of L ;
Compute the right hand side vector ∇x,λL ;
Initialize a monotone initial condition;
while ∇L(y)> ε do
Compute ∆x =−J−1L ×∇L ;
Compute the next solution yi+1 = yi +∆x ;
Check if z(yi) has changed the configuration and compute the new z function ;
end
Stop iteration when tolerance is met.
From the theory of quadratic programming problem (see [33]), we understand that when
the number of constraints m is greater than or equal to 1, the KKT matrix in (5.60) is
indefinite. To circumvent this problem and have effective direct method, we employ a
symmetric indefinite factorization to decompose the KKT matrix. This factorization in this
context is given by
P′JLP = LdB L′d d, (5.65)
where P is a permutation matrix, Ld is a unit lower triangular matrix, and Bd is a block
diagonal matrix with either 1× 1 or 2× 2 blocks. This allows us to solve for [∆x ∆λ]′
systematically as follows: Solve
Ld∆y =−P′∇x,λL(x,λ) (5.66)
114
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
for unknown ∆y, solve
Bd∆̂y = ∆y (5.67)
for unknown ∆̂y, and solve
L′d∆y = ∆̂y (5.68)
for unknown ∆y. Finally, set  ∆x= P∆̂y. (5.69)
∆λ
Also for some specific problems, if B(x,λ) in (5.60) is positive definite and well condi-
tioned, then the Schur-Complement method can be employed (see [33] for details). How-
ever, B(x,λ) cannot be assumed to be positive definite in general.
The algorithm 8 begins by allocating equal capacities to construct initial z configuration
for the problem. It then goes further to construct the associated Lagrangian function. Due to
the possible congestion parameter η, the algorithm employs initial solution vector that has
no equal elements while ensuring a monotone allocation towards the sink. The algorithm
also incorporates a template to monitor the change in the configuration of z at each iteration.
5.4.2 The Hessian approximation
The basic Newton SQP may encounter convergence difficulty if the initial solution x0
is chosen too far from the true solution. Moreover, when x0 is far from x⋆, the sequence
of Hessian approximations of the Lagrangian may not be positive definite on the required
subspace which may hinder the existence of solutions of the quadratic subproblems. Based
on the above reasons, the class of Hessian approximations discussed in Subsection 2.5.6
and Subsection 2.5.7 are employed. This class of approximations is referred to as secant
approximation and has been found efficient for unconstrained non-linear problems. Basi-
cally, the Hessian approximations satisfy the secant equation, which is a version of bounded
deterioration property and have some attractive convergence properties. We describe the se-
115
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
cant approximation below. Recall from Taylor’s theorem that if the objective function L is
continuously differentiable and given p ∈ Rn, then we have (see for example [33])
L(x+p) = L(x)+∇L(x+ tp)′p, (5.70)
for some t ∈ (0,1). And if L is twice differentiable and continuous, we have
∫ 1
∇L(x+p) = ∇L(x)+ ∇2L(x+ tp)pdt. (5.71)
0
Now, adding and subtracting ∇2L(x)p in Eq.(5.71), we have
∫ 1
∇L(x+p) = ∇L(x)+∇2L(x)p+ [∇2L(x+p)−∇2L(x)]pdt. (5.72)
0
The above equation reduces to
∇L(x+p) = ∇L(x)+∇2L(x)p+O(∥p∥) (5.73)
due to the continuity of ∇L . Letting x = xk and p = xk+1 −xk, we have
∇L(x ) = ∇L(x )+∇2k+1 k L(xk)(xk+1 −xk)+O(∥xk+1 −xk∥), (5.74)
and with the assumption that the iterates are close to the solution, we have the secant equa-
tion
∇L(xk+1)−∇L(xk)≃ ∇2L(xk)(xk+1 −xk). (5.75)
Note that ∇2L = HL the Hessian of the Lagrangian. The secant approximations now find
Bk that satisfies
Bk+1(xk+1 −xk) = ∇L(xk+1,λk+1)−∇L(xk,λk+1). (5.76)
116
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
5.4.3 Numerical experiment
We show the results of our algorithm on 2-edge length, 3-edge length and 4-edge length
path networks. We discovered that the results satisfy the monotone allocation condition
even when the sink is located in the middle. The last rows of the tables present the optimal
costs and solution vectors.
w1 = 10 x1 w2 = 18 x2
v S1 l1 = 3 v2 l2 = 4
Figure 5.10: Path network with two edges.
z-value x1 x2
161.65 6.6667 13.333
161.62 7.4287 12.571
161.6 7.2103 12.79
161.6 7.144 12.856
161.6 7.1429 12.857
Table 5.1: Table showing the solution for 2-edged network with fixed budgeted capacity
c = 20 and 4 iterations.
z-value x1 x2
181.6 3.3333 6.6667
181.34 3.815 6.185
181.2 3.5622 6.4378
181.2 3.5714 6.4286
181.2 3.5714 6.4286
Table 5.2: Table showing the solution for 2-edged network with fixed budgeted capacity
c = 10. We reduce the budgeted capacity to let the initial z indicate a congestion scenario.
The optimal solution however indicate no congestion. The algorithm takes 4 iterations to
complete.
117
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
w1 = 10 x w2 = 18 x w3 = 151 2 x3
v S1 l1 = 3 v2 l2 = 4 v3 l3 = 3
Figure 5.11: Path network with three edges. Fixed capacity at c = 20.
z-value x1 x2 x3
321.55 3.3333 6.6667 10
317.23 4.7522 8.3469 6.9009
317.53 4.8909 7.5867 7.5224
317.52 4.6281 7.6863 7.6856
317.52 4.6365 7.6818 7.6818
317.52 4.6365 7.6818 7.6818
Table 5.3: Table showing the solution for 3-edge network with 5 iterations.
w1 = 10 x1 w2 = 18 x2 w3 = 15 x3 w4 = 10 x4
v1 l1 = 3 v2 l2 = 4 v3 l3 = 3 v4 l4 = 2
S
Figure 5.12: Path network with four edges.
118
5.4. SOLUTION BY SEQUENTIAL QUADRATIC PROGRAMMING
z-value x1 x2 x3 x4
472 2 4 6 8
448.36 3.4176 5.8985 6.4379 4.246
448.32 3.9184 6.4114 5.0478 4.6224
449.37 3.7532 5.7909 5.2318 5.2241
450.03 3.5245 5.6524 5.4125 5.4106
450.44 3.4614 5.5617 5.4885 5.4884
450.63 3.4264 5.5256 5.524 5.524
450.63 3.4255 5.5248 5.5248 5.5248
450.63 3.4255 5.5248 5.5248 5.5248
Table 5.4: Table showing the solution for 4-edged network with c = 20. The algorithm
takes 8 iterations to complete.
w1 = 10 x1 w2 = 18 x2 x3 w3 = 15 x4 w4 = 10
v1 l1 = 3 v2 l2 = 4 S l3 = 3 v3 l4 = 2 v4
Figure 5.13: Path network with four edges and a sink in the middle.
z-value x1 x2 x3 x4
308.68 3.3333 6.6667 6.6667 3.3333
308.31 3.7691 6.1067 6.4399 3.6843
308.13 3.6199 6.494 6.3946 3.4915
308.13 3.6163 6.5093 6.389 3.4853
308.13 3.6163 6.5094 6.389 3.4853
Table 5.5: Solution for the case where the sink is in the middle of the path. The solutions
are monotone from either sides. The algorithm takes 4 iterations to complete.
119
Chapter 6
Conclusion and future work
In this thesis we have developed and analysed novel algorithms and datastructure suitable
for the evacuation planning process of sparse or remote communities. We first of all focus
on the algorithm that allocates capacities optimally with the objective of minimizing the
minmax (or completion time) objective function. When the topology is a path network, we
are able to reduce the problem to finding an optimal allocation to the furthest edge network,
and with an efficient binary search query through our constructed parametric equation and
segment tree datastructure to compute the optimal capacities for other edges efficiently.
We extend the same idea to the case when the topology is a star network and the sink is
located at the centre of the network. However, when the sink is located at a leaf of the star
network, different approach proposed in [10] called the Klinz theorem is employed, only
for a star network with three edges. For an extended star network with a sink on the leaf,
the problem becomes too complex and require further investigation of hidden combinatorial
properties.
In Chapter 5, we proceed to developing algorithm that allocates capacities optimally to
the edges of path networks with the objective of minimizing the total evacuation time of all
evacuees on the networks. We notice that the problem is more difficult than the minmax
problem. We discover that the associated objective function is convex and monotone allo-
cation in the optimal solution vector applies. We also discover that the objective function
may change configuration from a congestion scenario to a non-congestion scenario during
the optimization process. We are able to solve the problem by incorporating the idea of
120
6. CONCLUSION AND FUTURE WORK
sequential quadratic programming (SQP), together with monotone intial solution vector.
The SQP allows us to incorporate the associated constraints and also helps in tracking the
changing configuration of the objective function.
The extension of the sink location algorithms studied in this thesis can be taken as
foundation for studying the problems on complex topologies such as trees. We are excited
to have published the developed algorithms for capacity provisioning problems for path
networks (minmax criterion) in the networks journals. We also added the idea discussed
for star networks with a sink in the middle in the same publication. We were able to submit
our results for the minsum problem to the EAMMO conference at the time of writing this
thesis. We plan to have the results published in a reputable journal.
One could pose a question to know how the capacity provisioning problems could help
emergency responders with logistics for an effective evacuation of a remote region. We
consider the following steps:
(i) identify the location of sites of interest and the number of evacuees at each site
(ii) compute a steiner tree spanning the sites. This tree becomes our network. Any exist-
ing roads can also be added to the network.
(iii) decide and fix location of evacuation sink or sinks based on various logistical needs
(iv) solve a capacity provisioning problem over the designed network.
The computed capacities can then inform emergency responders in making strategic deci-
sions concerning the allocation of emergency vehicles or aircraft. We envision that solving
the problem on path and star networks seem relevant for practice, and we believe that the
extension to tree topologies will offer more in reality.
We propose to embark on the following work in future:
(i) Multiple sinks capacity provisioning problems: We have an idea that when there are
multiple sinks on path networks, then there exist some optimal split edges which need
to be discovered for efficient algorithms.
121
6. CONCLUSION AND FUTURE WORK
(ii) Minmax capacity provisioning on extended star network with the sink on one of its
leaf nodes: The application of the Klinz theorem gives us some ideas as we have
noticed for the case of three-edged star network. However, the influence of multiple
flow functions encountered in the extended star topology need a further investigation
(iii) Capacity provisioning on tree network in minmax criterion: We hope to extend our
study of the minmax capacity provisioning to tree topologies. We understand that
the traditional minsum sink location problem is limited to path network in literature,
while the minmax problems have been extended to general graph through approxima-
tion algorithm techniques in [35]. We hope our ideas for star networks would show
some light to hidden properties in tree networks.
122
Bibliography
[1] Niels Henrik Abel. “Démonstration de l’impossibilité de la résolution algébrique
des équations générales qui passent le quatrieme degré”. In: Journal für die reine
und angewandte Mathematik 1 (1826), pp. 65–96.
[2] Jasbir Arora. Introduction to optimum design. Elsevier, 2004.
[3] Guru Prakash Arumugam, John Augustine, Mordecai J Golin, and Prashanth Srikan-
than. “A polynomial time algorithm for minimax-regret evacuation on a dynamic
path”. In: arXiv preprint arXiv:1404.5448 (2014).
[4] Robert Benkoczi and Rajib Das. “The minmax sink location problem on dynamic
cycle networks”. unpublished manuscript. August 2018.
[5] Binay Bhattacharya, Mordecai J Golin, Yuya Higashikawa, Tsunehiko Kameda, and
Naoki Katoh. “Improved algorithms for computing k-sink on dynamic flow path net-
works”. In: Workshop on Algorithms and Data Structures. Springer. 2017, pp. 133–
144.
[6] Binay Bhattacharya, Yuya Higashikawa, Tsunehiko Kameda, and Naoki Katoh. “An
O (n2 log2 n) Time Algorithm for Minmax Regret Minsum Sink on Path Networks”.
In: 29th International Symposium on Algorithms and Computation (ISAAC 2018).
Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik. 2018.
[7] Binay Bhattacharya and Tsunehiko Kameda. “Improved algorithms for computing
minmax regret sinks on dynamic path and tree networks”. In: Theoretical Computer
Science 607 (2015), pp. 411–425.
123
BIBLIOGRAPHY
[8] Paul T Boggs and Jon W Tolle. “Sequential quadratic programming”. In: Acta nu-
merica 4 (1995), pp. 1–51.
[9] Charles George Broyden, John E Dennis Jr, and Jorge J Moré. “On the local and su-
perlinear convergence of quasi-Newton methods”. In: IMA Journal of Applied Math-
ematics 12.3 (1973), pp. 223–245.
[10] Hoppe Bruce and Tardos Éva. “The quickest transshipment problem”. In: Mathemat-
ics of Operations Research 25.1 (2000), pp. 36–62.
[11] Richard H Byrd and Jorge Nocedal. “An analysis of reduced Hessian methods for
constrained optimization”. In: Mathematical Programming 49.1 (1990), pp. 285–
323.
[12] Di Chen and Mordecai Golin. “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.
[13] Di Chen and Mordecai J Golin. “Minmax Centered k-Partitioning of Trees and Ap-
plications to Sink Evacuation with Dynamic Confluent Flows”. In: arXiv preprint
arXiv:1803.09289 (2018).
[14] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein. Intro-
duction to algorithms. Chapter 8.2. MIT press, 2009.
[15] Daniel Dressler and Martin Strehler. “Capacitated confluent flows: complexity and
algorithms”. In: International Conference on Algorithms and Complexity. Springer.
2010, pp. 347–358.
[16] Roger Fletcher. “The sequential quadratic programming method”. In: Nonlinear op-
timization. Springer, 2010, pp. 165–214.
[17] J Flum and M Grohe. Parameterized Complexity Theory. 2006. Texts in Theoretical
Computer Science. An EATCS Series, 2006.
124
BIBLIOGRAPHY
[18] Greg N Frederickson and Donald B Johnson. “Finding kth paths and p-centers by
generating and searching good data structures”. In: Journal of Algorithms 4.1 (1983),
pp. 61–80.
[19] Philip E Gill and Elizabeth Wong. “Sequential quadratic programming methods”. In:
Mixed integer nonlinear programming. Springer, 2012, pp. 147–224.
[20] Mordecai J Golin, Hadi Khodabande, and Bo Qin. “Non-approximability and poly-
logarithmic approximations of the single-sink unsplittable and confluent dynamic
flow problems”. In: arXiv preprint arXiv:1709.10307 (2017).
[21] Shih-Ping Han. “Superlinearly convergent variable metric algorithms for general
nonlinear programming problems”. In: Mathematical Programming 11.1 (1976), pp. 263–
282.
[22] Shih-Ping Han. “A globally convergent method for nonlinear programming”. In:
Journal of optimization theory and applications 22.3 (1977), pp. 297–309.
[23] Yuya Higashikawa, Mordecai J Golin, and Naoki Katoh. “Minimax regret sink lo-
cation problem in dynamic tree networks with uniform capacity”. In: International
Workshop on Algorithms and Computation. Springer. 2014, pp. 125–137.
[24] Yuya Higashikawa, Mordecai J Golin, and Naoki Katoh. “Multiple sink location
problems in dynamic path networks”. In: Theoretical Computer Science 607 (2015),
pp. 2–15.
[25] Yuya Higashikawa and Naoki Katoh. “A Survey on Facility Location Problems in
Dynamic Flow Networks”. In: The Review of Socionetwork Strategies 13.2 (2019),
pp. 163–208.
[26] Bruce Hoppe and Éva Tardos. “Polynomial Time Algorithms for Some Evacuation
Problems.” In: SODA. Vol. 94. 1994, pp. 433–441.
125
BIBLIOGRAPHY
[27] Naoyuki Kamiyama, Naoki Katoh, and Atsushi Takizawa. “Theoretical and practical
issues of evacuation planning in urban areas”. In: The Eighth Hellenic European Re-
search on Computer Mathematics and its Applications Conference (HERCMA2007).
2007, pp. 49–50.
[28] Van Marc Kreveld, Otfried Schwarzkopf, Mark de Berg, and Mark Overmars. Com-
putational geometry algorithms and applications. Springer, 2000.
[29] Ford Jr R. Lester and Fulkerson Delbert Ray. “Constructing maximal dynamic flows
from static flows”. In: Operations research 6.3 (1958), pp. 419–433.
[30] Satoko Mamada, Takeaki Uno, Kazuhisa Makino, and Satoru Fujishige. “An O (n log2 n)
algorithm for a sink location problem in dynamic tree networks”. In: Exploring New
Frontiers of Theoretical Informatics. Springer, 2004, pp. 251–264.
[31] Nimrod Megiddo. “Combinatorial optimization with rational objective functions”.
In: Math and Operation Research 4 (1979), pp. 414–424.
[32] Rolf Niedermeier. Invitation to fixed-parameter algorithms. Vol. 31. OUP Oxford,
2006.
[33] Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science &
Business Media, 2006.
[34] UM Palomares and Olvi L Mangasarian. “Superlinearly convergent quasi-Newton al-
gorithms for nonlinearly constrained optimization problems”. In: Mathematical Pro-
gramming 11.1 (1976), pp. 1–13.
[35] Belmonte Rémy, Higashikawa Yuya, Katoh Naoki, and Okamoto Yoshio. “Polynomial-
time approximability of the k-Sink Location problem”. In: arXiv preprint arXiv:1503.02835
(2015).
[36] Benkoczi Robert, Bhattacharya Binay, Higashikawa Yuya, Kameda Tsunehiko, and
Katoh Naoki. “Minsum k-sink problem on dynamic flow path networks”. In: Inter-
national Workshop on Combinatorial Algorithms. Springer. 2018, pp. 78–89.
126
6. CONCLUSION AND FUTURE WORK
[37] Benkoczi Robert, Bhattacharya Binay, Higashikawa Yuya, Kameda Tsunehiko, and
Katoh Naoki. “Minsum k-sink problem on path networks”. In: Theoretical Computer
Science (2019).
[38] Mamada Satoko, Uno Takeaki, Makino Kazuhisa, and Fujishige Satoru. “An O(n log2 n)
algorithm for the optimal sink location problem in dynamic tree networks”. In: Dis-
crete Applied Mathematics 154.16 (2006), pp. 2387–2401.
[39] F Bruce Shepherd and Adrian Vetta. “The inapproximability of maximum single-sink
unsplittable, priority and confluent flow problems”. In: arXiv preprint arXiv:1504.00627
(2015).
[40] Richard A Tapia. “Diagonalized multiplier methods and quasi-Newton methods for
constrained optimization”. In: Journal of Optimization Theory and Applications 22.2
(1977), pp. 135–194.
[41] Robert B Wilson. “A simplicial algorithm for concave programming”. In: Ph. D.
Dissertation, Graduate School of Bussiness Administration (1963).
127