2011Bc.LukáˇsRapant ApplicationofDynamicGraphsfortheFastestRouteSearchinaTransportNetworkVyuˇzit´ıdynamickýchgraf˚uprovyhledáván´ınejrychlejˇs´ıcestyvdopravn´ıs´ıti VˇSB–TechnicalUniversityofOstravaFacultyofElectricalEngineeringandComputerScienceDepar

(1)

Department of Applied Mathematics

Application of Dynamic Graphs for the Fastest Route Search in a

Transport Network

Vyuˇzit´ı dynamick ých graf ˚u pro vyhled áv án´ı nejrychlej ˇs´ı cesty v

dopravn´ı s´ıti

Diploma Thesis

2011 Bc. Luk ´aˇs Rapant

(2)

Ostrava6^thMay 2011 . . . .

(3)

support.

(4)

V práci budou k pˇredloˇzeny tˇri algoritmy ˇreˇs´ıc´ı tento problém: jeden pracuj´ıc´ı na bázi Dijkstrova algoritmu, druh ý pracuj´ıc´ı na bázi A* algoritmu a tˇret´ı pracuj´ıc´ı na bázi mra- venˇc´ıch algoritm ˚u.

Kl´ıˇcov ´a slova: silniˇcn´ı navigace, Dijkstr ˚uv algoritmus, mravenˇc´ı algoritmy, A* algoritmus, dynamick ´y graf, nejrychlejˇs´ı cesta

Abstract

This work deals with finding the fastest route in a dynamic graph representing road network. Its goal is to develop algorithms capable of finding such route. There will be three algorithms presented in this work: the first based on Dijkstra algorithm, the second based on A* algorithm and the third based on Ant colony optimization algorithm.

Keywords: road navigation, Dijkstra algorithm, ant colony optimization, A* algorithm, dynamic graph, fastest route

(5)

f(x) – a real function of a real variablex

R f(x)dx – an integral of the functionf with respect tox

ACO – Ant colony optimization

BFS – Breadth-first search

(6)

1 Introduction 5

1.1 Motivation . . . 5

1.2 Literature overview . . . 6

1.3 Thesis outline . . . 7

2 Theoretical background 8 2.1 Graph theory definitons . . . 8

2.2 Dijkstra algorithm . . . 10

2.3 A* algorithm . . . 12

2.4 Ant colony optimization . . . 14

3 Extension of graph theory to dynamic problems 17 3.1 Update method . . . 19

3.2 Function method . . . 20

4 Implementation of algorithms 26 4.1 Dijkstra-based algorithm . . . 26

4.2 A*-based algorithm . . . 36

4.3 ACO-based algorithm . . . 39

5 Results and their discussion 47 5.1 Comparison with literature . . . 55

6 Conclusion 56 7 References 57 Appendix 59 A Complete source codes 59 A.1 Dijkstra-based algorithm . . . 59

A.2 A*-based algorithm . . . 61

A.3 ACO-based algorithm . . . 64

A.4 Weighted A* algorithm . . . 68

A.5 Functions v(s) source codes . . . 69

B Other appendices 70

(7)

List of Figures

1 Progress of BFS algoritm . . . 11

2 Behavior of ants 1 . . . 15

5 Car traffic distribution by time of day on all roads . . . 19

6 Progress of update method 1 . . . 20

7 Progress of update method 2 . . . 20

8 Example of functionv(s) . . . 21

9 Example of average speed profile in time . . . 22

10 Example of functionφ(t). . . 23

11 Example of functionϕ(t). . . 23

12 Graph used in the example 4.1 . . . 30

13 Example of progress of Dynamic Dijkstra algorithm 1 . . . 31

20 Examples of sigmoid and inverse sigmoid . . . 38

21 Delimitation of A 1 . . . 40

22 Delimitation of A 2 . . . 40

23 Graph representing traffic network of Moravian-Silesian region . . . 48

24 Comparison of results of Dijkstra-based algorithm and ACO-based algorithm . . . 48

25 Error of results of ACO-based algorithm . . . 49

26 Number of expanded vertices by Dijkstra-based algorithm and A*-based algorithm . . . 50

27 Difference in number of expanded vertices by Dijkstra-based algorithm and A*-based algorithm . . . 50

28 Time difference between dynamic and static path changes during a week (in hours) . . . 53

29 3D plot of time difference between dynamic and static path changes during a week (in hours) . . . 53

30 Path difference between dynamic and static path (green denotes static path and red dynamic path) . . . 54

31 Magnified part of the Figure 30 showing differences in paths in a different time (green is a static path and red are dynamic ones) . . . 54

32 The fastest path from vertex 100 to vertex 500 at 00:00 on Monday . . . 55

(8)

1 Car traffic density distribution by time of day on all roads (in percent of average car trafic density) . . . 18 2 Structure and example of struct used to represent list of neigbors . . . 27 3 Comparison of results of dynamic algorithm with different discretization 51 4 Time difference between dynamic and static path changes during a week

(in hours) . . . 52

(9)

List of Source Codes

1 Dijkstra algorithm pseudocode [17] . . . 11

2 A* algorithm pseudocode [18] . . . 13

3 ACO pseudocode . . . 16

4 Dynamic Dijkstra algorithm part 1 . . . 27

9 Dynamic A* code . . . 37

10 Dynamic ACO code part 1 . . . 41

21 Dijkstra-based algorithm source code . . . 59

22 A*-based algorithm source code . . . 61

23 ACO-based algorithm source code . . . 64

24 Weighted A* algorithm . . . 68

25 Functions v(s) source codes . . . 69

(10)

1.1 Motivation

In last 15 years number of cars in Czech Republic almost doubled and is still steadily growing. It depicts importance of a road traffic in our country and a great potential car navigation systems posses. Let us outline how car navigation works and what features it has. When you turn on the navigation it receives information from the GPS satellites of its current position. Then you choose your destination and the navigation computes by some algorithm shortest or fastest route that will take you to your desired destination. GPS also continually fix to your position so you know exactly where on the route you are and navigation tell you in advance where you are to turn or in which lane you should drive.

Nowadays navigations also give you a lot of other information about petrol stations, restaurants and hotels and other points of interest. Some navigation devices also receive RDS-TMC channels that inform you about traffic problems along your route. Most of new smart cell phones already contain GPS and car navigation software, so car navigation is becoming more and more common.

