Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Reconfiguration strategies for mitigating the impact of port disruptions
(USC Thesis Other)
Reconfiguration strategies for mitigating the impact of port disruptions
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
RECONFIGURATION STRATEGIES FOR MITIGATING THE IMPACT OF PORT
DISRUPTIONS
by
Hwan Chang
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
August 2009
Copyright 2009 Hwan Chang
ii
Table of Contents
List of Tables iv
List of Figures v
Abstract vi
Chapter 1 Introduction 1
Chapter 2 Terminal level 6
2.1. Mitigating disruptions at terminal level 6
2.2. Formulations for Berth Allocation Problem 11
2.2.1. Relative Positioning Formulation for the BAP 12
2.2.2. Assignment formulation for the BAP 14
2.2.3. Mixed Integer Linear Programming Formulation for the BAP 15
2.2.4. Alternative choice for the objective function 16
2.3. Solution methods for the BAP 18
2.3.1. Lagrangian relaxation method 19
2.3.2. Simulated annealing method for the BAP 25
2.3.3. Heuristic method for an initial feasible solution 29
2.3.4. Heuristic method for an improved feasible solution 31
2.4. Computational Experiments for the BAP 35
2.4.1. CPLEX MIP: Computational experiments 35
2.4.2. Heuristic methods: Computational experiments 38
Chapter 3 Port Level 40
3.1. Mitigating disruption at port level 40
3.2. Multiple berth allocation problem 40
3.3. Solution method for the mBAP 43
3.4. Consideration of partially functional berth (terminal) 47
3.5. Computational experiments 51
3.5.1. Vessels with flexible berthing locations and departure times 52
3.5.2. Vessels with strict berthing locations and departure times 54
Chapter 4 Regional Level 57
4.1. Service network design 57
4.2. Service network for the U.S. west cost region 59
4.3. Mitigating disruptions at regional level 61
4.4. Computational experiments 65
4.4.1. Distributions of freights from LA port node (PART I) 69
4.4.2. Distributions of freights from all port nodes (PART II) 71
Chapter 5 Network link level 73
5.1. Mitigating disruptions at link level 73
5.2. The macroscopic traffic flow model 74
5.2.1. Variables and model of the METANET 74
iii
5.2.2. Extensions of the METANET model for speed limits 77
5.3. Model predictive control 78
5.4. ALINEA-base ramp metering and speed limit control strategy 79
5.5. Simulation experiments 81
5.5.1. Macroscopic simulation 81
5.5.2. Microscopic simulation 84
5.6. Integrated simulations for future works 86
Conclusions 89
References 91
iv
List of Tables
Table 2.1: Ranges of parameters 36
Table 2.2: Comparison between the BAP solutions using IP, SG, and SA 37
Table 2.3: Performance comparison between SG and SA methods 38
Table 3.1: Scenarios 1 and 2, experimental results 53
Table 3.2: Scenarios 3 and 4, experimental results 54
Table 3.3: Scenarios 5 and 6, experimental results 55
Table 3.4: Scenarios 7 and 8, experimental results 56
Table 4.1: Import Freight Distribution (ton/day) via Major West Coast Ports 66
Table 4.2: Legends and numbers of nodes describing FAF zones 66
Table 4.3: Activation and transportation cost of sea links 68
Table 4.4: Estimating capacity (ton) of port links 69
Table 4.5: Estimating capacity of ground links 69
Table 4.6: Freight distribution from LA port by activating sea links 70
Table 4.7: Freight distribution from all ports by activating sea links 71
Table 5.1: Performances under macroscopic environments 83
Table 5.2: Performance of the proposed DERAM strategy 86
v
List of Figures
Figure 2.1: Space-time diagram for the berth allocation problem.............................................. 7
Figure 2.2: Changes in the original allocation plan.................................................................... 9
Figure 2.3: Defining a vessel in the space-time diagram ......................................................... 12
Figure 2.4: Vessel-berth allocation in grid spaces.................................................................... 14
Figure 2.5: Possible optimal solutions under different cost functions ..................................... 18
Figure 2.6: Initial and final solutions by subgradient method.................................................. 24
Figure 2.7: Normalized objective values.................................................................................. 28
Figure 2.8: Average number of iterations and the number of minima ..................................... 29
Figure 2.9: Reduction by SWAPx resulting from Swap-Pair followed by Push-Down........... 33
Figure 2.10: Reduction by RIGHT resulting from Push-Right followed by Push-Down ........ 33
Figure 2.11: Reduction by LEFT resulting from Push-Left followed by Push-Down ............. 34
Figure 3.1: Space-time diagram for the Terminal Allocation Problem.................................... 41
Figure 3.2: Value of the block.................................................................................................. 45
Figure 3.3: Examples of partially functional berths ................................................................. 48
Figure 4.1: Service network for the U.S. west cost region....................................................... 59
Figure 4.2: Reconfiguration of the regional service network................................................... 62
Figure 5.1: Uni-directional freeway stretch divided into N sections........................................ 75
Figure 5.2: Fundamental flow-density diagram ....................................................................... 81
Figure 5.3: Architecture for an integrated simulation frameworks .......................................... 88
vi
Abstract
Currently, due to the continuing increases in demand, the U.S. freight transportation
system is running at capacity and disruptions are commonplace. However, by improving the
operational characteristics of the major resources in the network, the underutilization of
capacities can be reduced and, therefore, the impact of disruptions can be mitigated. The
purpose of this research is to study methods of modeling and evaluating the freight operations
and distributions and to develop mitigation strategies for reducing their impact at terminal,
port, and regional levels with special emphasis given to the ports of the U.S. west coast region.
The objective at the terminal level is to minimize the effect of disruption on the terminal
throughput. In order to effectively utilize the capacity of terminal, we focus on the scheduling
problem of calling vessels on a berth. It is the most essential planning since terminal operators
are required to coordinate the deployment of various resources so that freight are moved as
smoothly and as quickly as possible. We developed mathematical models and optimization
methods for the berth allocation problem (BAP). The objective at the port level is to route
ships to different terminals within the port so that the overall port throughput is affected less
by disruption. We focus on the scheduling problem of calling vessels on multiple berths which
is referred to as the multiple berth allocation problem (mBAP). We developed mathematical
models and optimization methods for the mBAP. The objective at the regional level is to
develop a mitigation strategy so that the regional throughput in moving goods is less affected
by disruptions. The regional service network is defined at a high level of aggregation, which
includes the major ports and aggregated zones representing broad geographical destinations
and intermediary locations. The regional service network under disruption is modeled as the
minimum cost flow problem with binary (or, integer) constraints. Furthermore, toward future
integrated evaluation works, we address the problem of improving the capacity utilization of a
vii
transportation link (a freeway stretch) using an integrated freeway traffic control system
comprised of ramp metering and variable speed limit strategies.
1
Chapter 1 Introduction
The elimination of international trade barriers, lower tariffs and shifting centers of global
manufacturing and consumption has led to new dynamics in intermodal shipping. Worldwide
container trade is growing at 9.5% annually, and the U.S. growth rate is around 6%. The
growth in containerized trade will continue as more and more cargo is transferred from break-
bulk to containers [46]. In today’s global economy, the oceans are becoming increasingly
important for international trade, since more than 80% of the world’s trade travels by water. A
very important component of this global economic chain is container transport, since about
half the world’s trade by value, and 90% of the general cargo, are transported in containers.
Shipping is the heart of the global economy, but it is vulnerable to disruptions. Trade
passes primarily through a small number of hubs spread around the globe. Close to 75% of the
world’s maritime trade and half of its daily oil consumption pass through a handful of
international straits and canals. Hence, the international commerce is at great risk from
disruptions at one of the major trading hubs or at one of a handful of strategic chokepoints
[12]. The adoption of a just-in-time delivery approach to shipping by most industries, rather
than stockpiling or maintaining operating reserves, means that a disruption or slowing of the
flow of almost any item can have widespread implications for the overall market, as well as
upon the national economy.
Disruptions to the maritime transportation system could be due to natural causes, (such as
hurricanes or earthquakes), or to man-caused activities (such as military surge or terrorism
acts). Moreover, the disruptions can be classified as predictable/anticipated (such as the
longshoremen strike at the port of Los Angles in 2002), or unpredictable/unanticipated (such
as a potential terrorist attack). Also, the point where the disruption occurs is a very critical
parameter in determining the disruption’s impact. For the purposes of our research, we will
2
concentrate on port disruptions with special emphasis given to the ports of the U.S. west coast,
because of their economic significance.
If we focus our attention on the four main container port complexes on the West Coast
(Long Beach/Los Angeles, San Francisco/Oakland, Portland, and Seattle/Tacoma), we see that
they handle almost 50% of the container traffic in the United States. In particular, the
combined ports of Long Beach/Los Angeles, the first and second larges container ports in the
US, handled roughly one-third of all import and export containers in the nation [1] This huge
volume moving through the local ports has very serious effects not only at the local and
regional levels, but on a national scale as well. In 1984, only one weekly eastbound double-
stack rail service existed from the whole West Coast (from POLA/POLB to Chicago). Less
than ten years later, in 1993 there were 241 train sets of weekly eastbound double-stack rail
services from all major ports in the West Coast to several metropolitan centers in the east [48]
A consequence of this growth is the fact that the national economy has become heavily
dependent on the smooth and reliable operation of the west coast ports.
Apparently, due to the continuing increases in demand, the US freight transportation
system is running at capacity and disruptions are commonplace. The situation can be
translated into various impact including unexpectedly increasing shipping costs and transit
times. However, considering the geographical location of the west cost ports, the scarcity of
lands and the limitation of budgets force the port and terminal operators to increase capacity
by applying more intelligent operational plans rather than by expanding its physical resources,
such as yards and equipments.
The purpose of this research is to study methods of modeling and evaluating the freight
operations and distributions and to develop mitigation strategies for reducing their impact at
the terminal, port, and regional levels.
3
Terminal level
The objective at the terminal level is to minimize the effect of disruption on the terminal
throughput.
In a container terminal, allocating berth spaces to calling vessels is the most essential for
efficient terminal utilization. The reason is that, based on a berth allocation planning, terminal
operators are required to coordinates the deployments of various resources within the entire
terminal so that containers are moved as smoothly and as quickly as possible.
In order to effectively utilize the capacity of terminal, we focus on the scheduling
problem of calling vessels on a berth which is referred to as the berth allocation problem
(BAP). We developed mathematical models and optimization methods for BAP. We consider
the continuous BAP as it is more diverse and practical. Since the continuous BAP cannot be
calculated in polynomially-bounded time, we develop and implement heuristic procedures
based on subgradient and simulated annealing optimization methods. For both optimization
methods, a set of systematic and efficient heuristics are implemented to find an initial feasible
solution and to update a current solution by exploring a sufficiently large solution space.
Port level
The objective at the port level is to route goods to different terminals within the port so
that the overall port throughput is affected less by disruption.
In order to effectively mitigate the impact of disruptions and to maintain terminal
operations and productivity, we focus on the scheduling problem of calling vessels on multiple
berths which is referred to as the multiple berth allocation problem (mBAP). We developed
mathematical models and optimization methods for the mBAP.
4
A heuristic set partitioning procedure is developed to generate a partition between
different berths (terminals) for solving the mBAP. Each sub-problem is repeatedly solved by
the previously developed simulation annealing method for BAP. Experimental scenarios are
prepared to evaluate the proposed procedures in the presence of a partially functional berth
(terminal) as well as some realistic requests issued by vessels.
Regional level
The objective at the regional level is to mitigation strategies so that the regional
throughput in moving goods is less affected by disruption.
We look into the US west coast region, consisting of multiple ports and the associated
traffic network used for moving the goods within and out of the region. The regional service
network is defined at a high level of aggregation, which includes the major ports and
aggregated zones representing broad geographical destinations and intermediary zones.
The regional service network under normal case is modeled as the minimum cost flow
problem in which the disruption level is low enough to be handled by ground transportation
modes. In a case of a disruption, the regional service network can be reconfigured. For
example, the LA port zone is rendered non-functional for a period of time, all services
emanating the zone either be discontinued or operating at a lower capacity, which will affect
its capacity characteristics. The reconfiguration of the service network is performed by
opening sea transportation links between port zones. The regional service network under
disruption is modeled as the minimum cost flow problem with binary (or, integer) constraints.
It is solved by branch-and-bound method and a LP relaxation is performed every leaf of the
branch-and-bound tree.
5
Network link level for future integration
Toward a future integrated work with the above described mitigation strategies, we briefly
address the problem of improving the operational characteristics of a link (a freeway stretch)
which is the most vulnerable link in the regional transportation network.
The increasing demand accompanied by the increasing throughput in maritime shipping
can also be translated into the increase in the port generated traffic. This traffic has emerged as
a major contributor to congestion on regional freeway stretches which are the most vulnerable
transportation links to capacity constraints and disruptions.
In order to effectively reduce congestion and improve safety on a freeway link, we
consider an integrated freeway traffic control system comprised of ramp metering and variable
speed limit strategies. A computationally feasible and practical controller is developed for the
dynamic traffic system to control flow rates of a ramp and maximum speeds of a mainline
section. The performance of the control system is evaluated both using a conventional second-
order macroscopic traffic flow simulation model and one of microscopic traffic flow
simulation models which have been recently used to support ITS applications.
6
Chapter 2 Terminal level
2.1. Mitigating disruptions at terminal level
A container terminal is a facility where containerized cargo is trans-shipped between ships
and land vehicles (trucks and trains). A terminal may have several wharfs (quays). Each quay
may consist of several berths, which in turn are divided into sections. Each berth corresponds
to a linear stretch of space in the terminal.
Container ships are moored at a berth of a terminal where they are unloaded and loaded
by gantry or quay cranes. Quay cranes traverse along the quay to position containers at any
point along the length of the ship. Quay cranes at berths load imported containers on in-yard
trucks, straddle carriers, or automated guided vehicles (AGVs). Quay cranes also unload
export containers from these vehicles into the ships. To maximize berth utilization and
minimize ship-turn-around time, ships should be optimally assigned and allocated to berths.
Hence, optimal allocation of berths to incoming ships will have a substantial impact on
logistics cost and level of service.
Usually, carriers inform the terminal operator of their estimated-time-of-arrival (ETA),
latest possible service completion (departure) time, and request a berthing time (called berth-
time-requested or BTR) several days in advance [42]. vessel is said to be berthed-on-arrival
(BOA) if the mooring operation commences within 2 hours of arrival. The BOA statistics is
often used as an indicator to gauge the quality of service provided by the port operator [39].
Based on the information received, a terminal operator tries to satisfy the requested
departure times of every vessel by allocating one or more sections on a berth to calling vessels
according to their ETAs, estimated departure times, and BTRs. However, in cases when the
arrival rate of vessels is high, or unexpected arrivals occur, or any type of disruption happens
7
at a terminal, it may not be possible to serve all the vessels before their requested service
completion time. Thus, departures of some vessels may be delayed past the requested due time.
Therefore, in order to effectively utilize the capacity of terminal, we focus on the
scheduling of calling vessels on a berth which is referred to as the berth allocation problem
(BAP).
Figure 2.1: Space-time diagram for the berth allocation problem
Figure 2.1 illustrates the space-time diagram of a berth allocation schedule. The
horizontal axis represents the berth length, while the planning horizon is represented by the
vertical axis. Note that an allocated vessel is represented by a rectangle in the diagram. The
length of the vertical side of each rectangle represents the duration of stay of a vessel at the
berth, while the length of the horizontal side represents the vessel length. Each calling vessel
is characterized by its own space-time rectangle. These rectangles cannot overlap either in the
space or in the time dimension. The Berth Allocation Problem (BAP) is to determine the
optimal locations of those rectangles without overlaps. In other words, the BAP consists of
8
optimally allocating and scheduling the berth space to calling carriers such that the carriers are
served within their time limits.
The BAP can be classified into the following two general categories.
Discrete vs. Continuous BAP
In the continuous BAP, the berthing can be done in a continuum of locations along the
berth [29],[30]. In contrast, in the discrete BAP, the entire quay is divided into a countable
number of berths. In the discrete problem, it is assumed that, at each instant of time, at most
one vessel can be served at each berth. It should be noted that a discrete BAP can be modeled
as an unrelated parallel machine scheduling problem, and the continuous BAP can be mapped
into a two-dimensional cutting-stock problem, which is an NP-Hard problem [11]. The focus
of this study is on the continuous BAP. The continuous BAP is more diverse and general. It
can address the growing trends in ship sizes, as there is a need for more flexible berth
allocation planning. For instance, mega-ships may sometimes moor across neighboring berths
in order to enhance berth usage.
Static vs. Dynamic BAP
In Imai et al. [27], the static BAP refers to the problem in which all the ships are assumed
to have arrived at the port prior to the beginning of the current scheduling. In contrast, the
dynamic BAP takes into account not only the ships that have already arrived at the time of
planning, but also those which will arrive later during the planning horizon. Furthermore, it is
assumed that the arrival times of all the ships are known a priori, hence re-planning is not an
issue. In a general planning problem, however, the dynamic environment refers to the events
whose occurrences may change in time. For example, in the work by Moorthy and Teo [39],
the authors deliberately induce delays in the port stay time of vessels and increase the number
of vessels to evaluate their dynamic policy in a rolling horizon framework. That is to say that
9
the notion of dynamic BAP used in [27] is very limited and, therefore, we will use the second
notion in which the occurrence of events may change in time.
Berth
Time
Berth
Time
Berth
Time
(a) (c) (b)
Figure 2.2: Changes in the original allocation plan
Figure 2.2 demonstrates possible disruptions in an original berth allocation. Figure 2.2a
shows the original assignment of ship. Disruptions may have occurred due to delays in vessel
arrival time (Figure 2.2b), or vessel service time (Figure 2.2c). If these unexpected delays
make the original allocation plan unsatisfactory, re-allocation of vessels may be needed to
maintain the same level of quality of berth performance. The re-allocation can easily be
implemented by using a rolling horizon framework in such a dynamic environment. It should
be noted that frequent re-planning is often undesirable and sometimes impossible as it has
adverse impact on other terminal resources. So, the re-planning of a berth allocation should be
carefully performed by considering the impact on all in-terminal resources.
Recently, berth allocation problems have been the focus of many research efforts. Imai et
al. proposed a mixed integer programming formulation of the discrete berth allocation
problem [27]. Two formulations are developed for static and dynamic variants of the problem.
A Lagrangian relaxation methodology equipped as a heuristic method was developed to
minimize the ships’ waiting and handling times.
Imai et al. in [28] developed a formulation and a solution methodology for the discrete
berth allocation problem with priority considerations. In this work, they extended their
10
previous work in [27] to serve calling vessels at various service priorities. A heuristic method
based on genetic algorithms was developed to approximately solve the problem with less
computational burden. In [29], Imai et al. considered a continuous berth allocation problem in
which a vessel can be moored across the designated quay boundaries. The authors developed a
heuristic solution based on their discrete BAP solution in [27]. They used a series of local
procedures to ensure the feasibility of the solution. The approach was based on the fact that an
optimal solution to the BAP provides an upper bound, when the berth length is set to the
maximum ship length, and provides a lower bound, when the minimum ship length is used.
Guan and Cheung developed a composite heuristic for a BAP whose objective is to minimize
the weighted completion time [22]. A batch was defined as a group of ships whose total size is
smaller than the overall berth space. As an exact solution method, a tree search procedure was
proposed to solve small-sized problems. In the composite heuristic proposed, a pair-wise
exchange was performed between batches. The tree search procedure was applied to enhance
the solution of each batch. Park and Kim in [42] addressed a BAP with a general objective that
minimizes the costs resulted from the vessels delayed departure times. The objective function
also consists of additional handling costs, which were resulted from deviated berthing
locations. The BAP formulation is very practical since the cost related to the delayed departure
times is explicitly minimized. In practice, penalties are imposed if the requested service cannot
be done by the requested departure time. A subgradient optimization technique was applied to
solve the proposed BAP formulation. A heuristic method was repeatedly used to update an
upper bound. The method attempted to move overlapping schedules to feasible directions
which yielded the minimum cost increase. Kim and Moon in [34] proposed a simulated
annealing based method to solve the BAP. Their heuristic resembled the tree search proposed
in [22] since it sequentially tried to locate each ship at lowest-cost point. Moorthy and Teo in
11
[39] studied the allocation of preferred berthing space to a set of vessels which arrive
periodically in a weekly basis. The authors used the concept of sequence pair for defining
search space. By defining the time and space constraint separately, cost estimation for each
dimension was provided. Several neighborhood searches were employed by a simulated
annealing method to modify or update a sequence pair. However, due to the periodicity of the
weekly arrivals, the authors focused on reducing the overlaps between unfinished vessels
rectangles in the current planning horizon and scheduled rectangles in the next horizon.
Frequently, the scheduling of quay cranes is simultaneously considered in a berth
scheduling [21],[43]. The problem can also be formulated as an integer programming problem
with objective of minimizing the costs resulting from the delays. In [21], authors used the
genetic algorithms to solve the problem. And, in [43], authors applied a subgradient
optimization technique for the first phase of the berth allocation and applied the dynamic
programming technique for the quay crane schedule based on a near-optimal solution of the
first phase.
2.2. Formulations for Berth Allocation Problem
In this section, we study the static continuous berth allocation problem which can be
represented by a space-time diagram where the horizontal and vertical axes represent the
berthing space and time, respectively. Figure 2.3 shows the representation of a berth in the
space-time diagram. A vessel k with length
k
l and width
k
h can be defined by a rectangle
in this diagram.
12
k
d
k
t
kk
th +
k
a
k
x
kk
x l +
vessel k
Berth
Time
Figure 2.3: Defining a vessel in the space-time diagram
The mathematical formulation for the BAP can be either represented by the relative
position of vessel rectangles, or by the space covered by the vessel rectangles.
2.2.1. Relative Positioning Formulation for the BAP
The relative positioning formulation, which hereafter is called formulation [PR],
considers the relative positioning of vessel rectangles. The objective function (1) minimizes
the penalty incurred by not satisfying the desired departure times requested by vessels.
[PR] min
1
()
K
kk k k
k
wt h d
+
=
+−
∑
(1)
s.t. [( ) ] 0
lk k kl
xx l α σ −+ − ⋅ ≥ ,, kl k l ∀ ≠ (2)
[( ) ] 0
lk k kl
tt h τ δ −+ − ⋅ ≥ ,, kl k l ∀ ≠ (3)
1
kl lk kl lk
σ σδ δ ++ + ≥ ,, kl k l ∀ ≠ (4)
1
kl lk
σ σ +≤ ,, kl k l ∀ ≠ (5)
1
kl lk
δ δ +≤ ,, kl k l ∀ ≠ (6)
[1, 1]
kk
xMl ∈−+ k ∀ (7)
13
[, 1]
kk k
ta T h ∈−+ k ∀ (8)
{0,1 }
kl
σ ∈ ,, kl k l ∀ ≠ (9)
{0,1 }
kl
δ ∈ ,, kl k l ∀ ≠ (10)
where
M Number of sections along the berth (length of the berth)
T Time horizon
K Number of vessels to be scheduled
k
a Estimated arrival time of ship (vessel) k
k
d Desired departure time of ship k
k
h Estimated handling time needed for ship k
k
w Priority factor (weight) assigned to ship k
k
l The length of ship k in terms of the number of sections along the berth
k
x The section of the berth for which the left bottom corner of ship k is assigned
k
t The time at which berthing of ship k starts
kl
δ The relative horizontal position of ship k with respect to ship l . It is 1 if ship k is
completely to the left of ship l , and 0 otherwise
kl
σ The relative vertical position of ship k with respect to ship l. It is 1 if ship k is
completely placed below ship l , and 0 otherwise
α The minimum space required between two ships moored at the same time at the berth (i.e.,
safety allowance)
τ The minimum time needed between a ship departure time and the next ship berthing time
at the same location (i.e., time stability factor)
14
The notation ()
+
⋅ represent the max operation, i.e., max{0, } yy
+
= . Constraints (2) are
to ensure the space requirement if a vessel is completely placed to the left of another vessel.
Constraints (3) are to guarantee the time requirement if a vessel is completely placed below
the other vessel. Constraints (4)-(6) ensure that no vessel rectangles are overlapping.
Constraints (7) are the space constraints, (8) are time constraints, and (9), (10) are binary
constraints.
2.2.2. Assignment formulation for the BAP
The assignment formulation, which hereafter is called formulation [P], is based on
partitioning the space-time diagram into grid tiles which are covered by vessel rectangles as
shown in Figure 2.4. A block (, ) ij is referred to as a grid space whose left bottom corner is
located at (, ) ij for [1 , ] iM ∈ and [1 , ] jT ∈ . As an example, the vessel rectangle k covers
the consecutive blocks from (, ) ij to (1, 1)
kk
il j h + −+ − .
k
d
k
a
rectangle k
(, )
kk
x t
i
j
k
l
M
T
k
h
1 2
2
k
il +
Figure 2.4: Vessel-berth allocation in grid spaces
Similar to (1), objective function (11) minimizes the penalty incurred by not satisfying the
desired departure times requested by vessels. However, in contrast to (1), the minimization in
15
(11) is not performed on a continuum but a discrete number of time periods (tiles) representing
the unmet requested departure time.
[P] min
11
11
()
kk
k
Ml Th K
ijk k k k
ki ja
zw j h d
−+ − +
+
== =
⋅+ −
∑∑ ∑
(11)
s.t.
1
1
K
ijk
k
y
=
≤
∑
, ij ∀ (12)
11
0
kk
il j h
ijk k k mnk
mi n j
zlh y
+− + −
==
−≤
∑∑
1,, 1, ,, 1,
kk k
iMl ja Th k ∀=−+= −+ (13)
11
1
1
kk
k
Ml Th
ijk
ija
z
−+ − +
==
=
∑∑
k ∀ (14)
{0,1 }
ijk
y ∈ ,, ijk ∀ (15)
{0,1 }
ijk
z ∈ 1, , 1, 1, , 1,
kk
iMl j Th k ∀=−+= −+ (16)
where
ijk
z : 1, if the left bottom corner of rectangle k is located at ( , ) ij ; 0, otherwise
ijk
y : 1, if rectangle k is covering block ( , ) ij ; 0, otherwise
Constraints (12) imply that a block is covered by at most one vessel rectangle. Constraints
(13) ensure that blocks covered by a vessel rectangle must be consecutive. Constraints (14)
imply that each vessel has only one berthing coordinate. Constraints (15) and (16) are binary
constraints.
2.2.3. Mixed Integer Linear Programming Formulation for the BAP
The following formulation [PL] represents the mixed integer linear programming (MILP)
equivalent to the previous BAP formulation, which will be used by CPLEX MIP solver for
comparison purposes.
16
[PL] min
1
K
kk
k
w β
+
=
∑
(17)
s.t.
kk k k k
th d β β
+−
+− = − k ∀ (18)
() ( 1)0
lk k kl
xx l N α σ −+ − − − ≥ ,, kl k l ∀ ≠ (19)
() ( 1)0
lk k kl
tt h N τ δ −+ − − − ≥ ,, kl k l ∀ ≠ (20)
,0
kk
ββ
+−
≥ k ∀ (21)
and Constraints (4)-(10).
The objective function in (17) minimizes the penalty cost resulting from the delay in the
departure time of vessels. A method for dealing with ( )
kk k
th d
+
+− variables is to introduce
new variables
k
β
+
and
k
β
−
, constrained to be nonnegative, and let
kk k k k
th d β β
+−
+− = − . It is
intended to have
kk k k
th d β
+
+− = or
kk k k
th d β
−
+ −=− , depending on whether
kk k
th d + −
is positive or negative. Then the original problem [PR] can be formulated as a MILP by adding
a corresponding constraint (18) to (21).
For an optimal solution to the problem, and for each value of k , we must either have
0
k
β
+
= or 0
k
β
−
= . That is true since otherwise we could have reduced both
k
β
+
and
k
β
−
by
the same amount while preserving feasibility. In other words, we could have reduced the cost,
which is in contradiction to optimality. Constraints (19) and (20) are to enforce the definitions
of
kl
σ and
kl
δ , and N is a large positive number.
2.2.4. Alternative choice for the objective function
The objective function used in (1) minimizes the penalty associated with the violation of
the desired departure times requested by the vessels. A generalized alternative for this
objective function is one that rewards early service operations. The following function is a
17
possible choice for generalized objective function, which is a linear combination of penalized
late departure time and rewarded early service operation.
{}
12
2
1
(( )) (( ))
K
kk k k k k k k
k
wt h d w d t h
+ +
=
=+−− −+
∑
J (22)
where
1
0
k
w > is the penalty weight, and
2
0
k
w > is the reward weight. We require that
12
kk
ww ≥ .
The setback is that
2
J may assume both negative and positive values. To make
2
J
non-negative, the possible maximum value of the second term is added to
2
J as an offset.
Hence, a better choice for the generalized objective function would be
()
12 2
3
1
(( )) (( )) (( ))
K
k k k k k k kk k k kk
k
w t h d wd t h wd a h
++
=
=+−− −+ + −+
∑
J . (23)
The objective function
3
J can also be rewritten as
3
2
12
1
(), for ( ) 0
(( ) ) ( ( )), for ( ) 0
K
kk k k k k
k
k k k k k k kk kk k
wt a t h d
wt h d w d a h t h d
=
⎧−+−<
=
⎨
+− + − + + − ≥
⎩
∑
J (24)
3
J tries to minimize the delayed departure times while maximizing the temporal margin
after completion of the service. A very interesting special case of
3
J - which has been
considered by other researchers (such as [39], [29]) - is the one in which
12
kk
ww = . This
special case will result in
4
1
1
()
K
kk k
k
wt a
=
=−
∑
J (25)
The objective function
4
J minimizes the delayed mooring time only. Furthermore, by
adding a fixed offset
1
kk
wh to
4
J , the choice of the objective function will be as follows.
18
()
5
11 1
11
() ( )
KK
kk k k k k k k k
kk
wt a wh w t h a
==
=−+ = +−
∑∑
J (26)
The objective function
5
J minimizes the delayed mooring and handling times. This
choice of objective function was considered by Guan and Cheung [22] and Imai et al. [27].
Figure 2.5 shows two different final solutions, when
4
J or
5
J is used. The two
solutions are both optimal if
1
k
w is the same for all the ships. However, if we consider a more
general situation in which a delayed departure disturbs the successive assignments, the first
solution cannot be optimal. In such a case, we can introduce the generalized objective function
and can penalize the delayed departure time with a higher penalty weight to get an optimal
solution.
2
a
2
d
1
a
1
d
1
a
1
d
2
a
2
d
Figure 2.5: Possible optimal solutions under different cost functions
2.3. Solution methods for the BAP
Since the continuous BAP is an NP-hard problem, we develop two approximation
methods to find good solutions in a reasonable amount of time. Without loss of generality,
formulation [P] is used to illustrate the solution methods. In the sequel, we will focus on this
formulation.
19
2.3.1. Lagrangian relaxation method
Many hard integer problems can be viewed as an easy problem complicated by a
relatively small set of side constraints. Dualizing side constraints produces a Lagrangian
problem which is easy to solve and whose optimal value is a lower bound for minimization
problems on the optimal value of the original problem. This method was termed “Lagrangian
relaxation” by Geoffrion [20], who developed a systematic methodology to construct the
lower bounds as a means of exploiting special problem structure.
The problem [P] can be converted into a relaxed problem using the Lagrange multipliers
0
ij
π ≥ , , ij ∀ , as follows:
()
LB
z π =
11
11 1 1 1 1
min () 1
kk
Ml T hKMTK
ijk k k k ij ijk
ki j i j k
zw j h d y π
−+ − +
+
== = = = =
⎛⎞
⎛⎞
⋅+ − + −
⎜⎟
⎜⎟
⎝⎠
⎝⎠
∑∑ ∑ ∑∑ ∑
[] LR
π
s.t. Constraints (13)-(16).
where
ij
π π ⎡⎤ =
⎣⎦
, 1, , iM = and 1, , jT = , is a matrix of Lagrange multipliers.
The dual problem of the relaxed problem [ ] LR
π
becomes
LB
z = max ( )
LB
z
π
π , [] D
where
11
11 1 1 1 1
()min () 1
kk
Ml ThKMTK
LB ijk k k k ij ijk
ki j i j k
zzwjhd y ππ
−+ − +
+
== = = = =
⎛⎞
=⋅+−+ −
⎜⎟
⎝⎠
∑∑ ∑ ∑∑ ∑
s.t. Constraints (13)-(16), and 0
ij
π ≥ .
If we ignore the last constant term of the minimization problem in [D], then it is
equivalent to
11
11 1 1 1 1
min ()
kk
Ml T hKMTK
ijk k k k ij ijk
ki j i j k
zw j h d y π
−+ − +
+
== = = = =
⋅+ − +
∑∑ ∑ ∑∑ ∑
(27)
20
When 1
ijk
z = , we have 1
ijk
y = for , , 1
k
mi i l = +− and , , 1
k
nj j h = +− . Therefore,
for a vessel k , the second term becomes
11
11
kk
il j h MT
ij ijk mn
ij mi nj
y π π
+− + −
== = =
=
∑∑ ∑ ∑
(28)
Then, this implies that the minimization problem can be considered as
11 1 1
11 1
min ()
kk k k
Ml Th i l j h K
ijk k k k mn
ki j mi nj
zw j h d π
−+ − + + − + −
+
== = = =
⎛⎞
+− +
⎜⎟
⎝⎠
∑∑ ∑ ∑ ∑
(29)
s.t. Constraints (14)-(16)
The solution to [ ] LR
π
can be obtained since (29) is separable in k . The optimal
coordinate for each sub-problem is calculated by setting 1
ijk
z = .
11
,
min ( )
kk
k
il j h
kk k mn
ij a
mi n j
wj h d π
+− + −
+
≥
==
+− +
∑∑
(30)
To optimize dual functions in Lagrangian relaxation, we use the subgradient method (SG)
for separable integer programming problems. All sub-problems are solved optimally to obtain
a subgradient direction. The subgradient method is an adaptation of gradient methods, in
which gradients are replaced by subgradients. For further discussion on subgradient methods,
see [25].
Given an initial value 0
ij
π = , a multiplier is generated by the following rule:
( ) { }
1
max0, 1
K
ij ij ijk
k
sy ππ
=
= +−
∑
(31)
where
ijk
y is from an optimal solution to [] LR
π
and s is a positive scalar step size. We
use the following step size which has been commonly adopted in practice [16].
( )
2
1
()
|| 1||
UB LB
K
ijk
k
zz
s
y
λ π
=
−
=
−
∑
(32)
21
where λ is a scalar satisfying 02 λ < ≤ and
UB
z is the best known feasible (upper-bound)
solution value obtained by applying a heuristic to formulation [P].
We use a general rule, which is to set 2 λ = for some fixed number of iterations. This
number is called maxIter hereafter. At each iteration, we successively halve both the values of
λ and maxIter until the value of maxIter reaches some threshold (here, number 4). Note that,
alternatively, the procedure can be stopped when λ reaches a threshold. To describe our
developed subgradient optimization procedure BAPsg, the following notation is used.
LB
z maximum lower bound from Lagrangian relaxation
UB
z
minimum upper bound
1 H
z initial upper bound found by Heuristic H1 (it is described in section 3.2.3)
2 H
z
updated upper bound found by Heuristic H2 (it is described in section 3.2.4)
minUB
best upper bound so far
maxLB
best lower bound so far
maxIter maximum number of iterations
BAPsg
1 Calculate the initial upper bound using H1,
1 UB H
zminUB z = =
Set 0
LB
zmaxLB == , 2 λ = , 0
ij
π = , and 2 maxIter K =
2 If 0.005 λ < , then stop. Otherwise, continue.
3 Update the lower bound
LB
z .
If
LB
zmaxLB > , then
LB
maxLB z = , 2 λ = , and 1 Iter = . Otherwise, 1 Iter Iter =+ .
4 Update the minimum upper bound using H2,
2 UB H
zz =
If
UB
zminUB < , then
UB
minUB z = .
22
If 0
LB
z = , then stop. Otherwise, continue.
5 If Iter maxIter > , then 0.5 λ λ = × , max{4, 0.5} maxIter maxITer = × , and 1 Iter = .
6 Update the step size s and the multiplier
ij
π . Go to step 2.
In Step 3 of the BAPsg procedure, the lower bound
LB
z is acquired by solving each sub-
problem (30). Then, the heuristics H1 and H2 are applied to find a feasible solution in Step 4.
Details about heuristics H1 and H2 will be provided in sections 3.2.3 and 3.2.4. In the end, the
procedure returns the best feasible solution (the minimum upper bound) and the maximum
lower bound. At the end of iteration of the procedure, an optimal solution to [D] may be still
infeasible for [P], which means that vessel rectangles may be overlapping. However, as the
multipliers corresponding to vessels with infeasible berthing are increased, the cost of [D] is
increased. This ensures that the solution associated to
LB
z gets close to a feasible solution of
the original problem [P].
Figure 2.6 demonstrates an instance of a BAP with initial and final solutions which were
obtained by the above procedure. In the figure, each batch is represented by a distinct color. A
batch is defined as a set of rectangles in the space-time domain, where the sum of their lengths,
including safe distances between them, is less than or equal to the berth length M. The batch is
a fundamental building block in our heuristic methods which will be described later. The left
figure is a berth allocation plan represented by color-filled rectangles. The figure on the right
shows the corresponding plan along with time constraints. The time constraints of each vessel
are represented with a surrounding dotted rectangle. The height of a dotted rectangle is
determined by the vessels’ arrival and desired departure times and the width is equal to the
vessel’s length. If a vessel rectangle is not fully enclosed by its (dotted) constraint rectangle, it
will results in a cost increase. As shown in Figure 2.6a, the initial feasible solution indicates
that two rectangles (A and B) are not completely enclosed by their constraint rectangles.
23
Figure 2.6b shows the final lower bound solution, which is still infeasible, since several
rectangles overlap. Finally, the final upper bound solution is shown in Figure 2.6c. In this
specific example, the procedure yields an optimal solution to [P] without penalty cost
associated with waiting times.
Various port disruption events can be graphically described by using a vessel rectangle
and its constraint rectangle. The disruption can result in changes to either the vessel rectangle
or to the constraint rectangle. For example: (a) If, because of equipment break-down, the
service time
k
h is delayed by ω hours, then the height of the vessel rectangle is lengthened
by ω hours. That is, the vessel rectangle could change. (b) If the arrival time
k
a is delayed
by ω hours, because of delays due to an unexpected event, the constraint rectangle should be
adjusted and, also, the berthing time should be adjusted so that
kk
ta ω ≥+ . That is, both the
vessel rectangle and its constraint rectangle could change. In both cases above, we can
graphically detect whether the changed rectangles (resulting from a disruptive event) disturb
an original allocation plan or not. If the original allocation plan is disrupted, and the
disturbance occurs in batch b, we will reallocate vessel rectangles in the later batches
(1, 2, bb ++ ).
24
200 400 600 800 1000 1200
20
40
60
80
100
120
140
160
berth [m]
time [h]
Initial feasible soultion
200 400 600 800 1000 1200
20
40
60
80
100
120
140
160
berth [m]
time [h]
Initial feasible soultion with time constraints
(a) Initial feasible solution
200 400 600 800 1000 1200
20
40
60
80
100
120
140
160
berth [m]
time [h]
Lower bound soultion
200 400 600 800 1000 1200
20
40
60
80
100
120
140
160
berth [m]
time [h]
Lower bound soultion with time constraints
(b) Lower bound solution (infeasible)
200 400 600 800 1000 1200
20
40
60
80
100
120
140
160
berth [m]
time [h]
Upper bound soultion
200 400 600 800 1000 1200
20
40
60
80
100
120
140
160
berth [m]
time [h]
Upper bound soultion with tim constaints
(c) Upper bound solution
Figure 2.6: Initial and final solutions by subgradient method
25
2.3.2. Simulated annealing method for the BAP
The second methodology proposed is based on the simulated annealing (SA) to find good
solutions to the BAP. The SA methodology was independently presented by Kirkpatrick et al.
[35] and Cerny [5]. The SA algorithm resembles the annealing process in metallurgy. In each
step of the algorithm, the current state (solution) is replaced, with some probability, by a
randomly generated neighboring state. This probability depends on the difference between the
current state energy (cost) and generated neighbor energy, as well as, on a global parameter T
(temperature) which is gradually decreased using a scheduled parameter r (cooling rate). The
SA algorithm is such that it makes the system ultimately move to lower energy states (called
downhill movements). Also, uphill movements prevent the algorithm from staying at local
minima.
Before describing the methodology in detail, we start by defining the variables and
parameters, neighbor search function, and acceptance probability function for our SA
procedure. The following notation is used to describe the variable and parameters in our SA
method.
sc current state
sn new state
sb best state
ec energy (cost) of current state, ec=E(sc)
en energy of new state, en=E(sn)
eb energy of best state, eb=E(sb)
0
T initial temperature
r cooling rate
26
The following neighbor (candidate move) search procedure, BAPneighbor, is defined for
out SA. A new state (solution) sn is obtained by exchanging a pair of consecutive vessel
rectangles in the current state sc . The current state (1, , ) scK = is defined as a sequence of
vessel rectangles.
BAPneighbor
1 Choose a rectangle k , where { 1, , 1 } kK ∈ −
2 Exchange the order of the pair ( , 1) kk + in the current state sc
3 Update to a new state sn by applying heuristics H1 and H2
Heuristics H1 and H2 are described later. A probability of making the transition from the
current state to a candidate state is specified by an acceptance probability function. We adopt
the following acceptance probability (AP) function, as described in [35].
()
1,
(, , )
exp( )/ , otherwise
en ec
AP ec en T
ec en T
< ⎧
⎪
=
⎨
−
⎪
⎩
(33)
where parameter T is a designed parameter (temperature in metallurgy).
Then, BAPsa implements the simulated annealing heuristic for the BAP, starting from the
state generated by the initialization routine BAPinit. The call BAPneighbor(sc) generates a
randomly chosen neighbor for a current state sc . The annealing schedule is defined by calling
(, , ) AP ec en T , which use the temperature to apply, given the fraction r of the time budget
that has been expended so far.
BAPinit
1 Sort vessels according to their arrival times so that
1 kK
aa a ≤ ≤≤ ≤
2 Generate a current state sc by applying H1 and H2
3 set ( ) ec bc E sc ==
27
BAPsa
1 set 1 i = , 1 flag = , and
0
TT =
2 while
max
ii ≤ and 1 T ≥
3 for 1:2 jK =
4 ( ) snsc = BAPneighbor , () en E sn =
5 if en eb < ,
6 then sbsn = , eb en = , and 1 flag =
7 elseif ( , ,) (1) AP ec en T rand > ,
8 then scsn = and ec en =
9 end
10 end
11 set TrT =
12 if 1 flag = ,
13 then 1 ii =+ and 0 flag =
14 end
15 end
The choice of the initial temperature
0
T and the cooling rate r affects the quality of the
solutions obtained by the simulated annealing procedure. To be able to determine the right
values for
0
T and r , we generated ten instances of the BAP problem whose sizes are
randomly varied from 16 to 21 vessels. Details about generating an instance of the BAP are
given in a later section. Each BAP was solved 700 times using different combinations of
0
T
and r . In our experiments,
0
T was varied from 10 to 130 by a step of 10, whereas r was
changed from 0.5 to 0.95 by a step of 0.05. For each
0
T and r , the results were averaged
28
and shown in Figure 2.7. Each normalized value represents the ratio of the corresponding
objective value to the lowest objective value. On average, the minimum objective value was
found at the initial temperature 70 and the cooling rate 0.6. These values are used in the
succeeding experiments.
0 20 40 60 80 100 120 140
1.06
1.07
1.08
1.09
1.1
Temperature
Objective value
Average normalized objective value
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95
1.04
1.06
1.08
1.1
1.12
Cooling rate
Objective value
Figure 2.7: Normalized objective values
Figure 2.8 shows the average number of iterations to reach final solutions (a) and the
number of minima (b). Note that the minima of each combination were selected over 10
repeated tests as described above. As seen from Figure 2.8a, a combination of a higher
temperature and a higher cooling rate, generally, yields better solutions. However, such
combination requires higher number of iterations too, which means more running times. The
actual required running time of a certain combination can be deduced from the figure. It shows
that our choice (
0
70 T = and 0.6 r = ) needs less than 300 iterations to reach to final solutions.
Also, the quality of solutions with respect to the parameters can also be estimated from Figure
2.8b. It indicates that our chosen parameters yield a higher number of minima (about 6.5).
Therefore, the chosen parameters are expected to produce good solutions within reasonable
time.
29
Cooling rate
Temperature
Average number of iterations to reach final solutions
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9
20
40
60
80
100
120
180
200
220
240
260
280
300
320
340
Cooling rate
Temperature
The number of minimums
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9
20
40
60
80
100
120
3.5
4
4.5
5
5.5
6
6.5
7
7.5
Figure 2.8: Average number of iterations and the number of minima
2.3.3. Heuristic method for an initial feasible solution
Although methods based on the subgradient optimization with Lagrangian relaxation and
methods based on simulated annealing have been used in the literature to solve the BAP and
similar allocation problems, our heuristic procedures are quite different from any procedures
reported previously in the literature. As already mentioned, our heuristics implement a set of
systematic and efficient methods in order to find an initial feasible solution and to update a
current solution by exploring a sufficiently large solution space.
To find and update feasible solutions for both subgradient optimization (SG) and
simulated annealing (SA) algorithms, two heuristic methods are developed. The first method
is called Heuristic H1 and is developed to generate an initial feasible solution. Using H1, each
vessel rectangle is sequentially placed on the solution space according to its order.
30
Heuristic H1
1 set
1
1 x =, 1 p = , and { 1 }
p
S =
2 for 2: kK =
3
11 kk k
xx l α
−−
=+ +
4 if
kk
x lM +≤ , then
pp
SS k = ∪
5 else, 1
k
x =, 1 p p =+ , and create { }
p
Sk =
6 end
7 for 1: kK =
8 if
1
kS ∈ , then
kk
ta =
9 else, max{ , }
kkss
tath τ =++ for all
1
k
p
s S
−
∈
10 end
where
p
S the set of rectangle indices in a batch p , where 1, , p B = , and | |
pp
Sq =
1
k
p
S
−
the subset of
1 p
S
−
, satisfying
sk k
xx l α < ++ and
ssk
xlx α + +> for
1 p
s S
−
∈ and
p
kS ∈ . That is, s represents a rectangle in
1 p
S
−
whose berthing space overlaps with
rectangle k .
We define a batch as a set of rectangles where the sum of their lengths, including safe
distances is less than or equal to the berth length M . In 1-6, a batch of rectangles is
constructed by assigning their berthing locations. In 7-10, a feasible berthing time is assigned
with the earliest possible berthing time.
31
2.3.4. Heuristic method for an improved feasible solution
Heuristic H2 is developed to repeatedly improve the feasible solution, using pair-wise
swaps between neighboring rectangles within a batch or over neighboring batches. Note that a
batch is defined as a set of rectangles where the sum of their lengths including safe distances
between them is less than or equal to the berth length M. Under any combination of rectangle
movements and swapping operations, Heuristic H2 tries to create rooms (temporal spaces) for
rectangles to be berthed to the earliest possible time (pushed down). That is, H2 uses a series
of rectangle movement functions, which consist of all possible directional movements and
swapping operations.
Heuristic H2
1 set
2 H
z with the minimum of known feasible solution value
2 find
2
new
H
z by running SWAPb, SWAPt, RIGHT, SWAPb, SWAPt, and LEFT
3 update and return
2 H
z if
22
new
HH
zz ≤
The rectangle movement functions are designed to generate a better solution by moving
each rectangle down to its allowable limit which is defined by the arrival time or the
neighboring rectangles.
SWAPb (SWAPt) function exchanges a pair of neighboring rectangles in the same batch,
()
p
Sq and ( 1)
p
Sq + , which starts from the bottom batch (or, from the top batch). These
functions eventually try to swap all the possible combinations of rectangle pairs using the
repetition scheme while Swap-Pair operator is designed to exchange of the positions of
consecutive rectangles.
1 for 1: p P = (:1 p P = )
2 for 1: 1
p
qq = −
32
3 if p P < , then
4 Push-Up
1 p
kS
+
∈ , Swap-Pair q and 1 q + , Push-Down
1 pp
kS S
+
∈ ∪
5 else,
6 Swap-Pair q and 1 q + , Push-Down
p
kS ∈
7 if the cost increases, then restore their positions
8 end
9 end
RIGHT (LEFT) function pushes a rectangle ( )
p
Sq to the allowable limit or to the next
nearest rectangle on the right (or, on the left) in the same batch.
1 for 1: p P =
2 for 1
p
qq = −
3 if p P < , then
4 Push-Right q (Push-Left q ) and Push-Down
1 pp
kS S
+
∈ ∪
5 else,
6 Push-Right q (Push-Left q ) and Push-Down
p
kS ∈
7. if the cost increases, then restore their positions
8 end
9 end
The following rectangle movement operators define all the possible movements of a
rectangle in time-space and are the basis of the rectangle movement functions.
Push-Up Moves a rectangle k up to the allowable limit
Push-Down Moves a rectangle k down to the allowable limit
Push-Right Moves a rectangle q right to the allowable limit
33
Push-Left Moves a rectangle q left to the allowable limit
Swap-Pair Swap consecutive rectangles q and 1 q + (
1 qq
x x
+
< ) if space is allowable
As some examples, the following figures show how the combination of several rectangle
movement operators can reduce the allocation cost. Each allocation plan for a vessel is
represented by a shaded rectangle along with its arrival time and requested departure time
(surrounding dotted box). Also, darker shade represents the region which causes a penalty cost
due to the violation of the desired departure time.
()
p
Sq
(1)
p
Sq +
1
(1)
p
Sq
+
′ +
1
()
p
Sq
+
′
1
()
p
Sq
+
′
1
(1)
p
Sq
+
′ +
()
p
Sq
(1)
p
Sq +
Figure 2.9: Reduction by SWAPx resulting from Swap-Pair followed by Push-Down
()
p
Sq
(1)
p
Sq +
1
(1)
p
Sq
+
′ +
1
()
p
Sq
+
′
()
p
Sq
(1)
p
Sq +
1
(1)
p
Sq
+
′ +
1
()
p
Sq
+
′
Figure 2.10: Reduction by RIGHT resulting from Push-Right followed by Push-Down
34
()
p
Sq
(1)
p
Sq +
1
(1)
p
Sq
+
′ +
1
()
p
Sq
+
′
()
p
Sq
(1)
p
Sq +
1
(1)
p
Sq
+
′ +
1
()
p
Sq
+
′
Figure 2.11: Reduction by LEFT resulting from Push-Left followed by Push-Down
35
2.4. Computational Experiments for the BAP
A disruption at a berth may result in changes in the space-time diagram as compared to
the baseline plan. The disruption may alter the diagram in the time dimension, the space
dimension, on in both dimensions. Usually, disruptions in the time domain are caused by
delayed arrival times, delayed berthing times, or extended service times. Disruptions in the
space domain are less frequent yet their impact is more severe on the port operations. They
may be the result of anticipated events (e.g., construction, scheduled maintenance, pre-planned
military surge, etc.), or unanticipated events (e.g., terrorist acts, earthquakes, hurricanes, etc.).
In this section, we will consider disruptions caused by delays, i.e., disruptions affecting
only the time dimension of the space-time diagram. Disruptions affecting the space dimension
will be studied in later section, where we introduce partially operational berths and the
multiple bertha allocation problem. Here, we assume that the delays are caused by an increase
in the number of calling vessels on the berth. This increase in the number of vessels to be
serviced could be for example a direct consequence of physical break-downs in adjacent
berths within the same terminal/port, due to unanticipated events. Because of the increased
number of vessels, we may not be able to serve all the vessels within their requested time
frames. The objective is to minimize the total delay incurred for all the vessels. We will show
that the developed optimization methods are able to deal with this type of disruptions by
finding the best allocation for the calling vessels such that the total delay is minimized.
2.4.1. CPLEX MIP: Computational experiments
The CPLEX MIP solver is used to find exact solutions for small size problems, enabling
us to evaluate the quality of the solutions generated by our methodologies for such problems.
The computational experiments are conducted using realistic data obtained from Park and Kim
36
[42], and Kim and Moon [34]. In order to generate random instances, we use a discrete
uniform distribution whose cumulative distribution function (CDF) is given by
1
1
(; , ) ( )
1
n
i
i
Fx ab H x x
ba
=
= −
−+
∑
(34)
where the Heaviside step function ( )
i
Hx x − is the CDF of the degenerate distribution
centered at
i
x . Table 2.1 shows the range of parameters a and b assumed in our
simulation experiments to generate
k
l ,
k
h ,
k
a , and
k
d . For instance,
k
h is generated
randomly using the given CDF for 7 a = and 23 b = . Upon choosing
k
h , the value of
k
a is
randomly generated using the same CDF with parameters 1 a = and 1
k
bT h =− + . Note that
k
d is chosen so that 2
kk k k
hd a h ≤− ≤ . In our simulation experiments, we set 1200[ ] M m = ,
168[ ] Th = , 20[] m α = , and 2[ ] h τ = .
Table 2.1: Ranges of parameters
Parameter
k
l [m]
k
h [h]
k
a [h]
k
d [h]
a 170 7 1
kk
ah +
b 290 23 1
k
Th − + 2
kk
ah +
We assume that the weight
k
w in (1) is equal to the area covered by the vessel rectangle,
i.e.,
kk k
wl h =⋅ . This choice of
k
w imposes higher penalties to be paid to ships with possibly
higher workloads, if the requested departure times are not met. Using the specifications given
in Table 2.1, we generate 500 random instances of the BAP. Table 2.2 compares the
performance of the CPLEX MIP solver (IP), the subgradient optimization method (SG), and
the simulated annealing method (SA). The CPLEX MIP solver is used to find the exact
solutions, which enable us to evaluate the quality of the solutions generated by the developed
heuristic methods.
37
Table 2.2: Comparison between the BAP solutions using IP, SG, and SA
Cost Average time [s] Maximum time [s]
K r
OC
z
IP
z
SG
z
SA
z
IP
/z
SG
z
IP
/z
SA
t
IP
t
SG
t
SA
m
IP
m
SG
m
SA
11 0.19 10.7 32.4 14.3 3.03 1.34 0.0 2.0 2.1 10.6 115.4 178.0
12 0.21 23.5 46.7 41.5 1.99 1.77 0.2 3.4 6.6 25.4 130.8 254.1
13 0.22 34.0 112.8 52.8 3.31 1.55 0.8 5.5 8.5 236.5 106.2252.0
14 0.24 71.2 187.1 115.9 2.63 1.63 1.7 11.2 20.0 468.8 183.8444.6
15 0.26 89.8 310.0 170.1 3.45 1.89 28.8 13.7 21.5 6781.0 238.1367.8
16 0.27 89.7 299.2 137.6 3.34 1.53 16.7 19.1 31.0 1486.2 234.2481.8
17 0.29 230.4 606.4 307.6 2.63 1.33 156.8 30.5 61.8 24739.6 260.4 639.6
18 0.31 244.0 683.4 389.7 2.80 1.60 361.1 38.7 76.9 38798.5 239.3 617.6
In Table 2.2,
IP
z ,
SG
z , and
SA
z are, respectively, the average cost of the best solutions
found by the IP, SG and SA. The average and the maximum running times for different
methodologies are also shown in the table. It indicates that the cost ratio of SA to IP (i.e.,
/
IPSA
zz ) is rather small which means that SA yields good quality solutions. The cost ratio
/
IPSG
zz is much higher than /
IPSG
zz . To improve the solution quality of SG, one can
increase the number of iterations of H2.
We define the occupation ratio of a berth,
OC
r , as the sum of the vessel rectangles divided
by the entire space-time domain. We use
OC
r as an indicator of the berth utilization. As
indicated by the results,
OC
r is increased, the time to find the best solutions increases
exponentially
1
()
K
OC k k
k
rhlMT
=
=
∑
(35)
The IP can solve the problem instances up to K=18. However, comparing to the heuristics,
as K increases, the IP needs considerable amount of time to find the exact solution. We noticed
that for K>17, we will not be able to find the exact solutions in a reasonable amount of time.
38
This limitation will be more pronounced as we move toward the dynamic BAP and the
multiple BAP, which are discussed in the following sections.
2.4.2. Heuristic methods: Computational experiments
In order to further evaluate the developed heuristic method, we consider instances of the
BAPs with 17 21 K ≤≤ . This range is similar to the number of vessels considered in real
practices [42],[34]. We also adopt the generalized objective function
3
J in (24) with
1
10
kk
wl = and
2
kk
wl = . These values are chosen so that we can impose higher penalties on
delayed departure times. The values of
1
k
w and
2
k
w are directly proportional to the length of
the vessel k ,
k
l . Hence, higher priority is given to larger vessels which, most probably, carry
higher volume of loads. Table 2.3 shows the performance of the SG and SA optimization
methods. Each row of the table is the averaged result of 20 independent trials whose initial
feasible solutions are non-zero.
Table 2.3: Performance comparison between SG and SA methods
Subgradient optimization Simulated Annealing
K r
OC
z
LB
z
UB
g
SG
t
SG
c
SG
z
SA
t
SA
c
SA
z
SA
/z
UB
17 0.29 145.2 380.3 1.6 590.4 227.1 322.1 133.3 148.1 0.85
18 0.32 403.2 737.4 0.8 742.5 281.6 532.8 183.1 195.4 0.72
19 0.33 525.8 897.9 0.7 711.0 268.6 569.3 154.7 183.4 0.63
20 0.35 428.6 899.7 1.1 845.4 304.8 545.9 261.1 289.0 0.61
21 0.36 729.5 1951.2 1.7 822.8 268.61409.1248.5 261.4 0.72
In the table,
OC
r
is the occupation ratio of a berth,
LB
z is the maximum lower bound
from Lagrangian relaxation model,
UB
z is the minimum upper bound (best feasible solution
by the SG method),
SG
g is the duality gap between
LB
z and
UB
z , defined as
()/
SG UB LB LB
gz z z =− ,
SG
t is the computational time of the SG method,
SG
c is the number
39
of iterations of H2 in the SG method,
SA
z is the best feasible solution by the SA method,
SA
t
is the computational time of the SA method, and
SA
c is the number of iterations of H2 in the
SA method. The duality gap
SG
g implies the maximum deviation of the final lower bound
from the best feasible objective value.
The duality gap ranges from 70% to 170% which means that the SG produces good
approximations of optimal solutions. In order to improve the solution quality of the SG
methodology, we may increase the number of iterations of H2. It is noted that step 4 in BAPsg
doesn’t need to be executed every iteration. We noted that the step should be executed when
LB
zmaxLB < or
LB
zmaxLB ≤ . We can also observe that, on the average, the SA algorithm
still yields a non-inferior solution, which ranges from 61% to 85% of the corresponding upper
bound found by the SG algorithm. This might have resulted from the fact that the range for
swapping operations in the SA method is little wider. Considering the superiority of the
solution quality and the running time of the SA methodology, in the following sections we will
only use the SA methodology to find the best solutions to other variants of the BAP namely
the dynamic BAP, and multiple BAP.
40
Chapter 3 Port Level
3.1. Mitigating disruption at port level
The objective at the port level is to route goods to different terminals within the port so
that the overall port throughput is affected less by disruption. When a disruption takes place,
the terminals might be unable to meet the expected demand, due to partial or total loss of
operational capabilities, or to a sharp rise in the demand (e.g. military surge). Furthermore,
such a disruption could affect several terminals within a particular port. It is therefore critical
to allocate berths to ships in such a way as to meet all demand, minimize the vessel berthing
time and maximize berth utilization. In order to mitigate the impact of disruptions, methods
will be developed to re-route ships to different berths within the terminal, or to different
terminals within the port so that the overall port throughput is affected as little as possible.
The question we answer in this level is how we can reassign ships to berths/terminals
within the same port, such that the overall port throughput is maintained. Therefore, we focus
on the scheduling problem of calling vessels on multiple berths which is referred to as the
multiple berth allocation problem (mBAP). We will see that the mBAP can be viewed as a set
partitioning problem, and hence it is an NP-hard problem.
3.2. Multiple berth allocation problem
Usually certain berthing locations (home berths) are preferred due to the long-term
contracts with carriers, the depth of water, differing wave heights, etc [42]. In this section, we
assume that for some predicted or unpredicted scenarios (disruptions) a calling vessel cannot
moor at its home berth location and that other berth locations (within the same terminal or
adjacent terminals) can accommodate the vessel. This leads to a more complex, yet general,
41
variant of the BAP which, hereafter, is referred to as the multiple berth allocation problem
(mBAP).
We assume that we have N disjoint berths. The disjoint berths could belong to the same
or to different terminals, and because they are disjoint, a vessel can utilize only one such berth
(i.e. a vessel cannot moor across two or more such berths). The disjoint property will enable us
to model a partially functional terminal (e.g., a terminal consisting of N disjoint berths, of
which
1
N N ≤ have been rendered non functional due a disruptive event, hence there are
only
1
N N − remaining functional berths). We assume that the length of berth n is
n
M
and that the time horizon for all berths is T . Figure 3.1 illustrates the TAP graphically in a
space-time diagram.
T
1
M
T
2
M
T
N
M
• ••
Figure 3.1: Space-time diagram for the Terminal Allocation Problem
The mBAP can be stated as follows: Determine the least-cost assignment of K vessels
to N disjoint berths such that each vessel is assigned to exactly one berth and no two vessels
are overlapping (in the space-time domain), while all the vessels’ constraints (arrival times,
service times, etc.) are met. Therefore, the mBAP consists of two intertwined problems: (a) A
set partitioning problem (SPP), and (b) A number of individual berth allocation problems
(BAPs).
Here, we briefly visit the definition of the set partitioning problem (SPP). Given a
collection of feasible subsets of a certain ground set, one can formulate the problem of finding
42
the best collection of subsets such that the cost associated with these subsets is minimized. Let
j
y be 1 if berth j is functional and at least one vessel is assigned to it, and zero otherwise.
Let
j
c be the cost of berth j which is computed as the aggregated costs of serving all the
vessels assigned to that berth. Let also
ij
a be 1 if vessel i is assigned to berth j , and zero
otherwise. The SPP can be formulated as follows:
min
1
N
jj
j
cy
=
∑
(36)
s.t.
1
1
N
ij j
j
ay
=
=
∑
, 1, , iK = (37)
{0,1 }
j
y ∈ , 1, , j N = (38)
Almost every heuristic approach, for solving general integer programming problems, has
been considered and applied to the set partitioning problem. For instance, in [17], some greedy
algorithms and interchange approaches were applied to the SPP problem. Other heuristic
approaches, such as genetic algorithm, probabilistic search, simulated annealing, and neural
networks, have also been studied and applied to the problem. In addition, heuristics are
embedded within exact algorithms so that one iteratively tightens the upper bound and at the
same time attempting to get a tight approximation to the lower bound for the problem [26],[3].
Unfortunately, there has not been a comprehensive and comparative test across such methods
to determine under what circumstances a specific method might perform best.
It should be noted that the SPP is known to be NP-hard [19]. As described above, the TAP
can be mapped into the SPP; hence, the mBAP is an NP-hard problem too. Moreover, the
mBAP consists of a collection of BAPs, where each BAP is an NP-hard problem [11].
Therefore, as the number of the berths and vessels in the mBAP increases, finding an optimal
solution will be computationally too expensive. In the following, we develop an
43
approximation methodology to find good solutions to the mBAP in a reasonable amount of
time.
3.3. Solution method for the mBAP
The approach is based on a typical set partitioning procedure. Assuming that the terminal
consists of N disjoint berths, the TAP is partitioned into N BAPs. Each BAP is solved
separately using the previously developed heuristics. First, we define a vessel rectangle k ,
denoted by
k
v , by a 9-tuple as follows,
( )
12
,, , , , , , ,
kkkkkk k kkk
vnlhadwwxt = (39)
where
k
n is the assigned berth index for vessel rectangle k
k
l is the length of vessel k
k
h is the estimated handling time of vessel k
k
a is the estimated arrival time of vessel k
k
d is the requested departure time of vessel k
1
k
w is the penalty weight assigned to vessel k
2
k
w is the reward weight assigned to vessel k
k
x is the assigned mooring location to vessel k
k
t is the assigned mooring time to vessel k
Let
n
B be the set of vessels assigned to berth n . That is,
( ) { } { } | , , , , 1, , , 1, ,
nk k k k k
vv n t n n k K n N == = = ∈ …… B
where K is the number of vessels to be scheduled, and
nn
K = B is the cardinality of the set
44
representing the number of vessels assigned to berth n. A state s is formed by aggregating
all
n
B in a single set, { }
1
,,
N
= s BB .
The cost of the state s is defined as the sum of the costs of all the individual berths
n
B
forming s . Recall that the cost of each berth is defined in (1) (or in (23) for the generalized
cost). The proposed SA algorithm uses the following stats and costs:
ec cost of the current state sc, () E = ec sc
en cost of the newt state sn, () E = en sn
eb cost of the best state sb, () E = eb sb
In the set partitioning algorithm, we construct a single aggregated state sc to easily
perform an interchange operation. The new state sn is obtained from the current state sc
by exchanging a pair of vessel rectangles across different berths. After interchanging vessels
across different berths, the optimization procedure is sequentially applied to the sub-problems.
To choose the candidate vessel rectangles, overlapping indices are introduced. Recall
from Figure 2.4 that a block ( , ) ij is referred to as a grid space whose left bottom corner is
located at the point ( , ) ij . The value of block ( , ) ij at berth n, denoted by
ijn
c , is defined
as the number of overlapping vessels at block ( , ) ij at that berth. These overlaps are resulted
from the occupancy of more than one vessel rectangle when the assigned berthing times are
replaced with the estimated arrival times.
Let
{ }
max,
,
max
nijn
ij
cc = , 1, ,
n
iM = and 1, , jT = . The overlapping index of vessel k ,
denoted by
k
IV , is defined as the number of blocks whose values are
max,n
c within the vessel
rectangle k . The overlapping index of berth n , denoted by
n
IB , is defined as the sum of the
block values in that berth.
45
11
n
M T
nijn
ij
IB c
==
=
∑ ∑
(40)
Figure 3.2 illustrates a simple example to demonstrate how to find a vessel’s and a berth’s
overlapping indices. The right figure shows berth n containing two overlapping vessel
rectangles k and k ′ . For this example
max,
1
n
c =, 2
n
IB = , and 2
kk
IV IV
′
= = .
k
a
k
t
k
d
k
k ′
Figure 3.2: Value of the block
The policy for choosing the pair of vessels to be exchanged is represented in lines 2 to 5
of the following procedure.
NewPartition
1 order the rectangles of a current state sc by arrival times.
2 choose (randomly when tie) a berth
1
b having max
n
IB .
3 choose (randomly when tie) a vessel
1
k having max
k
IV where
1
b
k ∈ B .
4 choose the first
2
k such that
12
kk
aa < and
12
kk
nn ≠ .
5 if
2
k cannot be found, go to step 2
6 else, exchange the berth indices of
1
k
v and
2
k
v .
7 divide into sub-states { } |
nk k
vn n ′ = = B and sort by arrival times.
8 generate an initial solution for each BAP by applying H1 and H2.
9 update each state ( )
nn
′ = BAPsa BB and return a new state { }
1
,,
N
= sn BB .
46
The proposed solution procedure for the mBAP is referred to as TAPproc. In the
procedure, the sum of
n
IB s under ith partition is denoted by
i
sIB . This procedure is
terminated when the best objective function value of the mBAP or the sum of the berth
overlapping indices
i
sIB reaches zero.
TAPproc
1 initialize = bc sc by TAPinit, () E = = eb ec sc, 1 i = , and 0 r = .
2 while
max
ii ≤
3 ( ) = sn sc NewPartition , () E = en sn , and
in
n
sIB IB =
∑
4 if < en eb , then = sb sn and = eb en ;
5 elseif 0 = eb or 0
i
sIB = , then stop
6 for 1: 1 ji = −
7 if
j i
sIB sIB = , then 1 rr = +
8 end
9 if
max
rr = , then stop
10 end
This algorithm starts by a feasible solution generated by TAPinit. Vessels are alternately
assigned to berths according to their arrival times. However, in the presence of a partially
functional berth (which will be defined later), the initial subsets will be assigned based on the
ratio of the operational region in the space-time domain.
TAPinit
1 sort the order of rectangles of a current state sc by arrival times
2 assign a berth index to each vessel (mod )
k
nk N ≡
3 divide into sub-states ( ) |
nk k
vn n ′ = = B
47
4 generate an initial solution by applying H1 and H2
5 update each state ()
nn
′ = BAPsa BB
3.4. Consideration of partially functional berth (terminal)
A berth is said to be functional (or fully functional) if the entire berth is available to be
used by the calling vessels for the entire planning horizon. A berth is partially functional if
some sections along the berth are unavailable for some periods of time. These non-operational
regions may be caused by anticipated events (e.g., construction, scheduled maintenance, pre-
planned military surge, etc.), or by unanticipated disruptions (e.g., terrorist acts, earthquakes,
hurricanes, etc.). A partial region in the space-time diagram is termed a cell. A non-operational
region is a non-operational cell. The remaining region in the space-time diagram can be
divided into several operational cells based on its shape. These cells are all mutually disjoint
and exhaustively cover the space-time diagram. To decompose the time-space diagram into
operational and non-operational cells, we employ a simple coverage type algorithm [10].
Since both types of our cells are rectangular in shape, hence we can easily employ a grid-
based coverage algorithm.
Subsequently, our developed heuristics are modified to incorporate the terminal allocation
problem when a berth is only partially functional. A cellular decomposition scheme is
incorporated into the heuristic H1, which was developed to generate an initial feasible solution.
The modified heuristic is referred to as the heuristic mH1. Once the vessels are assigned to
batches, the heuristic H2 will be deployed to improve the current solution.
The modified heuristic for an initial feasible solution is done in two stages. At first, the
solution space is roughly decomposed into cells that allow some overlapping with neighboring
operational cells. Figure 3.3 shows some examples of the decomposed cells from the first
stage. The space-time diagrams are divided into several operational and non-operational cells.
48
The cells are mutually disjoint and exhaustively covering the space-time diagram. An
operational cell r is defined by the coordinate of its lower-left corner and the projected
distances of its upper-right corner along both axes: ( ,, ,
rr r r
x tM T ), where the total number of
operational cells is R .
1
M
2
M
1
T
2
T
Operational
cell 1
Operational cell 2
Non-
operational
cell
M
T
1
2
M
T
2
3
1
Non-
operational
cell
Non-
operational
cell
1
2
11 1 1
(, , , ) x tM T
1
1
x
Figure 3.3: Examples of partially functional berths
In the second stage, any overlapping is removed since the upper cells are adjusted by the
allocation plan of lower cells. Then, the cells are mutually disjoint and exhaustively cover the
solution domain. In the modified procedure for an initial feasible solution, the cells are
decomposed and allocated in two stages by the following procedure in a bottom-up fashion.
mH1
1 define the operational cells with ( , , , )
rr r r
x tM T for 1, , rR =
2 set
1
1 k = and 1 p =
3 for 1: rR =
4 set
r
kk = ,
r
k
x x = , and { }
p
Sk =
5 for 1:
r
kk K =+
6 set
11 kk k
xx l α
−−
=+ +
7 if
r
kk
x lM +≤ , then { }
pp
SS k = + ;
49
8 else,
r
k
x x =, 1 p p = +, {}
p
Sk =
9 end
10 for :
r
kk K =
11 if
1
kS ∈ , then
kk
ta = ;
12 else, max{ , }
kkss
tath τ =++ for all
1
k
p
s S
−
∈
13 end
14 for :
r
kk K =
15 if
r
kk
tl T τ ++ > , then { }
pp
SS k ′ = − for all kk ′ ≥, 1 p p =+ ,
1 r
kk
+
=
16 end
17 end
where
1
k
p
S
−
is a subset of
1 p
S
−
, satisfying
sk k
xx l α < ++ and
ssk
xlx α ++ > for
1 p
sS
−
∈ and
p
kS ∈ .
In lines 4-9, each batch of rectangles is constructed by assigning their berthing locations.
The first rectangle of operational cell r is denoted as
r
k and is allocated at the beginning
section of the cell,
r
x . In lines 10-13, a feasible berthing time is assigned as the earliest
possible berthing time. Finally, in lines 14-16, a batch that overrides with the next cell is
determined as the last batch of the current cell. The contour of the next cell is adjusted
according to the allocation plan of the current cell. Eventually, the entire functional region can
be decomposed into independent cells occupied by a feasible berthing plan. If an overlapping
takes place in the last batch, it implies that the remnant of the current scheduling will stay for
the next planning horizon. However, in a rolling horizon framework, we can easily generate a
feasible solution by considering the next planning horizon as a neighboring cell.
50
Therefore, in the case of a partially functional berth, the heuristic mH1 replaces H1 in
procedures BAPneighbor, BAPinit, and BAPsa. Subsequently, the procedures for the mBAP
should be properly modified to handle a situation where a partially functional berth exists
among multiple berths. We assume that all ships can be allocated to any berth, possibly,
allowing for extra penalty if the ship does not moor at its home berth. Again, in the presence
of a partially functional berth, the heuristic mH1 replaces H1 in the procedure NewPartition.
The initialization procedure for the TAPproc, is modified so that vessel rectangles are
distributed according to the proper occupation ratio in the initialization procedure, which is
referred to as mTAPinit. In the presence of a non-operational region, the occupation ratio of
berth n, ()
OC
rn , is redefined as the sum of the vessel rectangles and the non-operational area
of
n
A divided by the entire domain.
()
( ) ( )
kn
OC k k n n
kv
rn hl A MT
∈
=+
∑
B
(41)
mTAPinit
1 sort a current state { }
1
,,
N
= sc BB by arrival times
2 assign a berth index to each vessel: (mod )
k
nk N ≡
3 randomly substitute berth indices so that
()
1
1()
1()
OC
n N
OC
n
rn
KK
rn
=
−
=
−
∑
where
1
N
n
n
KK
=
=
∑
4 divide into sub-states { } |
nk k
vn n ′ = = B
5 generate a feasible solution by applying mH1 and H2
6 update each state ( )
nn
′ = BAPsa BB
51
3.5. Computational experiments
Several computational experiments are performed to evaluate our developed methodology
for the mBAP. Our computational experiments consist of both functional and partially
functional berths. Here, we consider vessels with both flexible and inflexible berthing
locations and departure times. To evaluate our develop techniques, under different
experimental scenarios, we use the cost function
6
J defined as
3
63
1
min
K
kk k
k
wn m
=
⎡ ⎤
=+ −
⎢ ⎥
⎣ ⎦
∑
JJ , (42)
where
3
J is the cost function defined in (23),
k
m the desired (home) berth location for
vessel k , and
3
k
w the spatial penalty for vessel k . The penalty
3
k
w is applied to a vessel
which cannot be moored at its desired berth. The last term in
6
J , in fact, penalizes the
difference between the assigned and the desired berthing location for vessel k . Therefore, we
redefine a vessel rectangle k , denoted by
k
v , by a 11-tuple as follows,
( )
12 3
,, , , , , , , , ,
kk kkkkk k k kkk
v n m l h a d w wwx t = (43)
In the following computational experiments, and without loss of generality, we assume
the following values for penalty weights:
1
, if 1 and
10 , otherwise
d
kkkk
k
k
Qf thd
w
l
⎧ = +>
⎪
=
⎨
⎪
⎩
(44)
2
kk
wl = (45)
3
, if 1 and
2, otherwise.
b
kkk
k
k
Qf nm
w
l
⎧ =≠
⎪
=
⎨
⎪
⎩
(46)
where
Q is a big positive real number
52
1, if the departure time of vessel is not flexible
0, otherwise
d
k
k
f
⎧
=
⎨
⎩
1, if the berthing location of vessel is not flexible
0, otherwise
b
k
k
f
⎧
=
⎨
⎩
3.5.1. Vessels with flexible berthing locations and departure times
For the first category of scenarios, it is assumed that the calling vessels have flexible
berthing locations and flexible departure times (i.e., 0
d
k
f = and 0
b
k
f = ). By flexible
berthing locations, we mean that the vessels do not have any restrictions on the section of the
berth they will be assigned to. However, these vessels may be willing to moor at any other
berth within the terminal or may only choose a specific berth. For the remaining experiments,
and without loss of generality, we assume that 3 N = , 1200
n
MMm = = for 1, ,3 n = ,
and 168 Th = . We also assume that the number of vessels mooring at each berth is initially
uniformly distributed between 16 to 21.
Scenario 1: In this scenario, all three berths are functional and the vessels are not flexible
to change their berths. In fact, this scenario can be viewed as three disjoint BAPs.
Scenario 2: In this scenario, the vessels in Scenario 1 are flexible in changing their berths.
This scenario, in fact, represents an instance of the mBAP with three berths.
Twenty instances of the problem are randomly generated with general specifications as
described above. Each instance is comprised of three BAPs such that at least one BAP has a
non-zero initial objective function value. In order to evaluate the proposed heuristics for
Scenarios 1 (BAPs) and 2 (mBAP), a series of BAPprocs and TAPproc are successively
applied to each instance. The results are averaged and summarized in Table 3.1.
In Table 3.1,
OC
r is the occupation ratio of the terminal which is defined as the statistical
mean of ( )
OC
rn in (41).
Si
z is the value of the best solution in Scenario i ,
Si
t the
53
computational time in seconds,
Si
c the number of iterations of procedure H2 by calling
BAPneighbor,
Si
p the number of iterations of the set partitioning by calling NewPartition,
and
ij
g the cost reduction as a solution gap between Scenarios i and j , i.e.,
()/
ij Si Sj Si
g zz z =− . Due to the flexibility of vessels in choosing different berths, the total
allocation cost of Scenario 2 is reduced by 99% compared to that of Scenario 1.
Table 3.1: Scenarios 1 and 2, experimental results
Scenario 1 Scenario 2
r
OC
z
S1
t
S1
c
S1
z
S2
t
S2
c
S2
p
S2
g
12
0.31 2796.0 967.7 491.8 27.3 7700.7 4558.0 13.5 0.99
The following scenarios 3 and 4 are similar to scenarios 1 and 2 under the additional
condition that one of the berths in the terminal is now partially functional.
Scenario 3: Consider Scenario 1 again. Here, berths 2 and 3 are fully functional while
berth 1 is partially functional. We assume that berth 1 has a non-operational cell between 600-
1200m in space and 1-100h in time in the space-time diagram.
Scenario 4: Consider Scenario 2 again. Similar to Scenario 3, we assume that berths 2
and 3 are fully functional while berth 1 is non-operational between 600-1200m in space and 1-
100h in time.
Using the same problem instances generated for Scenarios 1 and 2, Scenarios 3 and 4 are
created by assuming a non-operational cell in first berth, respectively. Comparing Table 3.1
and Table 3.2, we can easily observe that the non-operational cell increases the occupation
ratio from 0.31 to 0.41 and it increases the cost of BAPs by more than 10 times. It indicates
that, even in the presence of a partially functional berth, the allocation cost can still be reduced
by using multiple berths. It is found that the cost is decreased by 93%.
54
Table 3.2: Scenarios 3 and 4, experimental results
Scenario 3 Scenario 4
r
OC
z
S3
t
S3
c
S3
z
S4
t
S4
c
S4
p
S4
g
34
0.41 32232.0 1241.0 504.8 2266.8 41036.8 26338.5 44.7 0.93
3.5.2. Vessels with strict berthing locations and departure times
For the second category of scenarios, we assume that the departure times and/or the
berthing locations of some vessels cannot be violated (i.e., 1
d
k
f = and/or 1
b
k
f = for some
vessels). Accordingly, the set of vessel rectangle
k
v assigned to berth n ,
n
B , is divided into
two subsets:
() { }
(){} {}
|,, , , 0, 1,, , and
| ,, , , 1, 1, , 1, ,
free b
nkk k k k k
fixed b
nkk k kk k
vv n t n n f k K
vv n t n n f k K n N
== = = =
== = = = ∈
……
…… B
B
(47)
The procedures for generating an initial feasible solution and an updated partition are
respectively modified to handle the inflexible requests by vessels. In mTAPinit, the rectangles
from the first subsequence
free
n
B are distributed over N berths according to the initial
occupation ratio ()
OC
rn and then the assigned rectangles are combined with the rectangles
from the second set
fixed
n
B . In NewPartition, only the rectangles from the set
free
n
B are
selected for swapping operation.
Scenario 5: We assume that all three berths are functional and that vessels are not flexible
in changing their berths. Additionally, we assume that 30% of the vessels in each berth have
strict departure time requests. In fact, this scenario represents three BAPs having strict
berthing time requests.
Scenario 6: In this scenario, we assume that some vessels in Scenario 5 are flexible in
changing their assigned berths, However, 30% of the vessels in each berth have requested
strict berth location. In fact, this scenario represents a mBAP having strict requests on berthing
times and berthing locations.
55
Each of 20 random instances is comprised of three BAPs such that at least one BAP has a
non-zero initial objective function value generated by mH1 or H1. The averaged results are
summarized in Table 3.3. Since the BAPs in Scenario 5 have the strict requests on departure
times, some vessels’ rectangles, probably with lower priority, may not be allocated in the
current time horizon (possibly placed over two consecutive planning horizons). The number of
these undesirable allocations in Scenario i is denoted by
Si
e . According to the results, if
vessels are allocated to the single integrated terminal with multiple berths, the allocation cost
can be reduced by 63% even with the strict requests on berthing locations and departure times.
Note that the average number of vessels having the inflexible departure times is 5.55 for 3
berths. These vessels having higher priority could make some allocations undesirable.
However, by implementing the mBAP, the number of undesirable allocations can also be
reduced from 0.27 to 0.09.
Furthermore, additional scenarios are prepared to compare the BAPs and mBAP in the
presence of a partially functional berth. Using the same problem instances generated for
scenarios 5 and 6, scenarios 7 and 8 are prepared by adding a non-operational cell in the first
berth, respectively.
Table 3.3: Scenarios 5 and 6, experimental results
Scenario 5 Scenario 6
r
OC
z
S5
t
S5
c
S5
e
S5
z
S6
t
S6
c
S6
p
S6
e
S6
g
56
0.31 31258.1 666.4 489.9 0.27 11442.2 8052.1 6916.6 10.8 0.09 0.63
Scenario 7: Consider Scenario 5 again. Assume that berth 1 has a non-operational cell
between 600-1200m in space and 1-100h in time while berths 2 and 3 are fully functional.
Scenario 8: Consider Scenario 6 again. Similar to Scenario 7, we assume that berths 2
and 3 are fully functional while berth 1 has a non-operational cell between 600-12000m in
space and 1-100h in time.
56
Table 3.4 shows the benefit of implementing the single integrated terminal in the presence
of a partially functional berth when vessels have both strict requests on berthing locations and
departure times. A comparison between Table 6 and Table 7 indicates that, although the
presence of the non-operational cell increases the cost of BAPs by more than 700%, the
allocation cost is still reduced to 57%. In addition, the number of undesirable allocations is
reduced from 1.91 to 0.82, on average.
Table 3.4: Scenarios 7 and 8, experimental results
Scenario 7 Scenario 8
r
OC
z
S7
t
S7
c
S7
e
S7
z
S8
t
S8
c
S8
p
S8
e
S8
g
78
0.41 221176.1 872.3 599.2 1.91 95818.3 21724.5 21444.5 33.5 0.82 0.57
57
Chapter 4 Regional Level
4.1. Service network design
Logistics service providers or carriers consolidate their freight in a network of hubs and
terminals and build up regular services. The design of such services requires decisions about
the frequency, mode, and route of the service and the corresponding schedule and routing of
freight. Planning decisions for such networks are made on a tactical level and have a direct
impact on customer services and costs. These kinds of tactical planning problems are
generally referred as a service network design [50].
Railways, less-than-truckload (LTL) motor carriers, mail/package delivery services,
airlines, intermodal shipping lines are typical examples of such systems that require tactical
planning of operations. Service network design is increasingly used to designate the set of
main tactical issues and decisions relevant to the described carriers: the selection and
scheduling of the services to operate, the specification of the terminal operations, the routing
of freight, the empty balancing strategy, etc. The corresponding models usually takes the form
of network design formulations, a class of the mixed-integer network optimization problems
for which no efficient exact solution method exists, except for special variants [13].
In this work, we describe a general service network design problem containing decisions
regarding the service selection and the traffic distribution. In the service selection, the routes
which services are offered and the characteristics of services are defined. In the traffic
distribution, the paths between origins and destinations are determined to move every demand.
Network design formulations are defined on graphs containing nodes connected by links.
The objective is to select links in a network, along with capacities, eventually, in order to
satisfy the demand for transportation at the lowest possible system cost computed as the total
58
fixed cost of the selected links, plus the total variable cost of using the network. The fixed cost
service network design can be formulated as a minimum cost flow problem. Consider a graph
(, ) = GNA which represents a physical network, where N is a set of nodes and A is a
set of links.
[PS] Minimize
(, ) ( , ) (, )
pp
ij ij ij ij
ij ij ij
fycx
∈∈∈
+
∑∑∑
AAP
(48)
Subject to
pp p
ij ji i
jj
x xd
∈∈
−=
∑∑
NN
i ∈ N , p ∈ P (49)
p
ij ij ij
p
x uy
∈
≤
∑
P
(, ) ij ∈ A (50)
0
p
ij
x ≥ (, ) ij ∈ A , p ∈ P (51)
{0,1, 2, }
ij
y ∈ (, ) ij ∈ A (52)
where
ij
f is the fixed cost incurred by opening service on link ( , ) ij ,
ij
y is the discrete
decision variable which represents number of service units (e.g. trucks) offered on link ( , ) ij ,
p
ij
c is the transportation cost per unit flow of product p on link ( , ) ij ,
p
ij
x is the flow
decision variable indicating the amount of flow of commodity p using link ( , ) ij ,
p
i
d is
the demand of product p at node i , and
ij
u is the capacity of link ( , ) ij . The objective
function (48) is to minimize the cost of all relevant decisions, i.e., the service selection and the
traffic distribution. Constraint (49) is the flow conservation and demand satisfaction
requirements. Constraint (50) is the capacity restriction of the flow on each link. Constraints
(51) and (52) ensure the variables are nonnegative and the design variable is an integer.
59
4.2. Service network for the U.S. west cost region
For the regional level, we look into the US west coast region, consisting of multiple
marine ports and its ground transportation network to move goods within the region and in
and out of the region. At this level, we are interested in answering questions such “what is the
impact of appropriate mitigation measures if one of the ports in the region becomes partially
or fully non-functional for a specific period of time?”
109
85
9
8
12
59
60
11
10
84
110
POo
West Coast
SEi
R
SEo
SAo
LAo
POi
SAi
LAi
O
N
Z
N
I
N
(a) FAF’s zoning (b) constructed network
Figure 4.1: Service network for the U.S. west cost region
The existing service network is defined at a high level of aggregation, which includes the
major ports and aggregated terminals (zones) representing broad geographical destinations and
intermediary terminals. To define the regional network at an appropriate level of resolution,
the 2002 Origin-Destination data from the Freight Analysis Framework (FAF) of U.S. DoT
FHWA (Federal Highway Administration) is used [18]. Using the FAF data, the service
60
network of the west coast region, including the state of Nevada, is constructed. This region
consists of 11 zones. In addition, the remainder of the U.S. is considered as one zone and 4
major combined ports (Los Angeles and long beach, San Francisco, Portland, and Seattle) are
also added in the network. Then, the regional service network consists of 20 nodes as shown
in Figure 4.1. The zoning of the FAF is represented in Figure 4.1(a) and the corresponding
network is shown in Figure 4.1(b).
Sine each zone in the FAF data represents an aggregated location of demands, it is
modeled as a node (
Z
N ). A port of disembarkation is not represented as a zone in the data.
However, each port node is added to supply the zones. Also, in order to assign the capacity to
a port node, it is split into an inbound node (
I
N ) and outbound node (
O
N ). In this level of
aggregation, a link does not imply a specific freeway or railway between two corresponding
nodes. Rather, it represents an aggregated path by possible ground transportation modes
(trucks and trains). This distribution network for import freights is assumed to handle an
aggregated single commodity. Then, N is a set of node comprised of port nodes and zone
nodes (
I OZ
=∪ ∪ NN N N ) and A is a set of links comprised of port and ground links
(
PG
=∪ AA A ),
where {(,)| , }
PIO
ij i j =∈ ∈ ANN and {( , )| , , , }
GOZZ
ij i j i j = ∈∈ ≠ ANNN .
If all links are activated to move freight, the regional service network can be modeled as a
minimum cost flow problem.
[PN] Minimize
(, )
ij ij
ij
cx
∈
∑
A
(53)
Subject to
ij ji i
jj
x xd
∈∈
− =
∑∑
NN
i ∈ N (54)
61
ij ij
x u ≤ (, ) ij ∈ A (55)
0
ij
x ≥ (, ) ij ∈ A (56)
where
ij
c is the transportation cost per unit flow in link ( , ) ij ,
ij
x is the flow decision
variable indicating the amount of flow using link ( , ) ij ,
i
d is the demand at node i , and
ij
u is the capacity of link ( , ) ij .
Classical pseudopolynomial-time algorithms, such as, cycle-canceling, primal-dual, out-
of-kilter, and relaxation algorithms can be applied to solve this problem. These algorithms
have some commonalities in that they all repeatedly solve shortest path problem. They
frequently provide the essential building blocks and core ideas used in more efficient
algorithms. Also, it can be solved by a linear programming method. Note that a linear program
with 0-1 incident matrix can be transformed into a minimum cost flow problem [2].
4.3. Mitigating disruptions at regional level
When a disruption occurs at the regional level, the regional service network may need to
be reconfigured. For example, consider the U.S. west coast region presented in Figure 4.1(a).
If the LA port node is rendered non-functional for a period of time, all services originating
from the node will either be discontinued or be operated at lower capacity. The reconfiguration
of the service network is performed by installing a sea transportation mode between ports.
Figure 4.2 shows the possible sea transportation links (dotted arrows) from the LA port node
to other port nodes.
The choice of sea transportation links will determined by an optimization procedure,
which minimizes a generalized cost while satisfying the transportation demands. Our solution
provides the data about the amount of goods that can be distributed by sea links and the
corresponding alternative ground links. Due to capacity constraints of links, some amounts of
62
goods may not be distributed although all the available sea links are activated. The amount of
goods may remains at the disruptive port waiting to be distributed in the next planning horizon.
14
13
6
5
9
10
11
8
7
12
15
3
West Coast
20
16
4
2
1
19
18
17
O
N
Z
N
I
N
Figure 4.2: Reconfiguration of the regional service network
The set of links A are differentiated into port operation links
P
A , ground
transportation links
G
A and sea transportation links
S
A , where
{(,)| , , }
SII
ij i j i j =∈ ∈ ≠ ANN . Then, the service network formulation which can mitigate
disruptions is modeled as a minimum cost flow problem with binary constraints.
[PD] Minimize
(, ) ( , )
S
ij ij ij ij
ij i j
fycx
∈∈
+
∑∑
AA
(57)
Subject to
ij ji i
jj
x xd
∈∈
− =
∑∑
NN
i ∈ N (58)
ij ij ij
x uy ≤ (, )
S
ij ∈ A (59)
ij ij
x u ≤ (, ) ,
P G
ij ∈AA (60)
63
0
ij
x ≥ (, ) ij ∈ A (61)
{0,1 }
ij
y ∈ (, )
S
ij ∈ A (62)
where
ij
f is the activation cost of sea link ( , ) ij ,
ij
y is the discrete decision variable of sea
link (, ) ij ,
ij
c is the transportation cost per unit flow in link ( , ) ij ,
ij
x is the flow decision
variable indicating the amount of flow using link ( , ) ij ,
i
d is the demand at node i , and
ij
u is the capacity of link ( , ) ij .
If the maximum flow saturates all the links from the source, the minimum cost flow
problem is feasible. That is to say, when the supplies/demands at the nodes satisfy the
condition 0
i
i
d
∈
=
∑
N
, the minimum cost flow problem has a feasible solution. To make this
problem feasible by ensuring 0
i
i
d
∈
=
∑
N
, we add a set of slack transportation links between
port-inbound nodes and zone nodes. The set of slack transportation links defined as
{(,)| , }
KIG
ij i j =∈ ∈ ANN . Then, the set of links are redefined as
PG S K
=∪ ∪ ∪ AA A A A . These slack links will provides information how much of goods
cannot be delivered to each demand nodes.
The formulation can also be rewritten by separating the set of slack transportation links.
[PE] Minimize
(, ) ( , ) , , (, )
SPGS K
ij ij ij ij ij ij
ij i j ij
fycx x ω
∈∈ ∈
++
∑∑ ∑
AAAA A
(63)
Subject to
ij ji i
jj
x xd
∈∈
− =
∑∑
NN
i ∈ N (64)
ij ij ij
x uy ≤ (, )
S
ij ∈ A (65)
ij ij
x u ≤ (, ) , ,
P GK
ij ∈AA A (66)
64
0
ij
x ≥ (, ) ij ∈ A (67)
{0,1, 2, }
ij
y ∈ (, )
S
ij ∈ A (68)
where
ij
ω is the penalty to undistributed goods, which could be greater than the cost related
to the alternative sea link choice.
This problem is solved by a branch-and-bound (B&B) method and a linear programming
(LP) relaxation is performed on every leaf of the branch-and-bound tree. The Simplex
algorithm is applied to solve the relaxed LP problem. The Simplex algorithm is defined by the
pivot rule. This rule defines the way that decides with vertex of polyhedron is selected when
there are many basic feasible solutions (BFSs) to choose from [15]. Suppose that the pivot rule
is always to move the adjacent to BFS which there is at least increase in the objective function
value. Under this pivoting rule, the Simplex algorithm requires 21
n
− pivoting steps before
terminating. More precisely, in a standard form, the time complexity of LP is
(2) (2) (2)
nn n
OmnOnO == where m is the number of rows and n is the number of
columns of the incident matrix. In a worst case, since a B&B could generate all leafs on the
tree and it solves a relaxed LP on every leaf, it takes exponential time (( 1)2)
n
a
Op p ⋅− ⋅ to
solve our problem, where
a
p the number of ports in which activated sea links are originated
and p is the number of all ports in the network.
As described above, the Simplex algorithm has exponential time complexity. Its average
behavior and worst case behavior have been studied and explained by Borgwardt [4] and Klee
and Minty [36], respectively. There is no deterministic pivot rule under which the Simplex
algorithm is known to take a sub-exponential number of iterations. However, the numerical
behavior is conflict with theoretical analysis [47]. That means that it is efficient in practice,
while having no polynomial time worst-case complexity, although there are no satisfactory
65
theoretical explanations of its excellent performance. Therefore, the running time in solving
our problem is governed by the term (1)
a
pp ⋅ − .
The B&B method explores the set of feasible integer solutions. However, it uses bounds
on the optimal cost to avoid exploring certain part of the feasible set. The developed algorithm
follows the typical B&B strategies, with few implementation strategies which are listed below.
- It always branches a non-convergent sub-problem unless it is a leaf of the B&B tree.
- To break a sub-problem, it chooses a variable
i
x which is not integer and creates two sub-
problems by adding either of the constraints
*
ii
x x ⎢ ⎥ ≤
⎣ ⎦
or
*
ii
x x ⎢ ⎥ ≥
⎣ ⎦
.
- The lower bound of the optimal cost of a sub-problem is obtained by the linear programming
relaxation using the simplex method.
- The upper bound is updated with the best feasible integer solution so far.
- As a way of choosing an active sub-problem, it uses depth-first search with back-tracking.
4.4. Computational experiments
As we described earlier, our regional service network is constructed based on 2002 FAF’s
zoning and data. The service network of the US west coast region including the state of
Nevada is constructed, which consists of 11 zones for the interested region, 1 zone for the
remainder of the US, and 4 imaginary zones for the major combined port areas (Los Angeles,
San Francisco, Portland, and Seattle). In order to assign the capacity to a port zone, each port
node is split into an inbound node and outbound node. Therefore, the regional service network
consists of 20 nodes as shown in Figure 1b.
The freight demands (in tons) are extracted from multi-dimensional matrices which
consist of origin, destination, commodity, mode, and port of disembarkation. Originally,
origins and destinations consist of 114 regions and 7 international regions. We extracted the
66
freights from the international regions to the interested ports in the west coast regions. Among
these freights, we selected freights whose transportation modes are truck, train, and a
combination of truck and train. Based on this information, the west cost service network is
constructed as shown in Figure 1(b) and the corresponding demands are defined in Table 4.1.
Table 4.1: Import Freight Distribution (ton/day) via Major West Coast Ports
node z8 z9 z10 z11 z12 z59 z60 z84 z85 z109 z110 zR supply
LA 113667 5651 444 2457 4926 1369 61 239 159 757 233 12969142932
SA 385 8 454 105081112 24 59 32 26 76 68 681 13433
PO 321 44 58 34 120 3 25 6643 661 526 1325 935 10695
SE 710 3 31 77 83 1 1 994 67021861 2525 281929775
Table 4.2: Legends and numbers of nodes describing FAF zones
Legend FAF acronym Zone Node
LA CA Los A Los Angeles/Long Beach port areas 17, 1
SA CA San J San Francisco/Oakland port area 18, 2
PO OR Portl Portland port area 19, 3
SE WA Seatt Seattle/Tacoma port areas 20, 4
z8 CA Los A Los Angeles-Long Beach-Riverside, CA 5
z9 CA San D San Diego-Carlsbad-San Marcos, CA 6
z10 CA Sacra Sacramento-Arden-Arcade-Truckee, CA-NV 7
z11 CA San J San Jose-San Francisco-Oakland, CA 8
z12 CA rem Remainder of California 9
z59 NV Las V Las Vegas-Paradise-Pahrump, NV 10
z60 NV rem Remainder of NV 11
z84 OR Portl Portland-Vancouver-Beaverton, OR-WA 12
z85 OR rem Remainder of Oregon 13
z109 WA Seatt Seattle-Tacoma-Olympia, WA 14
z110 WA rem Remainder of Washington 15
zR - Remainder of West Coast Zone 16
It summarizes the supplies and demands of 4 independent distribution networks
originating each combined port complex. First column represents ports of the west coast
region and first row represent zones in the region. Columns 2 to 13 represent demand
XX
i
d of
zone i supplied by port XX, independently. Then,
LA SA PO SE
ii i i i
dd d d d =− − − − . The last
67
column represents the supply to each port. Table 4.2 table describes the legends and indices
used in above table and figures.
The activation cost
ij
f of link (, )
S
ij ∈ A can be defined by ‘the running cost’ of a ship
during the entire trip including the loading at
I
i ∈ N , and the unloading at
I
j ∈ N . The
transportation cost
ij
c of link ( , )
S
ij ∈ A is defined as ‘the price per ton’ since the distance is
already considered in the opening cost. Following relations are used to calculate the opening
cost and transportation cost of sea links, especially from LA to other ports.
Fuel cost (FC) = Loading and Unloading time × Fuel per day at port + Transit time × Fuel per
day at sea
Running Cost = Daily running cost × Total time + Fuel cost + Port cost
Price per ton = (Daily profit × Total times + Running Cost) / Total tons
Besides the distance information between port nodes, several assumptions are made to
calculate these costs in Table 4.3.
- Average ship speed is 25 knot.
Transit time (17-18) = 368/25/24 = 0.6133 hour
- Loading/unloading time is 1 hour.
- Fuel consumption per day at port and at sea is 2 ton and 34 ton, respectively.
- Fuel price per ton is $100.
Fuel cost (17-18) = (1 × 2+0.6 × 34) × 100 = 2240
- Daily running cost of a ship is $5000 and each port cost is $3000.
Activation cost (17-18) = 5000 × 2.6+2240+3000 × 2 = 21240 = f
17,18
- Daily profit required is $1000.
68
- Maximum load of a ship is $10000.
Price per ton (17-18) = (1000 × 2.6+21240)/10000 = 2.384 = c
17,18
Table 4.3: Activation and transportation cost of sea links
Sea link Distance Transit time activation cost Price per ton
i-j (nautical miles) (h) f
ij
c
ij
17-18 368 0.6 21240 2.384
17-19 976 1.7 25480 2.818
17-20 1127 1.9 27160 3.006
According to [45], the price per ton-mile of ground transportation modes (truck and train)
is roughly 22 times higher than that of sea transportation mode. Since our price per-tone mile
of sea links is about $0.004, we can proportionally assumed that the price per ton-mile of our
ground transportation link is about $0.1. Then, the ground transportation cost per unit ton is
0.1
ij ij
cm = , where
ij
m is distance in miles of link ( , )
G
ij ∈ A .
The maximum capacity of link is assumed to be high enough to distribute current
demands found in the FAF data. First, the maximum capacity of port link (, )
P
ij ∈ A is
defined by increasing the current optimal flow which is an optimal solution to [PN].
*
/
m
ij ij p
ux r =, (69)
where
m
ij
u is the maximum capacity,
*
ij
x is the optimal solution to [PN], and
p
r represents
the capacity utilization ratio (0 1
p
r < ≤ ). Without any disruption, a higher
p
r is given to the
LA port node. However, we will generate a disruption in LA node by reducing its
p
r in the
experiments.
Table 4.4 shows the maximum capacity of each link and the residual capacity
r
ij
u which
will be used in PART I of the experiments, later.
69
Table 4.4: Estimating capacity (ton) of port links
Port Link Optimal flow r
p
Maximum
capacity
Residual capacity
17-1 142932 0.9 158813 158813
18-2 13433 0.5 26866 13433
19-3 10695 0.5 21390 10695
20-4 29775 0.6 49625 19850
Secondly, the capacity of ground link ( , )
G
ij ∈ A is defined by increasing the current
optimal flow which is an optimal solution to [PN].
*
/
m
ij ij g
ux r =, (70)
where
m
ij
u is the maximum capacity,
*
ij
x is optimal solution to [PN], and
g
r represents the
capacity utilization ratio. A higher
g
r is given to links whose destination has more demand as
shown in Table 4.5.
Table 4.5: Estimating capacity of ground links
Zone r
g
Links
LA Metropolitan area 0.9 1-5,1-6
Other Metropolitan area 0.6 2-7,2-8,3-12,3-14,4-12,4-14
Other than above 0.5 Other than above
4.4.1. Distributions of freights from LA port node (PART I)
In this part, it is assumed that, if some of freights which cannot be processed in the LA
node, they are distributed through other port nodes. Therefore, higher priority is given to the
supplies of those nodes. That is, the distributions between the other nodes to zones are
performed a priori. That is,
r
ij ij
uu = for ( , )
P
ij ∈ A and
m
ij ij
uu = for ( , )
G
ij ∈ A . Assume
that a disruption occurs in LA port. As a result of the disruption, we assume that the capacity
of link 17-1 is reduced from 90% to 70% in 5% decrements. Table 4.6 shows the performance
70
of our mitigation strategy summarized in terms of the cost of distribution and the amount of
remaining freights.
Table 4.6: Freight distribution from LA port by activating sea links
Capacity Capacity Base Case [PN] Reconfiguration Case [PE]
(%) (ton) RGoods PCost SA PO SE RGoods PCost
90 142932 0 3.53e+6 2962 1388 0 0 3.36e+6
85 134991 7941 1.91e+6 6457 1484 0 0 3.42e+6
80 127050 15882 6.74e+5 6457 4658 4767 0 3.61e+6
75 119110 23822 5.11e+5 6457 8499 6101 2765 3.78e+6
70 111169 31763 4.45e+5 6457 8499 6101 10706 3.72e+6
The first column represents the current capacity of LA port link with respect to the
maximum capacity. The second column shows the corresponding capacity in tonnage. PCost is
the optimal cost of distributed goods supplied to the LA node. RGoods is the amount of
remaining goods at LA node due to capacity insufficiency of the overall service network.
Suppose that the optimal solution to [PD] or [PE] is
**
(, )
ij ij
x y . Then,
**
(, ) ( , ) , ,
PCost
SPGS
ij ij ij ij
ij ij
fycx
∈∈
=+
∑∑
AAAA
.
*
RGoods , ( , )
ij K
xij =∈ A
When the current capacity is down to 90%, all the freights from LA node can be
distributed using only ground links. However, by activating two sea links, the computational
result indicate that the distribution cost originating LA node can be further reduced under our
assumptions regarding costs. If the current capacity of LA node is reduced to less than 90%,
the part of freights supplied to LA node cannot be distributed in [PN]. When the capacity is
reduced to 85% or 80%, the freights supplied to LA zone can still be distributed by opening
71
some of the sea links in [PD]. However, if the capacity of LA node is reduced to 75%, we have
excess freights of 2765 ton even though all the possible sea links from LA are activated.
4.4.2. Distributions of freights from all port nodes (PART II)
In this part, the distributions from all port nodes are performed at the same time and, also,
all the possible sea transportation links between ports would be activated. That is, we have
m
ij ij
uu = for ( , ) ,
PG
ij ∈AA without any disruptions. Again, assume that a disruption occurs
in LA port. As a result of the disruption, we assume that the capacity of link 17-1 is reduced
from 90% to 70% in 5% decrements. Table 4.7 shows the performance of our mitigation
strategy summarized in terms of the cost of distribution and the amount of remaining freights
Table 4.7: Freight distribution from all ports by activating sea links
Base Case [PN] Reconfiguration Case [PE]
Activated links
Capacity
(%)
Capacity
(ton) RGoods PCost
17,18 17,19 17,20 18,19 19,18
RGoods PCost
90 142932 0 4.56e+6 1 0 4.55e+6
85 134991 7941 2.94e+6 1 1 0 4.71e+6
80 127050 15882 1.41e+6 1 1 1 0 4.91e+6
75 119110 23822 1.11e+6 1 1 1 2503 5.29e+6
70 111169 31763 1.05e+6 1 1 1 10444 5.23e+6
The first column represents the current capacity of LA port link with respect to the
maximum capacity. The second column shows the corresponding capacity in tonnage. PCost is
the optimal cost of distributed goods supplied to all ports nodes. RGoods is the amount of
remaining goods at ports due to capacity insufficiency of the overall service network. An
activated sea links is denoted as ‘1’ in the table.
In PART I, the distribution cost of Base Case is 5.01 × 10
6
(3.53 × 10
6
+ costs of other 3
ports). However, if an optimization is performed in this single distribution network, the cost is
reduced to 4.56 × 10
6
(91% of Base Case in PART I). Even though freights can be distributed
using only ground links at 90% capacity level, the distribution cost can be further reduced by
72
activating sea link 19-18. If the current capacity of LA node is reduced to less than 90%, some
of freights supplied to the region cannot be distributed in [PN]. However, if the capacity is
reduced to 85% or 80%, the freights supplied to the region can still be distributed by opening
some of the available sea links in [PD]. However, if the capacity of LA node is reduced to
75%, we have residual freights of 2503 ton mainly due to the excessive supplies to LA port.
Currently, the penalty to the remaining freights is set to high (
ij ij
f ω >> ). But, the related
cost (
*
(, )
K
ij ij
ij
x ω
∈
∑
A
) is not reflected in the table. Note that the related solution
*
ij
x ,
(, )
K
ij ∈ A indicates the unsatisfied demand at each node j, which is waiting for the next
planning horizon. Based on the experimental results, if more goods are required to be
distributed in the current horizon, the increases in the capacities of other ports and of grounds
links should be accompanied by the activation of sea links.
73
Chapter 5 Network link level
5.1. Mitigating disruptions at link level
Toward a future integrated work with the previously described mitigation strategies, we
briefly address the problem of improving the operational characteristics of a link (a freeway
stretch) which is the most vulnerable link in the regional transportation network.
It is found that the most of the freeway congestion is due to poor traffic management
rather than to demand exceeding capacity. Therefore, efficient management and control
supported by advanced communication and control technologies provide potential solutions to
fully utilize the capacity without significant changes to the existing infrastructure. Despite the
local feedback introduced by ramp metering, the current highway system is a dynamical
system that runs almost open loop. In an effort to smooth the speed and density profile along
the lanes, the infrastructure could issue variable speed limits. The variable speed limits
together with ramp metering could serve as the control inputs to close the loop of the
dynamical traffic flow system.
In order to effectively reduce congestion and improve safety on a freeway link, we
consider an integrated freeway traffic control system comprised of ramp metering and variable
speed limit strategies [49],[31]. A computationally feasible and practical controller is
developed for the dynamic traffic system to control flow rates of a ramp and maximum speeds
of a mainline section. The performance of the control system is evaluated both using a
conventional second-order macroscopic traffic flow simulation model and one of popular
microscopic traffic flow simulation models which have been recently used to support ITS
applications.
74
5.2. The macroscopic traffic flow model
The analogy between traffic flow and fluid flow led to the first traffic flow model
proposed by Lighthill and Whitham [38]. This model has only one state variable, which is
traffic density, and shows poor transient behavior. Payne [44] overcame this by adding another
differential equation representing the dynamics of the mean speed. Later, his model is
extended by Cremer and May [41], Papageorgiou et al. [40], and Karaaslan et al. [33]. The
macroscopic model in Papageorgiou et al. [40] which is called as the METANET model is
known to be suitable for describing free flow, critical, and congested traffic condition. This
second-order traffic flow model is adopted for the description and evolution of freeway traffic
flow. A full description of the METANET model can also be found in Kotsialos et al. [37], and
an extended model for speed limits can be referred to Hegyi et al. [23].
5.2.1. Variables and model of the METANET
The freeway network is represented as a directed graph with the links corresponding
freeway stretch. Each stretch has uniform characteristics. The nodes of the graph are placed
where a major change in road geometry occurs, including on-ramp, off-ramp, and junctions.
The macroscopic description of traffic flow implies the definition of variables expressing
the average behavior of traffic at certain times and locations. Consider a link m in Figure 5.1
which is divided into
m
N sections of equal length
m
L . Each segment i of link m is
characterized by following three macroscopic variables: traffic density
,
()
mi
k ρ , mean speed
,
()
mi
vk , and traffic volume or outflow
,
()
mi
qk , where k indicates time instance tkT = and
T is the discrete time step used for the simulation of traffic flow.
75
,1
()
m
k ρ
,1
()
m
vk
,1
()
m
qk
,
()
mi
k ρ
,
()
mi
vk
,
()
mi
qk
,
()
mN
k ρ
,
()
mN
vk
,
()
mN
qk
m
L
Figure 5.1: Uni-directional freeway stretch divided into N sections
Then, the traffic variables for each segment i of link m at each time step k are
calculated by the following equations. However, for the sake of notational simplicity, the link
index m is omitted hereafter.
() () ()
ii i
qk kv k ρ λ = (71)
()
1
(1) () () ()
ii i i
T
kk qkqk
L
ρρ
λ
−
+= + − (72)
() ()
1
1
( 1) () ( ())() () ()()
() ()
()
ii i i i i i
ii
i
TT
vk vk V k vk vk v k vk
L
kk T
Lk
ρ
τ
ρρ υ
τρ κ
−
+
+= + − + −
−
−
+
(73)
free
crit
() 1
(()) exp
m
a
i
i
m
k
Vk v
a
ρ
ρ
ρ
⎡⎤
⎛⎞
⎢⎥ =⋅ −
⎜⎟
⎢⎥
⎝⎠
⎣⎦
(74)
where λ denotes the number of lanes, τ , υ , κ , and
m
a are model parameters, the free
flow speed
free
v is the average speed that drivers assume if traffic is freely moving, and the
critical density
crit
ρ is the density at which the traffic flow is maximal.
The density equation (72) describes the relation between the density at time ( 1) kT + and
kT by considering the vehicles entering and leaving section i . The traffic volume leaving
section i (71) is equal to the density multiplied by the mean speed and the number of lanes
on that section. The mean speed at time ( 1) kT + (73) is described by the mean speed at time
kT , the relaxation term, convection term, and the anticipation term. The relaxation term
76
describes the convergence of the mean speed to its equilibrium values at a rate determined by
the time constant τ . The convection term expresses the speed changes by the incoming
vehicles which have a different speed than the mean speed of that section. Finally, the
anticipation term describes the driver response to the downstream density. If the downstream
density is lower, drivers tend to speed up and vice versa. The parameter κ insures that the
last term is less influence in low density situation.
Also, we use different values for the sensitivity parameter υ of the anticipation term in
(73) as appeared in Hegyi et al. [24]. This was introduced to express the different reactions of
drivers to varying downstream densities.
11
21
if () ()
if () ()
ii
ii
kk
kk
υρ ρ
υ
υρ ρ
+
+
≥ ⎧
=
⎨
<
⎩
(75)
Origin links which receive the traffic demand are modeled with a simple queue model.
() (1) () () ()
oo o o
wk wk T d k q k += + − (76)
where ()
o
dk is the traffic demand and ( )
o
qk is the outflow.
The outflow ()
o
qk is determined by either the current origin demand, the ramp capacity,
or the reduced capacity based on the mainstream condition.
max 1
0
max crit
() ()
()min () , (),
o
oo oo
wk k
qk d k Q r k Q
T
ρρ
ρρ
⎧⎫ −
=+ ⋅
⎨⎬
−
⎩⎭
(77)
where ()
o
Qk is the on-ramp capacity under free-flow condition, ( )
o
rk is the metering rate
(( ) 1
o
rk < ), and
max
ρ is the maximum density of the link to which the on-ramp is connected.
If there is an on-ramp, then, the negative is added to (73) to account for the speed drop
caused by the on-ramp merging.
77
1
1
() ()
(() )
o
Tq k v k
Lk
δ
λρκ
−
+
(78)
where δ is a model parameter.
5.2.2. Extensions of the METANET model for speed limits
Since the METANET model does not describe the effect of speed limits, it is extended so
that speed limits are incorporated into the mean speed equation [23]. It is assumed that the
desired speed is determined by either the desired speed based on the experience density or the
desired speed caused by the speed limit displayed.
free ctrl,
crit
() 1
(( )) min exp , (1 ) ( )
m
a
i
i i
m
k
Vk v v k
a
ρ
ρα
ρ
⎧⎫ ⎡⎤
⎛⎞
⎪⎪
⎢⎥ =⋅ − +
⎨⎬ ⎜⎟
⎢⎥
⎝⎠
⎪⎪
⎣⎦ ⎩⎭
(79)
where
ctrl,
()
i
vk is the speed limit imposed on segment i and α is introduced to express the
non-compliance factor since drivers target speed is usually higher than what is displayed.
A mainline origin link o is modified from (77) because the outflow can be limited by an
active speed limit or by the actual speed on the first segment.
lim,1
()
() min () , ()
o
oo
wk
qk d k q k
T
⎧⎫
=+
⎨⎬
⎩⎭
(80)
where
lim,1
() qk is the maximum outflow determined by the limiting speed on the first
segment which is
lim,1 lim,1
( ) min{ (), ()}
k
vk v k vk = .
1
lim,1
lim,1 crit lim,1 crit
lim,1 free
crit crit lim,1 crit
()
() ln if () ( )
()
() if () ()
a
vk
vk a v k V
qk v
VvkV
λρ ρ
λρ ρ ρ
⎧
⎡⎤ ⎛⎞
⎪
⋅⋅ − < ⎪
⎢⎥ ⎜⎟
=
⎨
⎝⎠ ⎣⎦
⎪
⋅⋅ ≥
⎪
⎩
(81)
78
It is assumed that the boundary conditions for the upstream speed equals to the speed of
the first segment. It is also assumed that the density of the leaving link equals to the density of
the last segment in free flow or equals to the critical density in congested flow.
01
() () vk v k = (82)
crit
1
crit crit
() if ()
()
if ()
NN
N
N
kk
k
k
ρ ρρ
ρ
ρ ρρ
+
< ⎧
=
⎨
≥
⎩
(83)
5.3. Model predictive control
The problem of the coordinated ramp metering and speed limits is formulated as a
discrete-time non-linear optimal control problem with constrained control variables which can
be solved numerically over a given-time horizon. For the optimal control problem, the state
vector consists of the density
i
ρ and mean speed
i
v of every segment i of the link and the
queue
o
w of every origin. The control vector consists of the ramp metering rate
o
r of every
on-ramp and speed limit
lim,i
v of every segment i of the link. The disturbance vector
consists of the demand
o
d at every origin o .
To solve the problem, we apply the model predictive control scheme in a rolling horizon
framework proposed in Hegyi et al. [23],[24]. This approach is chosen for comparing the
performance of our proposed approaches which will be introduced later. The chosen objective
function contains the total time spent (TTS) which is a combination of the total travel time
(TTT) and the total waiting time (TWT) at origins. The second term penalizes abrupt
variations in the ramp metering and speed limit control signals. It is assumed that the
controller time step is an integer multiple M of the simulation time step.
79
{}
()1
2
1
ctrl, ctrl, 2
ramp ctrl
free
( ) () ()
() ( 1)
(() ( 1))
p
c
Mk N
io
io
jMk
kN
ii
oo
lk o i
Jk T j L w j
vl vl
arlrl a
v
ρλ
+−
=
+−
=
=+
⎧⎫
−− ⎛⎞ ⎪⎪
+−−+
⎨⎬
⎜⎟
⎝⎠
⎪⎪
⎩⎭
∑∑ ∑
∑∑ ∑
(84)
where
ramp
a and
ctrl
a are non-negative weighting parameters for ramp and speed control
signals,
p
N is a prediction horizon, and
c
N is a control horizon.
The utilized numerical algorithm to solve the optimization problem is the MATLAB
implementation of the SQP (sequential quadratic programming) algorithm, fmincon.
5.4. ALINEA-base ramp metering and speed limit control strategy
ALINEA is a simple local ramp metering strategy based on pure integral control action.
Different studies have demonstrated that ALINEA is non-inferior to sophisticated coordinated
approaches under recurrent traffic congestion [41].
This strategy developed in [31] is referred to as the DEnsity based RAmp Metering
(DERAM) strategy. The ramp metering strategy is described as follows:
()
( )
()
()
max 1 max
1min 1 min
1
, if
, if
, otherwise
j
jj
j
RRnTR
RnT R RnT R
RnT
⎧ >
⎪
⎪
=<
⎨
⎪
⎪
⎩
(85)
where () ( ) () ()()
11 1 00
1
1
1+ 1
c
N
jj rd j c
m
c
R nT R n T K n N T mT
N
ρρ
=
⎛⎞
=− − − +
⎜⎟
⎝⎠
∑
and
r
K is a
positive controller parameter.
For the mainline speed control strategy, section i is regarded as a virtual on ramp of
section 1 i + and the same DERAM strategy as in the case of ramp metering is applied to
regulate the flow rate
i
Q from section i to section 1 i + , i.e.,
80
()
( )
()
()
max 1 max
1 min 1 min
1
, if >
, if
, otherwise,
i
ii
i
QQnTQ
QnT Q QnT Q
QnT
⎧
⎪
⎪
=<
⎨
⎪
⎪
⎩
(86)
where
() ( ) () ()()
11 2 1 00
1
1
1+ 1
c
N
ii vd i c
m
c
QnT Q n T K n NT mT
N
ρρ
+
=
⎛⎞
=− − − +
⎜⎟
⎝⎠
∑
and
v
K is a
controller parameter.
The above equations provide the regulation of the flow at a particular section of the
highway. This mainline traffic flow control strategy cannot be implemented as done in the
case of ramp metering, because in this case the control variable is the speed. Therefore, in
order to regulate traffic speed instead of the traffic flow rate, we map a flow command into a
speed command using the flow-speed relationship as described by the fundamental flow-
density diagram shown in Figure 5.2. It is only reasonable to set upper and lower bounds on
the speed commands of section i ,
1
()
i
VnT , where
min 1 max
()
i
VVnT V ≤ ≤ ,
max
V is the
maximum speed allowed, which often set to be the traditional speed limit, and
min
V is the
lowest dynamic speed limit we want to apply. Therefore, we set
min
Q as the flow
corresponding to
min
V . Because the capacity flow is usually not achieved at the maximum
speed allowed, we set
max
Q as the flow corresponding to the critical density which is the
capacity flow. Denote the speed corresponding to the critical density as
crit
V and we have
min max crit
VV V ≤≤ . The curve between point A and point B is a mapping from
min max
[, ] QQ to
min
[, ]
crit
VV . In this paper, we simplify this curve to a straight line AB, and denoted this
mapping as ( ) f Q . The point A and B on the flow-density relationship of every section could
be estimated either offline or online, therefore, ( ) f Q could be estimated.
81
max
Q
min
Q
critical
V
min
V
c
ρ
max
V
Figure 5.2: Fundamental flow-density diagram
Combined DERAM strategy with the mapping ( ) f Q , we have the desired speed limit
for section i,
( )
11
(( ))
ii
VnT f Q nT = (87)
However,
i
V generated by (87) may lead to unsafe changes of speed limits. For practical
purposes, we use the following speed limit
i
V which is smoother.
()
() ( )
() ( ) ()
()
() ()
()
1
11
111
11 1
1
1,
if 1
,
if
, otherwise,
iv
ii v
ii v
ii v
i
Vn T c
VnT V n T c
VnT V nT c
VnT V nT c
VnT
+
+
⎧ −−
⎪
≤ −− ⎪
⎪
⎪
=+
⎨
⎪
≥+
⎪
⎪
⎪
⎩
(88)
where
v
c is a design constant.
5.5. Simulation experiments
5.5.1. Macroscopic simulation
A benchmark network for evaluating both MPC and PI control methodologies is
considered. The network consists of two origins: a mainline and an on-ramp. The freeway
82
mainline has five lanes, which is a 5km long link consisting of 10 segments of 500m each. A
single-lane onramp is attached to section 9 after a 4 km long undisturbed link. The sections 5
to 8 are equipped with VMS (variable message signs), which are the candidate sections where
speed limits are applied. This network is carefully prepared to make our simulation setup
similar to one shown in Hegyi et al. [23]. Subsequently, we use the model parameters as found
in Hegyi et al. [23],[24].
For the MPC methodology, the prediction horizon 42
p
N = is used to reflect the typical
travel time through the 4 km link. The control horizon for the ramp metering only case and the
coordinated ramp metering and speed limit case are set to be 3
c
N = and 5
c
N = ,
respectively. The on-ramp queue is set to 100 veh/h. The other parameters are
500 m, L = 10 s, T = 18 s, τ = 40 veh/km/lane, κ =
2
2
30 km h, υ =
max
180 veh/km/lane, ρ = 0.0122, δ =
ramp ctrl
0.4, aa = =
crit
33.5 veh/km/lane, ρ = 1.867, a =
free
102 km/h, v = and 0.1 α = .
To examine the effect of the ramp metering and the variable speed limits, a simple, but
typical, demand scenario is considered. The mainline demand has a constant value of 1800
veh/km/lane for 1.5 h close to the capacity (2000 veh/km/lane) while the on-ramp demand
reaches to 1500 veh/km/lane for 1 h. This setup is prepared to reflect the on-ramp burst during
peak hours. Both demands later decreases to 500 veh/km/lane so that the final states are the
same for all cases considered. This enables us to properly compare the TTS of the different
cases. Under the given demand scenario, the following five different cases are compared.
Case 1: No control
Case 2: MPC ramp metering only
Case 3: MPC ramp metering and speed limits
83
Case 4: ALINEA-based ramp metering only
Case 5: ALINEA-based ramp metering and speed limits
Table 5.1 summarizes the performances of the control strategies in terms of their TTS.
The simulation results for no control case (Case 1) shows that, as demands increases on both
the mainline and the on-ramp, the congestion is created and propagates through the mainline.
This congestion yields a queue outside the network of slightly more than 150 veh. The TTS of
the case 1 is 1892.3 veh/h. The MPC ramp metering effectively suppresses the onset of the
congestion before the on-ramp queue reaches to the maximum. It eventually reduces the
magnitude and the duration of the later congestion and yields the 7.2% improvement (1755.3
veh/h). However, from the point that the on-ramp queue reaches to the maximum, the speed
limits is more actively applied to reduce the inflow of the near upstream sections. The TTS of
the case 3 is 1629.7 veh/h, which is an improvement of 13.9% from Case 1. This improvement
in terms of the TTS coincides with the results shown in Hegyi et al. [23]. Therefore, this could
be considered as the upper bound of the improvement by Cases 4 and 5 which are ultimately
designed for the microscopic (or real-time) traffic simulation environments.
Table 5.1: Performances under macroscopic environments
Case Description TTS Reduction
1 No control 1892.3 -
2 MPC ramp metering only 1755.3 7.2 %
3 MPC ramp and speed limits 1629.7 13.9 %
4 ALINEA ramp metering only 1763.7 6.8 %
5 ALINEA ramp metering and speed limits 1711.2 10.0 %
The performance of the ALINEA ramp metering only (Case 4) is almost close to the case
2. Due to the lack of the prediction of the traffic evolution and the lack of information about
the demand flows, the improvement in Case 5 is less, but slightly, than Case 3. The remaining
congestion outside the controlled sections is little higher than that of Case 3. This causes a
84
queue at the mainline of approximately 50 vehicles for a while. However, it substantially
improves the network performance by suppressing the onset and the propagation of the
congestion.
5.5.2. Microscopic simulation
Microscopic traffic simulations have become common tools to support Intelligent
Transportation Systems (ITS) applications. To design, analyze, and evaluate various needs,
these tools should meet a wide range of requirements. These include signal control design,
variable speed limit design, route planning, noise and emission assessment, etc. Sophisticated
tools, such as PARAMICS, AIMSUN, and VISSIM, are commonly used to model any
combination of surface transportation networks at a high level of detail. They support most
signal control and other operational strategies and provide various output formats for analysis
and comparison.
Due to its flexibility, VISSIM is selected to simulate various needs in the field of traffic
analysis. VISSIM is a discrete, stochastic, time step based microscopic traffic simulation
program developed to analyze traffic and transit operations. VISSIM uses the psycho-physical
driver behavior model based on the work of Wiedemann [51]. The basic idea of this model is
stochastic perceptual thresholds which replicate individual driver characteristics. Detailed
configuration of VISSIM for setting for our purposes is available in our previous work
[31],[49].
Our simulations in the microscopic traffic flow simulation tool require on-line data
collection and on-line command generation. Therefore, a simulation framework is developed
to access the VISSIM objects via the COM interface and to communicate with the proposed
controllers via the MATLAB external interface.
85
Two congestion scenarios were created on the benchmark network constructed in VISSIM.
In scenario 1, the demanding flow of the on ramp at section 9 is set to a high value (1100
veh/h/lane) from 0.5 hour to 1.5 hour during the 2.5 hour simulation run. Therefore, the
disturbance from the on ramp creates a bottleneck at section 9. In scenario 2, a vehicle stops in
lane 3 at section 9 from 800 sec to 1200 sec which creates a lane blockage for about 6mins. To
quantify the controller performance, Total Time Spent (TTS) in the network is used to assess
the average travel time, and standard deviation of density and speed are used to assess the
smoothness of traffic. Environmental effects and safety effects are related to the standard
deviation of density since the smoother the density of the segment is, the fewer acceleration or
deceleration events take place. Therefore, smaller density deviation is an indicator of lower
emission rates and lower possibility of accidents.
Table 5.2 summarizes the performance of the proposed control strategies in microscopic
simulation environments For scenario 1, if both DERAM ramp metering and speed limits
control are applied, the TTS is reduced by 5.5%, which is more than just applying the ramp
metering. When compared with the macroscopic simulation results, the TTS reduction is less
which is consistent with the difference between the two models. Moreover, the standard
deviation of density and average speed has been reduced by 22.5% and 17%, respectively.
These standard deviations are indicators of traffic smoothness and, therefore, indicators of
better safety and less pollution.
For scenario 2, ramp metering, speed limit control and lane change guidance are all active
after the accident. Lane change guidance starts at 2 sections upstream of the accident. During
the 2nd upstream section, 20% of the vehicles in the accident lane are assumed to change to
their left lanes and another 20% are assumed to change to their right lanes. During the
immediately upstream section, the third 20% of the vehicles in the accident lane are assumed
86
to change to their left lanes and the fourth 20% change to the right. The left 20% of the
vehicles in the accident lane are assumed to change lanes during the last 200 to 300 meters
before the accident spot. Due to the facts that the accident is unpredictable and one lane
blockage physically reduces the roadway capacity by at least 1/5, the TTS has not been
reduced much. The smoothing effect of the controller is showing in term of standard
deviations of density and speeds.
Table 5.2: Performance of the proposed DERAM strategy
DERAM
Scenario Measure No control
ramp only Integrated
TTS (-%) 1721.3 (0) 1650.3 (4.1) 1626.4 (5.5)
1: on ramp burst Std. of d (-%) 13.0 (0) 10.6 (18.2) 10.1 (22.5)
Std. of v (-%) 25.2 (0) 22.1 (12.2) 20.9 (17.0)
TTS (-%) 823.0 (0) - 806.8 (2.0)
2: disturbance Std. of d (-%) 5.8 (0) - 4.5 (22.9)
Std. of v (-%) 14.8 (0) - 13.1 (11.5)
5.6. Integrated simulations for future works
Higher efficiency within the terminals/ports, however, could imply more truck traffic in
and out of the terminal that will lead to congestion on the roadway network adjacent to the
terminals. This traffic congestion will, in turn have a negative effect on the efficiency of the
terminals, in addition to increased pollution. It is therefore clear that the efficiency and
capacity of the terminals cannot be decoupled from the effect on congestion on the traffic
network outside the terminals. That is, the problem of maintaining high terminal throughput
while at the same time managing congestion and maintaining traffic efficiency in the traffic
network outside the terminals should be viewed as an integrated problem.
For this reason, authors develop a simulation test bed to evaluate the various concepts
with respect to the performance and cost of terminals as well as the impact on the traffic
network [7]. It consists of macroscopic terminal operation model (TermSim), microscopic
87
traffic simulation model (TrafficSim), and cost analysis model (TermCost). The TermSim
module allows the simulation of container movements inside the terminal by modeling
terminal operations. The TrafficSim module enables the microscopic simulation of the traffic
flow on the traffic road network outside the terminals. The TermCost module enables the user
to perform a cost analysis of the impact on the terminal of any changes in the terminal
characteristics, volume of containers, new approaches etc.
In [6, 8, and 32], the authors studied the problem of reducing truck traffic by using
optimization and information technologies to develop methods of empty container reuse that
will require a smaller number of truck trips to and from the terminal. The analytical models
and optimization techniques show that empty container reuse can produce significant
reduction in truck miles traveled, which in turn reduces traffic congestion, pollution and the
associated negative impact on the surrounding communities. The simulation framework was
implemented to evaluate the impact of different concepts including the empty container reuse,
which have an effect on truck movements and volume of trucks on the terminals and adjacent
traffic network [7, 9].
Eventually, more decision making will be possible to increase the capacity utilization and
to mitigate the impact of disruptions if distribution operations is combined with scheduling
operations in the major nodes (terminals and ports) and the traffic control operations in major
transportation links (freeways). Toward a future integrated work with proposed mitigation
strategies, our ongoing works have been focused on developing more disaggregated
simulation models and on integrating more of such models. Due to the transition to
microscopic simulation models and the use of object-orient programming, a future simulation
framework for freight distribution systems will have the software architecture as showed in
Figure 5.3.
88
traffic simulation
model
coordinated
traffic controllers
client/operator
signal control
logic
run command
control
commands
terminal (port)
simulator
cost
analyzers
run command
traffic
information
resource information
traffic
information
traffic
information
control
commands
control
commands
resource allocation
schedulers
traffic
information
allocation
schedules
routing
schedules
advanced truck
movement
concepts
traffic simulator
Figure 5.3: Architecture for an integrated simulation frameworks
89
Conclusions
In this work, we investigate the ways to effectively utilize resources in a regional freight
transportation network. By improving the operational characteristics of the major resources in
the network, we showed that the underutilization of capacities can be reduced and the impact
of disruptions can be mitigated. To address strategies for reducing the impact of disruptions at
different levels i.e., terminal, port and regional levels, we studied methods of modeling and
optimizing the scheduling operation in terminal and port level and the distribution operation in
a regional transportation network.
In the terminal level, the berth allocation problem (BAP) is addressed as the essential
allocation problem in a container terminal. We developed mathematical models and
optimization methods based on a subgradient optimization and simulated annealing methods.
The proposed methods consist of a set of systematic heuristics to find an initial feasible
solution and update a current solution by exploring a sufficiently large solution space.
Computational experiments showed that the proposed heuristics can deal with disruptions
caused by delays by effectively finding the best allocation such that total delay is minimized.
In the port level, the multiple berth allocation problem (mBAP) is addressed more demanding
situations. When disruptions render a berth unable to meet the demands, or even make it
partially or totally non-functional, other berths within the terminal or within a port can be
called upon to assist in handling the excessive work. The mBAP was modeled as a set
partitioning problem and a heuristic procedure is developed to partition the set of vessels into
a finite number of berths. Computation experiments showed that the impact of disruptive
events can be significantly reduced under a wide range of simulation scenarios. In the regional
level, we focused on the service network in the US west coast region, consisting of multiple
ports and the associated transportation network. The regional service network is defined at a
90
high level of aggregation and is modeled as a minimum cost flow problem. Under a disruptive
event, the network is reconfigured by opening sea transportation mode between ports. This
problem is modeled as a minimum cost flow problem with binary (or, integer) constraints.
Computational experiments showed that the possible cost reduction acquired under the
maximum utilization of the network capacity.
Toward a future integrated work with proposed mitigation strategies, we addressed the
problem of improving the capacity of a freeway stretch. Specifically, to reduce congestion and
improve safety, we implemented an integrated freeway control system comprised of ramp
metering and variable speed limits. Eventually, more decision making will be possible to
increase the capacity utilization and to mitigate the impact of disruptions if distribution
operations is combined with scheduling operations in the major nodes (terminals and ports)
and the traffic control operations in major transportation links (freeways). Using more models,
to be precise, more disaggregated models, our ongoing works have been focused on
implementing integrated simulation frameworks which can be comprised of terminal
simulation models, terminal cost model, and traffic simulation models.
91
References
[1] AAPA (American Association of Port Authorities), 2009.
http://aapa.files.cms-plus.com/PDFs/North_American_Container_Traffic.pdf
[2] R. K. Ahuja, T. L. Magnanti, and J. B. Orlin, Network flows: Theory, algorithms, and applications,
Prentice-Hall, New Jersey, 2003.
[3] E. Balas and M. C. Carrera, “A Dynamic Subgradient-based Branch-and-Bound Procedure for Set
Covering,” Operations Research, vol. 44, no. 6, pp. 875-890, 1996.
[4] K. H. Borgwardt, The Simplex Method: A Probabilistic Analysis, Berlin: Springer–Verlag,
Heidelberg, 1987.
[5] V. Cerny, “A thermodynamical approach to the travelling salesman problem: an efficient
simulation algorithm,” Journal of Optimization Theory and Applications, vol. 45, pp. 41-51, 1985.
[6] H. Chang, H. Jula, A. Chassiakos, and P. Ioannou, “Empty Container Reuse in the Los
Angeles/Long Beach Port Area,” 1st Annual National Urban Freight Conference 2006, Long
Beach, CA, Feb. 2006.
[7] H. Chang, G. Valencia, A. Chassiakos, and P. Ioannou, “Simulation Test-bed and Evaluation of
Truck Movement Concepts on Terminal Efficiency and Traffic Flow,” METRANS report, 2007.
[8] H. Chang, H. Jula, A. Chassiakos, and P. Ioannou, “A heuristic solution for the empty container
substitution problem,” Transportation Research Part E, vol. 44, no. 2, pp. 203-216, 2008.
[9] H. Chang, H. Jula, A. Chassiakos, and P. Ioannou, “Integrated simulation of freeway traffic flow
and container terminal operation: framework and case study,” IFAC Symposium on Control in
Transportation Systems, Redondo Beach, CA, Sept. 2009.
[10] H. Choset, “Coverage for robotics – A survey of recent results,” Annals of Mathematics and
Artificial Intelligence, vol. 31, pp. 113-126, 2001.
[11] J.-F. Cordeau, G. Laporte, P. Legato, and L. Moccia, “Models and Tabu Search Heuristics for the
Berth-Allocation Problem,” Transportation Science, vol. 30, no. 4, pp. 526-538, Nov. 2005.
[12] D. Y. Coulter, “Globalization of Maritime Commerce: The Rise of Hub Ports,” in Globalization
and Maritime Power, Edited by Sam J Tangredi, Institute for National Strategic Studies, Dec.
2002.
[13] T. G. Crainic, “Service network design in freight transportation,” European Journal of Operations
Research, vol. 122, pp. 272-288, 2000.
[14] M. Cremer and A. D. May, “An extended traffic model for freeway control,” California PATH
Research Report UCB-ITS-RR-85-7, 1985.
[15] G. B. Dantzig, Linear Programming and Extensions, Princeton University Press, Princeton, 1963.
[16] M. L. Fisher, “The Lagrangian relaxation method for solving integer programming problems,”
Management Science, vol. 27, no. 1, pp. 1-18, 1981.
92
[17] M. L. Fisher and L. A. Wolsey, “On the Greedy Heuristic for Covering and Packing Problems,”
SIAM Journal on Algebraic and Discrete Methods, vol. 3, no. 4, pp. 584-591, 1982.
[18] Freight Analysis Framework (FAF), FHWA, U.S. DoT, 2002.
http://ops.fhwa.dot.gov/freight/freight_analysis/faf/index.htm
[19] M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-
Completeness, Freeman, San Francisco, CA, 1979.
[20] A. M. Geoffrion, “Lagrangian relaxation for and its uses in integer programming,” Mathematical
Programming Study, vol. 2, pp. 82-114, 1974.
[21] M. Golias, S. Theofanis S., and M. Boile, “Berth and quay crane scheduling: a formulation
reflecting start and finish service deadlines and productivity agreement.” Presented at the 2nd
Annual National Urban Freight Conference, Long Beach, CA, December 2007.
[22] Y. Guan and R. Cheung, “The berth allocation problem: models and solution methods,” OR
Spectrum, vol. 26, no. 1, pp. 75-92, 2004.
[23] A. Hegyi, B. D. Schutter, and H. Hellendoorn, “Model predictive control for optimal coordination
of ramp metering and variable speed limits,” Transportation Research Part C, vol. 13, pp. 185-209,
2005.
[24] A. Hegyi, B. D. Schutter, and H. Hellendoorm, “Optimal coordination of variable speed limits to
suppress shock waves,” IEEE ITS, vol. 6, no. 1, pp. 102-112, 2005.
[25] M. Held, P. Wolfe, and H. P. Crowder, “Validation of subgradient optimization, Mathematical
Programming,” vol. 6, pp. 62-88, 1974.
[26] K. L. Hoffman and M. Padberg, “Solving Airline Crew Scheduling Problems by Branch and
Cut,” Management Science, vol. 39, no. 6, pp. 657-682, 1993.
[27] A. Imai, E. Nishimura, and S. Papadimitriou, “The dynamic berth allocation problem for a
container port,” Transportation Research Part B, vol. 35, no. 4, pp. 401-417, 2001.
[28] A. Imai, E. Nishimura, and S. Papadimitriou, “Berth allocation with service priority,”
Transportation Research Part B, vol. 37, no. 5, pp. 437-457, 2003.
[29] A. Imai, X. Sun, E. Nishimura, and S. Papadimitriou, “Berth allocation in a container port: using
a continuous location space approach,” Transportation Research Part B, vol. 39, no. 3, pp. 199-
221, 2005.
[30] A. Imai, E. Nishimura, M. Hattori, and S. Papadimitriou, “Berth allocation at indented berths for
mega-containerships,” European Journal of Operational Research, 179, pp. 579–593, 2007.
[31] P. Ioannou, H. Chang, and Y. Wang, “Control of Highway Traffic Flow: Design and Evaluation,”
submitted to IEEE ITS, January 2008.
[32] H. Jula, A. Chassiakos, and P. Ioannou, “Port dynamic empty container reuse,” Transportation
Research Part E, vol. 42, no. 1, pp. 43-60, 2006.
93
[33] U. Karaaslan, P. Variya, and J. Walrand, “Two proposals to improve freeway traffic flow”,
California PATH research report UCB-ITS-PRR-90-6, 1990.
[34] K. H. Kim and K. C. Moon, “Berth scheduling by simulated annealing,” Transportation Research
Part B, vol. 37, no. 6, pp. 541-560, 2003.
[35] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by Simulated Annealing,” Science,
New Series, vol. 23, no. 4598, pp. 671-680, 1983.
[36] V. Klee and G. J. Minty, “How good is the Simplex algorithm?” in Inequalities III (Edited by O.
Shisha), Academic Press, New York, pp. 159-175, 1972.
[37] A. Kotsialos, M. Papageorgiou, C. Diakaki, Y. Pavlis, and F. Middelham, “Traffic flow modeling
of large-scale motorway networks using the macroscopic modeling tool METANET,” IEEE ITS,
vol. 3, no. 4, pp. 282-292, 2002.
[38] M. J. Lighthill and G. B. Whitham, “On kinematic waves ii: A theory of traffic flow on long
crowed roads,” Proc. R. Soc. Lond. A229, pp. 317-345, 1995.
[39] R. Moorthy and C.-P. Teo, “Berth management in container terminal: the template design
problem,” OR Spectrum, vol. 28, no. 4, pp. 485-518, 2006.
[40] M. Papageorgiou, J. M. Blosseville and H. Hadj-Salem, “Modeling and real-time control of
traffic flow on the southern part of Boulevard Peripherique in Paris: Part I: Modelling,”
Transportation Research Part A, vol.24A, no. 5, pp. 361-370, 1990.
[41] M. Papageorgiou, H. Hadj-Salem, and F. Middelham, “ALINEA local ramp metering: summary
of field results,” Transportation Research Record, no. 1603, pp. 90-98, 1997.
[42] K. T. Park and K. H. Kim, “Berth scheduling for container terminals by using a sub-gradient
optimization technique,” Journal of the Operational Research Society, vol. 53, no. 9, pp. 1054-
1062, 2002.
[43] Y . M. Park and K. H. Kim, “A scheduling method for Berth and Quay cranes,” OR Spectrum, vol.
25, pp. 1-23, 2003.
[44] H. J. Payne, “Models of freeway traffic and control,” Mathematical Models of Public Systems,
Simulation Council Proceedings, pp. 51-61, 1971.
[45] RITA, Bureau of Transportation Statistics, 2009.
https://www.bts.gov/publications/national_transportation_statistics/index.html
[46] N. K. Ryan, “The Future of Maritime Facility Designs and Operations,” Proceedings of the 1998
Winter Simulation Conference, pp. 1223-1227, 1998.
[47] R. Shamir, “The efficiency of the simplex method: A survey,” Management Science, vol. 33, no. 3,
pp. 301-334, 1987.
[48] M. J. Vickerman, “Planning for Future Transportation Technological Change,” 2005 Executive
Management Conference, May 2005.
94
[49] Y. Wang, H. Chang, J. Zhang, and P. Ioannou, “Integrated roadway/adaptive cruise control
system: safety, performance, environmental and near term deployment considerations,” California
PATH Research Report, UC Berkeley PATH, 2007.
[50] N. Wieberneit, “Service network design for freight transportation: a review,” OR Spectrum, vol.
30, no. 1, pp. 77-112, 2007.
[51] R. Wiedemann, “Modeling of RTI-Elements on multi-lane roads,” in Advanced Tematics in Road
Transport edited by the Comission of the European Community, XIII, Brussels, 1991.
Abstract (if available)
Abstract
Currently, due to the continuing increases in demand, the U.S. freight transportation system is running at capacity and disruptions are commonplace. However, by improving the operational characteristics of the major resources in the network, the underutilization of capacities can be reduced and, therefore, the impact of disruptions can be mitigated. The purpose of this research is to study methods of modeling and evaluating the freight operations and distributions and to develop mitigation strategies for reducing their impact at terminal, port, and regional levels with special emphasis given to the ports of the U.S. west coast region.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Optimum multimodal routing under normal condition and disruptions
PDF
Models and algorithms for the freight movement problem in drayage operations
PDF
Intelligent urban freight transportation
PDF
Strategies for effective rail track capacity use
PDF
Control of mainstream traffic flow: variable speed limit and lane change
PDF
Disruptive innovation’s impact on the staffing and recruitment industry: a case for capacity building for future technological adaptation
PDF
Staying ahead of the digital tsunami: strategy, innovation and change in public media organizations
PDF
Material and process development and optimization for efficient manufacturing of polymer composites
PDF
Integration of truck scheduling and routing with parking availability
PDF
Intelligent near-optimal resource allocation and sharing for self-reconfigurable robotic and other networks
PDF
Developing frameworks to quantify the operational and environmental performance of energy systems within the context of climate change
PDF
Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures
PDF
Measuring the impact of CDN design decisions
PDF
Beyond greenhouse gases and towards urban-scale climate mitigation: understanding the roles of black carbon aerosols and the urban heat island effect as local to regional radiative forcing agents
PDF
Quiet leading from the depths of pharmaceutical organizations: quality IT leader strategies for nurturing diverse project teams
PDF
Investigating the role of climate in affecting residential electricity consumption through high spatiotemporal resolution observations
PDF
The roles of surface and pore properties in wetting resistance for membrane distillation membranes
PDF
Negotiation strategies for women impacting the expanding gender earnings gap at midlife
PDF
Impact of exposure to brine/CO₂ on the mechanical and transport properties of the Mt. Simon sandstone
PDF
Characterization of black carbon: from source to evolution of physical and optical properties in the atmosphere
Asset Metadata
Creator
Chang, Hwan
(author)
Core Title
Reconfiguration strategies for mitigating the impact of port disruptions
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
07/17/2009
Defense Date
05/19/2009
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
heuristics,intelligent transportation systems,mixed integer linear programming,OAI-PMH Harvest,optimization for container terminals and logistics network,set partitioning,traffic control
Place Name
USA
(countries)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Ioannou, Petros A. (
committee chair
), Moore, James Elliott, II (
committee member
), Safonov, Michael G. (
committee member
)
Creator Email
chang921@chol.com,hwanchan@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m2373
Unique identifier
UC1166126
Identifier
etd-chang-2980 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-561404 (legacy record id),usctheses-m2373 (legacy record id)
Legacy Identifier
etd-chang-2980.pdf
Dmrecord
561404
Document Type
Dissertation
Rights
Chang, Hwan
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
heuristics
intelligent transportation systems
mixed integer linear programming
optimization for container terminals and logistics network
set partitioning
traffic control