This thesis deals with algorithms that find the fastest route from your position to chosen destination. These algorithms are based on the Graph theory. Traffic network is transformed into graph and some graph algorithm for finding the fastest/shortest route is applied. Resulting route is relayed to the driver. It may seem simple at the first glance but several issues arise. First of them is that graph generated from traffic network would be too large. Generally navigations are small and not tooled for computational power.

It is necessary to separate such a big graph to smaller ones. This is done by partitioning (see e.g. [1]). Graph representing traffic network is divided into several smaller sub- graphs which contain roads of the same level (for example motorways, 1st class roads, 2nd class roads etc.) connected by some boundary graph. This allows searching on rela- tively smaller sub graphs. Search is first done on higher level roads which are faster and then moves to lower levels as faster roads usually do not lead to hotels, places of business etc. Also some heuristics can be added to solve this problem. Second issue, which is main motivation for this thesis, is what graph algorithm for finding fastest/shortest route to use.

Today’s most used algorithms are Dijkstra shortest path algorithm (see part 2.4) and its variants like A* (see part 2.3). While these algorithms perform very well they have some problems. They take traffic network as network with static edge labeling. While it may seem correct because distances and quality of roads usually do not change it is not so. Travel times are influenced by time of the day and to a lesser extent by day of a week.

Indeed it is very different to drive on a motorway at midnight or at 16:00. This thesis will work with dynamically evaluated graphs. These graphs better represent real traffic network than static ones. Aim of this thesis is to propose and analyze algorithms for finding the fastest route on these graphs. We will not deal with shortest routes because they do not need any time information as distance is never influenced by time. They can be found by using common Dijkstra-type algorithms. Proposed algorithms should

(11)

be capable of finding the fastest possible route taking time of a day and corresponding traffic density into account.

1.2 Literature overview

Shortest/fastest path problem is studied since 1950s. One of the first and one of the most well-known was developed by Dijkstra [2] and so is called Dijkstra’s algorithm. On this algorithm is based another widely used algorithm that is called A* [3]. It adds heuristic function which estimates cost of the route from current node to the destination. If the estimation satisfies some conditions it can be used to effectively find shortest route.

Many other algorithms were proposed like Bellman-Ford algorithm [4] [5] which allows search for shortest path even in case that this graph contains negative edge weights and Floyd algorithm [6] which also accepts negative edge weights but can only solve all-pairs shortest paths and provides only distance and not the path.

Important part of traffic navigation is also speeding up the search process because Dijkstra algorithm works slowly on large graphs. Jung and Pramanik [7] describe method for graph partitioning. They divide large graph into sub graphs which are connected by boundary graph. Another method is using planar separators invented by Lipton and Tarjan [8]. This was practiced for example by Chen and Xu [9].

There were a lot of research activity on behalf of the dynamic graphs and routing on them. Harary and Gupta [10] layed down some basics of dynamic graphs and outlined some applications for them. They also study them as extension of static graphs. Kali- nowski and Wenk [11] present model where navigation receives periodic updates and computes only shortest path to the next update. It works on base of A* algorithm. This algorithm tries to be at high degree node when update comes so it has a lot of possi- bilities to react on the change in graph. They also show that these computations can be done locally on navigation or on some central computer connected to navigation. Deme- trescu, Finocchi and Italiano [12] have come up with similar approach. Dynamism is also represented by some periodic updates to the edges. Their all pairs shortest path recom- putes only shortest path which were invalidated by update. After the update they are recomputed by Dijkstra algorithm and by observation that shortest route is composed of other shortest routes. These approaches are unsuitable for our problem because they are not good at modeling periodical processes. Useful insights are also provided by Tong and Wong [13]. They point to the problem that updates of edge value must be done even when a car is already travelling on the updated edge and the algorithm must consider this change. Not considering it can provide results that are very different from the reality on long roads. This is because traffic is changing continuously even when you are driving on a road without exits or ramps. Flinsberg [2004] in her dissertation used modified A*

algorithm in conjunction with daily traffic speed profiles representing time dependence of the graph. However she computed speed on an edge as constant.

Ant’s colony optimization (ACO) and its application on problem of finding the shortest route was also discussed in some papers like Angus and Hendtlass [14]. They show that ACO can compete with deterministic algorithms but has several weaknesses of sto-

(12)

produce some solution.

1.3 Thesis outline

In the next part of this thesis we define all theoretical terms that we use in this thesis.

These include definitions from the Graph theory, mathematical analysis, and descriptions of some basic graph algorithms like graph searches and basic Dijkstra algorithm. We will also describe here Ant colony algorithms which we will use later.

After theoretical terms follows part that is about making static graphs dynamic. We will show some representations of time in dynamic graph and discuss their advantages.

Then we will decide which one we will use for our algorithms.

After that we will describe proposed algorithms. We will present their features and talk about their implementation. They will be demonstrated on some simple examples solved by our implemented algorithms and then manually for comparison of results.

Next part comprises of analyze and comparison of proposed algorithms. We will test all presented algorithms on real data obtained from local traffic network and do thorough analysis. Results of this analysis will be discussed and cons and pros of each algorithm presented.

In the final part we will summarize our results and achievements. This part will be followed by appendix containing proposed algorithms implemented in Matlab programming language.

(13)

2 Theoretical background

2.1 Graph theory definitons

In this part we present graph theory definitons needed for defining algorithms, that we use in this thesis [25].

2.1.1 Undirected Graph

Definition 2.1 Undirected graph is an ordered pairG = (V, E)comprising a setV of vertices or nodes together with a setE of edges or lines, which are 2-element subsets ofV.

2.1.2 Weighted Graph

Definition 2.2 Weighted grapg is graph that has numeric weightw(e)associated with each edge einE. Edge weights can be integers, rational numbers or real numbers. This number represents some property like length or cost.

2.1.3 Walk

Definition 2.3 A walk in a graphGis such sequence of edges a vertices (v0, e₁, v₁, e₂, v₂, ..., e_n, v_n),

thate_i has adjacent verticesv_i−1 andv_ifor eachi= 1,2,3, ..., n. Such walk is called a(v₀, v_n)- walk.

2.1.4 Trail

Definition 2.4 Draw is special kind of walk, in which no edge is repeated. We call it a(v0, v_n)- trail.

2.1.5 Path

Definition 2.5 Path is a trail, in which no vertex is repeated. We call it a(v₀, vn)-path.

2.1.6 Length of a path

Definition 2.6 Length of a path is determined by the by number of edges that are traversed on a path. In weighted graph is this number is denoted as sum of weights of edges

l(u, v) =

n

X

i=1

w(e_i), wherew(ei)is weight of an edgeei, that belongs to a path.

(14)

Definition 2.7 Distance from vertex v₀ to vertex v_n in graph G is a length of the shortest (v₀, vn)-path. We denote it asdist(v₀, vn). If there is no(v₀, vn)-path then we setdist(v₀, vn) =

∞.

2.1.8 Adjacency matrix

Definition 2.8 Adjacency matrixA(G) is representation of a graph G. It is symmetric matrix n×nwherenis number of vertices in graph G. Each elementa_ijhas value determined by

a_ij =n

1if vivj∈E(G) 0otherwise,

whereE(G)is set of edges of graphG. In case of weighted graph this is modified to aij =n

w(vivj)if vivj∈E(G) 0otherwise.

wherew(v_iv_j)is the weight of an edge betweenv_iandv_j andE(G)is set of edges of graphG. 2.1.9 Adjacency list

Definition 2.9 Adjacency lista[]is representation of a graph G. Every element of lista_i[]stores vertices that are adjacent tov_i, for example

a₁[] = [v₂, v₃, v₄].

In case of a weighted graph we add a second row to the matrix with adjacent vertices containing weights of corresponding edges

a₁[] =

v₂, v₃, v₄ 2,5,6

.

(15)

2.2 Dijkstra algorithm

Dijkstra algorithm is used for solving single-source shortest paths problem. Single-source shortest paths problem is defined as finding shortest path from vertexvito all other vertices in a graph. Dijkstra algorithm was developed by Esger Dijkstra in 1956 and published in 1958 (see [2]).

Dijkstra algorithm is base on the idea of greedy algorithms. A greedy algorithm is some algorithm that follows the problem solving heuristic of making the locally optimal choice at each step with the hope of finding the global optimum. This approach can be generalized into five properties that greedy algorithm must posses to work [15] [16]:

• A candidate set, from which a solution is created

• A selection function, which chooses the best candidate to be added to the solution

• A feasibility function, that is used to determine if a candidate can be used to contri- bute to a solution

• An objective function, which assigns a value to a solution, or a partial solution, and

• A solution function, which will indicate when we have discovered a complete solution

Also problems on which greedy algorithms can be used must have two properties. The first is that it must have Greedy algorithm property. It means that during sovling the problem algorithm never needs to reconsider its previous steps. The second property is optimal substructure. It means that optimal solution of the problem consists of optimal solutions of the subproblems. Problem on which greedy algorithms cannot be used is for example Change-making problem. Let us suppose we have coins with values of 25c, 10c and 4c. We want to use as few coins as possible to have 41c. Greedy problem uses 25c coin, 10c coin and then fails because it cannot combine any coins to produce 6c. Optimal solution is 25c coin and four 4c coins. On the other hand the single-source shortest path problem is typical problem with these properties.

Dijkstra algorithm uses breadth-first search so it will be useful to present it here. This search algorithm is used for visiting all vertices in a graph and is used by other algorithms that need to test every vertex for some condition. The algorithm starts at starting vertex and visits neighbor vertices. Then it checks these neighbor vertices one at a time and visits their neighbors creating another layer of vertices. After all vertices neighboring to the starting one has been checked the algorithm moves onto next layer. It countinues like this until all vertices in the graph has been visited. For clearness and better understanding see Figure 1.

Now we will show how Dijkstra algorithm works. Let us have an undirected and weighted graph. The first of all we will mark all vertices as unvisited and set shortest paths to all vertices as∞except for the starting one which is set to 0. Then we will start searching in starting vertex from which we are to find shortest routes. We will mark it as visited and do breadth-first search. For each vertex found by breadth-first algorithm

(16)

Figure 1: Progress of BFS algoritm

we will check whether distance from current vertex plus weight of the edge connecting these two vertices is smaller than current shortest path to found vertex. We call this as triangular inequality. If it is so, new shortest path is saved. After checking all vertices found by breadth-first search, algorithm will mark vertex with shortest distance from starting one as visited and repeat the procedure. Algorithm proceeds as long as there are any unmarked vertices. For better understanding we present pseudocode for Dijkstra algorithm [17]:

functionDijkstra (Graph, source):

for each vertex v in Graph: // Initializations

dist [v] := infinity // Setting starting shortest routes to infinity dist [source] := 0 // Setting distance to starting point C := theset of all nodes in Graph // Set of all unvisited nodes

whileC is not empty: // The main loop u := vertex in C with smallest dist []

if dist [u] = infinity :

break //all remaining vertices are inaccessible from source remove u from C // Mark vertex u as visited

for each neighbor v of u:

alt := dist [u] + dist between(u, v)

if alt <dist [v ]: // test triangular inequality dist [v] := alt

returndist []

Source code 1: Dijkstra algorithm pseudocode [17]

It is interesting, that every time we mark some vertex as visited, we have already found the shortest path from starting vertex to this marked one. This is because Dijkstra algorithm applied on the graph without negative edge weights allows greedy method mentioned earlier to work properly. Proof of this can be found for example here [15].

Now we will analyze time complexity of Dijkstra algorithm. We denote with n and m, the number of vertices and edges of the input graph G, respectively. We assume that the edge weights can be added and compared in constant time. It is also important to mention data structures used for storing the graph. For Dijkstra in most cases it is desirable to use

(17)

adjacency list because it allows us to go through all neighbor vertices in time equivalent to their number. We also represent set of unvisited vertices as a vector. This means that in our case search for the vertex with shortest distance form starting vertex it is a linear search through entire set of unvisited vertices. Taking all this into consideration we can define time complexity of basic Dijkstra algorithm with data structures mentioned above asO

|V|²+|E|

=O

|V|²

where V is set of vertices of a graph and E is set of edges of a graph.

2.3 A* algorithm

A* is best-first search algorithm. Best-first search means that it explores a graph by ex- panding the most promising vertex chosen according to a specified rule. A* finds the shortest path from starting vertex to one other vertex. It was proposed by Nils Nilsson [3] as extension of Dijkstra algorithm.

A* is basicly Dijkstra algorithm which adds some heuristic function f(x) to speed up the search. This function composes of functiong(x)which gives the length of a path from starting vertex to current vertex andh(x) which is heuristic estimate to the goal.

Functionh(x)must be andmissible heuristic.

Definition 2.10 Let us denotevto be some vertex of graphG,h(x)to be a heuristic function, h(v)to be cost indicated byh(x)to reach a goal fromvandC(v)to be the actual cost to reach a goal fromv. Then heuristic functionh(x)is admissible if it satisfies

∀vεV (G) :h(v)≤C(v)

Such admissible heuristic can for example be Euclidean distance to the goal as it is geome- trically shortest possible distance between two vertices. This also guarantees optimality of the result. But if we want to use closed set of visited vertices (i.e. we will store visited vertices in a closed set and do not reconsider shortest paths to them further in the algorithm) like the one in Dijkstra, the heuristic function h(x) must also satisfy additional condition given in the following following for a warranty of optimality.

Definition 2.11 Heuristic functionh(x)is consistent if

∀u, vεV (G) :h(u)≤d(u, v) +h(v), wheredis cost of the shortest path from vertexuto vertexv, holds true.

Algorithm as such is very similar to Dijkstra algorithm which we presented earlier.

When we use heuristich(x) = 0it turns into Dijkstra algorithm. It differs mainly in that distance is replaced by heuristic function f(x) as we denoted it before. However this holds true only for choosing next vertex to be explored, otherwise triangular inequality is tested as in Dijkstra algorithm. For clearness we present here pseudocode for A* algorithm [18].

(18)

functionA∗(start,goal)

closedset := the emptyset // Thesetof nodes already evaluated.

openset :=setcontaining the initial node // Thesetof tentative nodes to be evaluated.

came from := the empty map // The map of navigated nodes.

g score[ start ] := 0 // Distance from start along optimalpath.

h score[ start ] := heuristic estimate of distance ( start , goal)

f score [ start ] := h score[ start ] // Estimated total distance from start to goal through y.

whileopenset is not empty

x := the node in openset having the lowest f score [] value if x = goal

returnreconstruct path(came from, came from[goal]) remove x from openset

add x to closedset

foreach y in neighbor nodes(x) if y in closedset

continue

tentative g score := g score[x] + dist between(x,y) if y not in openset

add y to openset

tentative is better := true else if tentative g score <g score[y]

tentative is better := true else

tentative is better := false if tentative is better = true

came from[y] := x

g score[y] := tentative g score

h score[y] := heuristic estimate of distance (y, goal) f score [y] := g score[y] + h score[y]

return failure

functionreconstruct path(came from, current node) if came from[current node] isset

p = reconstruct path(came from, came from[current node]) return(p + current node)

else

returncurrent node

Source code 2: A* algorithm pseudocode [18]

There are some special versions of A* algorithm. In this work we use Weighted A*.

It is basically A* with inflated heuristic function hw(x) = α·h(x), α > 1. This speeds up a search (as fewer nodes are expanded) but sacrifices optimality. Shortest route found by this algorithm is at worst α times longer than optimal one. This algorithm is best

(19)

for solving problems where we need some rough estimation of shortest route as fast as possible.

Time complexity of A* algorithm cannot be easily defined as it depends on heuristic function used. Algorithm time complexity is polynomial in case when the search space is a tree, there is a single goal, and heuristic functionh(x)satisfies

|h(x)−h^∗(x)|=O(log(h^∗(x))),

whereh^∗(x) is the exact cost to get from x to the goal. Otherwise it is in worst case exponential.

2.4 Ant colony optimization

This is a class of algorithms based on behavior of colonies of ants. It was proposed by Marco Dorigo in 1992 in his PhD thesis [18] and further developed by him and V.

Maniezzo and A. Colorni. It is not a single algorithm with clear definition but entire spectrum of algorithm using the same principles. Ant colony algorithms are regarded as populated metaheuristics with each solution represented by an ant moving in the search space. Ants mark the best solutions and take account of previous markings to optimize their search.

Let us present natural behavior of ants on which these algorithm is based. We have an ant colony C and a food source F. Ants travel between these places and on the way back from F to C leave small doses of pheromone along their trail for navigating other ants.

Pheromones are chemicals capable of acting outside the body of the secreting individual to impact the behavior of the receiving individual. Other ants then follow trails with strongest concentration to reach F. This can be generalized into following points.

1. Ant leaves C in search of food.

2. As it is first ant, it searches randomly neigborhood of C.

3. When it finds F, it will go back to the C leaving pheromones along its path.

4. Other ants then search for food not entirely randomly, but follow strongest pheromone trails. However some ants will still sometimes follow other paths exploring for new and possibly better ways to the F.

5. As time passes pheromones vaporizes which highlights most used trails, where pheromone level is still high.

We will show some example. Let us have a situation depicted in Figure 2. We have two alternative routes and the first ant (represented by dot). It randomly chooses route to pass through this situation and returns the same route placing pheromones. Next ants choose their way accordingly to the pheromone levels with some ants choosing other routes as in Figure 3. After some time the pheromone trail starts to evaporate, thus reducing its attractiveness. The more time it takes for an ant to travel down the path and back

(20)

C F

Figure 2: Behavior of ants 1

C F

again, the more time the pheromones have to evaporate. This leads to situation showed in Figure 4 as most of the ants choose shorter route. All this was proven by biologists in their experiments.

C F

Now we will use this model of behavior to define computer algorithm based on it.

Neighborhood of ant colony C is represented by a graphG(V, E)where some vertex v represents ant colony C and some vertexurepresents food source F. Ants through edges and vertices of the graph and on the edges leave pheromone trails. Ant is prohibited to enter any vertex more than once to prohibit it going in circles. Potential circles are not useful because shortest path will not contain any and they would only disattract other ants from seeking optimal solution. If the ant is prohibited to go any further, it is removed.

At each vertex ant must choose his next way represented by edge leading from that vertex. We determine probability of entering edge between verticesuandvfrom a [19]:

p_k(u, v) = δ(u, v) [τ(u, v)]^β P

y∈Nuδ(u, y) [τ(u, y)]^β, (1) whereδ(u, v)is pheromone level on the edge(u, v),τ(u, v)is reciprocal value of length of the edge(u, v),N_u as set of neighboring vertices of vertexu andβ > 0 is coefficient determining influence of edge weight.

(21)

Next we will show how pheromone laying and vaporization is done. Both processes are done by following formula

δ_k(u, v) = (1−α)δ_k(u, v) +△δ_k(u, v), (2) wherekis a number of ant,α∈(0,1)is coefficient of vaporization and

△δk(u, v) =

( ₁

Lif A_k

0otherwise, (3)

whereA_kis set of vertices traversed byk-th ant andLis cost or length of its path. This update is done after each ant has completed its journey or after batch of ants has completed their journey.

Now we will define how typical ant algorithm proceeds. It is not surprising that it almost mimics behavior of true ant. Proceeding of the algorithm is shown in the following pseudocode.

functionACO(start,goal,number of ants)

best journey=Inf; // setting starting distance to infinity

while(n=!number of ants) // algorithm termination by number of ants used setAnt( start ) ;

while(goal not reached)

applyProbabilityRule() ; // determining probabilities for neghbouring edges moveToChosenVertex();

end

if length(ants journey)<length(best journey) best journey=ants journey; // saving bestpath end

updatePheromones();

end

returnbest journey;

Source code 3: ACO pseudocode

It is good to emphasize that ACO are heuristic (or stochastic) algorithms so even several thousand ants are not guaranteed to find solution in case of huge problem. In such case the only thing we can do is use more ants, what can be time demanding. Also we can modify algorithm for better suiting some concrete problem.

ACO algorithms have many uses. The most common is finding shortest path in a graph, which is for what we use them. With some modifications they can be used for graph coloring, travelling salesman problem, sequential ordering, etc.

In terms of time complexity, ACO are classed as algorithms for solving NP-hard problems which means that their time is at worst exponential.

(22)

Now let us move from known results to our own work. Graph theory can be used to solve many real problems and routing is one of the more important ones. Routing problem can be defined as aa search for the shortest route from the starting position in some network (in our case road network) to desired goal position. The road network is represented by a weighted graph where edges represent roads, vertices their intersections and weights represent lengths of individual roads between the intersections. It is important to note that this graph is usually planar and every vertex has its constant position. It also means that the spatial configuration of the graph is important and we can use it in our calculations. To solve such problem we can apply one of the algorithms presented in the previous section. However this is not exactly accurate way to model routing problem. Graph theory takes entire graph as a static structure with static parameters. In case of traffic network this is not true. Many parameters, for example number of cars in the network which affects trafficability of roads, change in time can seriously influence the optimal solution. See Table 1 and Figure 5 for a notion how much can these parameters change in time. These statistics are car traffic density distribution by time of day on all roads in Great Britain in 2009 [20]. Numbers are percents of average traffic density. Here we should note one important thing. For all our calculations we have used data from various sources regarding various states (i.e. Great Britain, Netherlands, Czech republic) from various years. We know that this is not optimal, but data regarding traffic network are not publicly acessible so we have to use the best data we have found. However we assume that traffic networks in the developed states will have more or less same charac- teristics so our calculations are based on real situation.

Difference between 5% and 222% is self explaining. Also another issue is that not only parameters of network change but sometimes even the network itself can change as for example some road or intersection can be completely blocked by accident. Because of these reasons we need to introduce some dynamism to our model to represent these factors. So we will introduce timetinto our model. This will allow us to make parameters and structure of traffic network variable dependant on time t which is continous and has period of one week or mathematically<0,168>. This is because traffic is not only influnced by time of a day but also by day of a week. All this is summarized in the following definition.

Definition 3.1 Let us denoteG^∗(V, E, t) to be dynamic graph whereV is set of vertices,E is set of edges andtto be time . This graph has following properties:

• In each timetgraphG^∗ can have different sets V and E.

• If graphG^∗is weighted graph, it can have differentw(e, t)∀e∈Ein each timet.

This has influence not only on precision of our results but redefines our very task. In fact we are not usually looking for shortest route in traffic network but fastest (i.e. least time demanding route). But with introduction of timet, we must clearly state that we are looking for fastest route. In case when we are looking for the shortest route, we will only have to take into consideration graph structure changes.

(23)

Monday Tuesday Wednesday Thursday Friday Saturday Sunday

00:00-01:00 16 12 12 13 15 23 27

01:00-02:00 9 7 7 8 8 13 16

02:00-03:00 6 5 6 6 7 9 10

03:00-04:00 7 5 5 6 6 8 8

04:00-05:00 14 10 9 10 10 10 8

05:00-06:00 37 29 28 28 27 18 12

06:00-07:00 94 87 85 82 76 34 21

07:00-08:00 177 179 177 172 159 60 34

08:00-09:00 187 197 197 190 175 98 55

09:00-10:00 143 148 149 145 139 138 96

10:00-11:00 139 129 132 131 144 176 143

11:00-12:00 144 128 132 134 157 194 172

12:00-13:00 145 131 136 139 169 191 178

13:00-14:00 145 135 141 144 179 182 172

14:00-15:00 148 143 149 153 189 168 170

15:00-16:00 161 162 169 172 204 161 175

16:00-17:00 192 202 207 208 219 159 183

17:00-18:00 207 218 222 221 222 157 172

18:00-19:00 156 165 171 174 183 134 149

19:00-20:00 99 104 112 118 140 101 125

20:00-21:00 69 69 74 83 101 72 98

21:00-22:00 51 53 56 61 70 53 72

22:00-23:00 37 41 43 44 51 45 47

23:00-00:00 22 24 26 27 35 38 29

Table 1: Car traffic density distribution by time of day on all roads (in percent of average car trafic density)

(24)

Figure 5: Car traffic distribution by time of day on all roads

It is also desirable to mention that dynamism must not only be represented by changes of parameters of the graph in time but also in space. For example maximum speed on a road is not the same for the entire road. We can have some restrictions like maximum speed in town, in passages with a lot of bends, etc.

Now we will show two ways how to represent dynamism in a graph representing traffic network. The first one is developed by Demetrescu, Finocchi and Italiano [12] and the second one is developed by ourselves.

3.1 Update method

In its core, this method is quite simple. Its principle can be described by other definition of dynamic graph used by Demetrescu, Finocchi and Italiano [12].

Definition 3.2 Dynamic graph is a graph that undergoes a section of updates. An update is an operation that deletes or inserts vertices or edges, or changes some attributes associated with these features like weight (or color).

Unlike in our definition, this approach does not take time into consideration. Updates do not even have to be periodic. Update only happens when reasons for update are met.

For better understanding of this method, see following figures. In the left figure of Figure 6 you can see for example some traffic situation. Weights of the edges represent time that it takes to traverse them. In the right figure some accident happens on edge marked with !!! and the edge is blocked. In the left figure of Figure 7 you can see update reflecting

(25)

8 10

5 6

8 10

!!! 6

Figure 6: Progress of update method 1

10 12

4

8 10

5 6

Figure 7: Progress of update method 2

situation from the right figure of Figure 6. After accident is cleared we run another update which returns graph representing this traffic situation to normal state.

However for our traffic model this kind of approach is not suitable. It is great in representing dynamic events that happen only seldom and not periodically. As it is evident from Table 1 certain processes which affect the traffic are periodic and are happening constantly. This approach is good to solve problems like car accident but unsuitable in representing other more important processes affecting traffic like car density change. So for our model we proposed different approach which will be described in the following part.

3.2 Function method

This method instead of updates works with edge weight represented as functions. These functions represent time and space dynamism. This approach does not remove any edges or vertices but only works with edge weights. Basic idea is that we recomputed edge weights to represent time it takes to traverse them in certain time. It is important to know in what time we start to traverse the edge because in the dynamic graph parameters can change dynamically in time. Let us start with some physical interpretations. Time to get from place A to place B is by well known physic defined as

t= s v,

wheretis a time,vis a speed andsis a length of trajectory betweenAandB. This simple equation is however not suitable for us as it needsv to be homogenous. Therefore we have to redefine this relation by using integral calculus. This yields

time=

l

Z

0

1

v(s)ds, (4)

(26)

It gives us method to compute time it takes to traverse some road when the speed is not constant. This is very important relation for us and we will use it in our model.

This equation gives us a tool to represent space dynamism of the weight of some edge. We have a length of an edge (in our case representing length of a road). Next we will have to have function representing maximum speed of a car on the road. In our model this is represented by functionv(s)denoting maximum speed on each part of the road. We assume that the driver is not violating the traffic rules and he also drives as fast as possible so he travels by maximum speed. Such function can look like one in the Figure 8 which is denoted as

v(s) =







90if s∈ h0,2300) 50if s∈ h2300,3000) 70if s∈ h3000,3500i.

Figure 8: Example of functionv(s)

This function can for example represent some second class road running through a village and a wood. Now with equation (4) in a mind we can computetimethat represents time to traverse this road. We can compute it numerically by using Rectangle Rule (see [22]).

This calculation transforms weight of an edge from its length to time needed to traverse it. It takes into consideration inhomogeneous speed limit on the road represented by the edge. However this time is not influenced by other cars. It is time that it takes

(27)

single car to drive through the road. Now we must represent effect of the traffic as it changes in time. First of all we need some way to find these effects. For this we can use average speed profile in time. This profile could look like one in the figure 9 representing average speeds on Netherlands highways in 2003 [1]. This profile shows how the average

Figure 9: Example of average speed profile in time

speed changes in time of the day and a day of a week. We will transform this profile into functionv(t)by using a spline interpolation (see e.g. [21]) which is dependent on timet and has 168 hours period. From this function we can get information how much speed changes by dividing this profile function by maximum speed on its particular road. In this case it is120^km_h and resulting function looks like one in Figure 10. We will denote this function to beφ(t).

Now we want to apply this functionφ(t)on time we get from equation (4) to receive time it takes to traverse the edge with taking changing traffic density into consideration.

However we need to multiply time we have obtained from equation (4) to represent slower speed in higher traffic density. Values of functionφ(t) are between 0 to 1 so we need to make reciprocal value of this function. This gives us functionϕ(t)whose example is in Figure 11 and which is denoted as

ϕ(t) = 1 φ(t).

In relation to equation (4) this function gives us information of how maximum speed on a road represented by functionv(s) changes in time due to changing traffic density.

Now we must take into consideration one more factor. Speed is not affected by traffic

(28)

Figure 10: Example of functionφ(t)

Figure 11: Example of functionϕ(t)

(29)

by same measure on all roads. Speed certainly drops more on highway than on some village road in the peak hours so we must reflect this fact in our functionϕ(t). This fact could be easily modelled by speed profile for each type of the road. However we only have profile for Netherlands highway so we have developed following approximation method. While we know that it is not absolutely accurately describing reality we believe it decsribes it adequatly. We will useϕ(t)as reference function. We need to find adequate functionϕ_v(s)(t)for each functionv(s)which reflects maximum speeds on the road. This will be done by transforming referenceϕ(t)intoϕ_v(s)(t). We know that values ofϕ(t)are greater than one. So if we want to change them it is suitable to lower all values of function ϕ(t)by one as any operation with its values without this reduction could produce value lesser than one. This would mean that traffic density would increase maximum speed of the road. This is however unrealistic and undesirable. So after we lower values ofϕ(t)we must apply effect of different maximum speed. Then we will add one to return function values to be greater than one. This is done by following formula

ϕ_v(s)(t) = ϕ(t)−1

prms wams(v(s))

+ 1, (5)

whereprmsis maximum speed on reference profile road which is in our case120^km_h and wams(v(s))is weighted average of maximum speed on road represented byv(s) weighted by lengths of road segments with different maximum speed. This is because speed functionv(s)is not always constant and we want to take into consideration not only the highest maximum speed on a road represented byv(s)but all maximum speeds present on this road. Average is weighted to reflect how much of the road has certain maximum speed as following result will better reflect reality. Now we have functionϕ_v(s)(t) for road represented byv(s)whose values are lesser thanϕ(t)if weighted average speed on road represented byv(s)is smaller than maximum speed of a road of profile function or higher if otherwise. We will call this function traffic density function

Now we can finally denote how we calculate final time it takes to traverse the road.

Thistime_{f inal}is defined as

time_{f inal}=

t0+time

Z

t0

ϕ_v(s)(t)dt, (6)

where t₀ is time we enter the road, time is time we receive from equation (4) and ϕ_v(s)(t)is a traffic density function adapted for this particular road from equation (5).

In case of accident like in the Figures 6 and 7 we will change value of the speed functionv(s)in place where accident happened to very small number so resulting integrals will yield number close to∞and the edge will not be traversable.

We believe this approach to the dynamization of graph representing traffic network to be most accurately describing reality. Each edge have its functions v(s) and corre- spondingϕ_v(s)(t)and by these functions and application of equations (4) and (6) we can compute time it takes to traverse it. This allows our model to be as close to the reality as

(30)

the fastest route in such dynamic graph.

(31)

4 Implementation of algorithms

In this part we will discuss developed algorithms and their implementation. Algorithm based on Dijkstra algorithm will be presented first, followed by algorithm based on ACO.

For implementation of our algorithms we use Matlab programming language because it provides us all tools needed for computation and visualization of our experiments.

However for real application these algorithm should be programmed in some more common an faster programming language like C#.

First of all we shall discuss features that both algorithms have in common. It is a choice of graph representation. Choosing graph representation is very closely tied to algorithm which we want to use it and what we want to do with it in the algorithm. We mainly use two kinds of graph representation: Adjacency matrix and Adjacency list. Both of them have their cons and pros and we will present them here.

Adjacency matrix stores graph as a matrix (or in other words as two-dimensional array). This reveals one of it flaws. In algorithm where we do some graph search (breadth- first search especially) we must test entire row of the matrix regardless of number of actual neighbor vertices. In very large graph this can be very slow. This makes this kind of graph representation not very useful in conjunction with Dijkstra based algorithms. On the other hand it allows very easy editing of the graph. To delete some vertex you only remove corresponding row and column of the adjacency matrix. Editing edge weight is also very easy as you need to access only two values (as adjacency matrix is diagonally symmetric). Some algorithms are even based on adjacency matrix such as Floyd’s algorithm (see [6]).

Adjacency list stores graph as a list of neighbors (which can be represented by dynamic lists, some form of struct or in the simplest case as array of arrays). This means that it is powerful tool for finding neighbor vertices and therefore it is very useful for algorithm performing search for neighbors. You only need to access as many elements as there are neighbors. The only exception is when represented graph is almost complete where it is almost equally demanding as search in the adjacency matrix. Small downside of this representation is that changing the graph structure is not as easy as in matrix form.

Although there is no doubt that for our algorithms adjacency list is much better, we have for completeness and comparison implemented both ways of graph representation and we will show their comparison in the next part containing results and their discussion. However in describing the algorithms we will only show ones that use adjacency lists.

4.1 Dijkstra-based algorithm

This algorithm is based on Dijkstra algorithm and uses Function method for dynamization of the graph representing traffic network. As we mentioned earlier we will use adjacency list for graph representation. This adjacency list is implemented as Matlab struct array. Each struct is list of neighbor vertices and contains information about weight of edges leading to these vertices and their type. Format of this struct is depicted in Table 2 containing designation of each row of the struct and some example values.

(32)

edge to neighbor vertex weight 5.26 type or identifier of the road 2

Table 2: Structure and example of struct used to represent list of neigbors

Now let us talk about meaning of the last row. It contains information about type or ID of the road. This information is used to determine which functionv(s) will be used for computation of dynamism of the graph.

Now let us move to description how algorithm works. It proceeds like Dijkstra. It initializes variables in the same way. The only differences are starting time which is time in which we start our route and density function which represents function ϕ(t)from which are later in the algorithm computed invidualϕ_v(s)(t)functions. Algorithm returns vertices that the fastest path passes through and time it takes to traverse this path.

function[path, lenght of path ] = dijkstraplus (graph representation, starting vertex , goal vertex, starting time , density function )

% inicialization of varibles

number of vertices=length(graph representation);

visited verticies (1:number of vertices) = 0;

distances(1:number of vertices) = inf ;

previous vertex(1:number of vertices) = number of vertices+1;

distances( starting vertex ) = 0;

t0=starting time ;

terminate=ones(1,number of vertices)∗Inf;

time period=168;%period of function v(s)

space discretization=0.1; %discretization of function v(s) time discretization =60;%discretization of function phi( t ) max profile speed=120;%maximum profile function speed

Source code 4: Dynamic Dijkstra algorithm part 1

Now we move into the loop of the algorithm. Next it looks for the closest vertex to expand. This is however not choosing it by its distance but by the time it takes to get there.

There is also included a condition to terminate algorithm in case that no other vertex can be expanded.

%loop which terminates after goal vertex is marked as tested while visited verticies (goal vertex)˜=1

%choosing next expanded vertex candidates=[];

for i =1:number of vertices if visited verticies ( i )==0

candidates=[candidates distances(i)];

else

candidates=[candidates inf];

end end

if candidates==terminate break

end

(33)

[u index u]=min(candidates);

starting time =min(candidates);

%marking chosen vertex as visited visited verticies (u)=1;

%retrieving neigbor verticies form graph representation current vertex=graph representation{u};

[kk,number of neigbors]=size(current vertex);

%computing of current time in current vertex current time=rem(starting time+t0,time period);

Next part of the algorithm contains dynamization of the graph. It computes time it takes to get to each neighbor of expanded vertex in time we arrive to expanded vertex. This time is computed by equations (4) and (6). In our code we obtain the first function from the third row of the struct representing list of neighbors. We have simply divided edges to 4 classes equaling classes of the roads. Each function then represents maximum speed allowed on this road. We do not take into consideration that these functions are not constant because we do not posses so much detailed data regarding speed limits on each road. So we streamlined these functions. However this part of the code can be rewritten to choose by road ID contained in the same row of the struct instead of the road class. This would be useful for real application. After choosing right function forv(s)we compute time it takes to traverse the edge without taking traffic density into account. After that we compute equation (6) using functionv(x)chosen by road class of the computed edge.

All integrals are computed numercally using Rectangular Rule. After computing this integral we finally have time it takes to get from expanded vertex to some of its neighbor one.

%testing each neighbor for i =1:number of neigbors

%declaring vectors fo computing of space dependant integral current vertex length=current vertex(2, i ) ;

space integral discretization =current vertex length / space discretization ;

space integral vector=space discretization∗ones(1,floor( space integral discretization )+1);

space integral vector (floor( space integral discretization )+1)=rem(

space integral discretization, space discretization ) ; road coefficient =0;

% choosing right functions for cpace and time according to road type switch current vertex (3, i )

case 1

space function=’f1s’ ; case 2

space function=’f4s’ ; end

% computing space dependant integral for j =1:length(space integral vector)

road coefficient = road coefficient +feval(space function,j∗space discretization) ;

(34)

space discretization) ; end

space integral=sum(space integral vector);

%adjusting traffic density function

road coefficient = road coefficient /length(space integral vector) ; phi function =density function ;

phi function =phi function−1;

phi function =density function /( max profile speed/road coefficient ) ; phi function =phi function+1;

%declaring vectors fo computing of time dependant

%integral

time integral discretization =space integral /(1/ time discretization ) ;

time integral vector =(1/ time discretization )∗ones(1,floor( time integral discretization )+1);

time integral vector (floor( time integral discretization )+1)=rem(

time integral discretization ,(1/ time discretization ) ) ; current time seconds=floor(current time/(1/ time discretization ) ) ;

%computing time dependant integral

phi function vector =zeros(1,floor( time integral discretization )+1);

for w=1:floor( time integral discretization )+1;

phi function vector (w)=phi function(rem(current time seconds+w,time period /(1/time discretization))+1);

end

%calculation of final time it takes to traverse edge time final =sum(phi function vector.∗time integral vector ) ;

Time final that we get from following part of the code is applied into triangular inequality with time it takes to get to expanded vertex. If this is smaller than current time it takes to get to neighbor vertex it is saved as new fastest path.

if(distances(u)+ time final )<distances(current vertex(1,i ) ) distances(current vertex (1, i ) )=distances(u)+time final ; previous vertex(current vertex (1, i ) )=u;

end end end

Normally this loop terminates after all vertices are expanded. However having in mind properties of greedy algorithms we can terminate this loop when target vertex is expanded because it already has found the fastest path to it. After the fastest path is found we reconstruct it and save indexes of vertices it passes through.

path= [goal vertex ];

%recontruction of the path while path(1) ˜= starting vertex

if previous vertex(path(1))<=number of vertices path=[previous vertex(path(1))path];

end end

lenght of path = distances(goal vertex) ;

(35)

Source code 8: Dynamic Dijkstra algorithm part 5 Now we will show simple example how our algorithm works.

Example 4.1

Let us have a traffic network represented by graph in Figure 12 wheresdenotes starting position andgdenotes our goal destination. We start our journey on Monday at 8:00. Our goal is to find the fastest way fromstog.

s 2.5/1 2/2 4/1a

3/2a 2/3

1/1 1/1

1.5/2b

g

1.5/3 5

1 2

3

4

Figure 12: Graph used in the example 4.1

Each edge has its weight represented by two numbers: the first is length of the edge in km and the second is type of the road which influences functionsv(s) and therefore ϕ_v(s)(t). We have three types of the functionv(s)for corresponding road types 1, 2 and 3

v₁(s) = 90, v₂(s) = 70, v₃(s) = 50.

We also have 3 special functionsv(s)representing edges with unhomogenous function value from Figure 12. These edges are denoted by1a,2aand2band are defined as

v_1a(s) =

(90if s∈<2,4>

60if s∈<0,2), v_2a(s) =

(70if s∈<0.5,3>

20if s∈<0,0.5), v_2b(s) =

(70if s∈<1,1.5>

40if s∈<0,1).

Traffic density functionsϕ_v(s)(t)are computed accordingly to eachϕ_v(s)(t)function.

Now we have defined all functions needed and can start the search. We will set all distances to∞except one to the starting vertex and start testing starting vertex (Figure 13).

(36)

s 0 1

1 ∞ 0

2 ∞ 0

3 ∞ 0

4 ∞ 0

5 ∞ 0

g ∞ 0

s g 1 2

3 4

5

Figure 13: Example of progress of Dynamic Dijkstra algorithm 1

We shall compute corresponding (4)and (6) equations for each edge. Starting time t= 8as we start at 8:00.

t= 8 time₁=

Z2.5

0

1

v₁(s)ds= 0.02777

time₂= Z3

0

1

v_2a(s)ds= 0.06

time_{f inal}₁ =

t+time1

Z

t

ϕ_v₁_(s)(t)dt= 0.0456

time_{f inal}₂ =

t+time2

Z

t

ϕ_v_2a_(s)(t)dt= 0.1006

We can narrowly update the table as the times we get are smaller than∞and choose new vertex to expand (Figure 14). Timetis now done by adding starting time plus time it takes to get to the current vertex.

t= 8 + 0.0456 = 8.0456 time₄ =

Z2

0

1

v₂(s)ds= 0.0285

(37)

index distance visited

s 0 1

1 0.0456 1

2 0.1006 0

3 ∞ 0

4 ∞ 0

5 ∞ 0

g ∞ 0

s g 1 2

3 4

5

time₃ = Z2

0

1

v₃(s)ds= 0.04

time_{f inal}₄ =

t+time4

Z

t

ϕ_v₂_(s)(t)dt= 0.0472

time_{f inal}₃ =

t+time3

Z

t

ϕ_v₃_(s)(t)dt= 0.0532

Again we can simply update the table because we have not reached new vertices before and we choose the nearest vertex for another expansion (Figure 15). Timetis again computed as in previous step.

s 0 1

1 0.0456 1

2 0.1006 0

3 0.0988 0

4 0.0928 1

5 ∞ 0

g ∞ 0

t= 8 + 0.0928 = 8.0928

(38)

s g 1 2

3 4

time_g =

4

Z

0

1

v_1a(s)ds= 0.0555

time_{f inal g} =

t+timeg

Z

t

ϕ_v_1a_(s)(t)dt= 0.1082

We have path to goal vertex and as it is smaller than∞ so we save it into the table, choose another vertex to expand (Figure 16) and calculate timet.

s 0 1

1 0.0456 1

2 0.1006 0

3 0.0988 1

4 0.0928 1

5 ∞ 0

g 0.2010 0

s g 1 2

3 4

5

t= 8 + 0.0988 = 8.0988 time₂ =

1

Z

0

1

v₁(s)ds= 0.0111

(39)

time₅ =

1

Z

0

1

v₁(s)ds= 0.0111

time_{f inal}₂ =

t+time2

Z

t

ϕ_v₁_(s)(t)dt= 0.0213

time_{f inal}₅ =

t+time5

Z

t

ϕ_v₁_(s)(t)dt= 0.0213

Now we must test triangular inequality for vertex 2. However0.0213 + 0.0988>0.1006 so we need not to adjust time it takes to get to the vertex. For vertex five we have not yet found any other path to it so we only save the distance into the table, proceed to next vertex and compute new timet(Figure 17).

s 0 1

1 0.0456 1

2 0.1006 1

3 0.0988 1

4 0.0928 1

5 0.1201 0

g 0.2010 0

s g 1 2

3 4

5

t= 8 + 0.1006 = 8.1006 time₅ =

Z4

0

1

v₃(s)ds= 0.03

time_{f inal}₅ =

t+time5

Z

t

ϕ_v₃_(s)(t)dt= 0.0396