Simulating when a Green Light turns Red

January 9th, 2010

In a previous post, I talked about the Stochastic Traffic Cellular Automaton. I am finally going to post some simulations of this method. If you want to see the code, please leave a comment.

In this simulation, I will be modeling the traffic flow of an intersection when a green light turns red. From personal experience, you know that traffic will slowly get backed-up behind the red light. This can be described as a shock wave.

In Chapter 5 in Mark Holmes book “Introduction to the Foundations of Applied Mathematics”, he sets-up the exact same problem. If you want more details on the problem, you can download the chapter for free on his website.

In my simulation, I have discretized the road into 1000 segments and evenly distribute 250 cars as the initial data. I then set the max speed to 2 and the probability to slow-down to .25.

The following picture shows you how the simulation ran. Each row (left to right) represents the entire road where the black spots represent cars. As you move down, you will see how the road changed for each time-step. It might take you some time to make sense of the picture.

As you can see, the traffic on the right hand side starts to pile-up more and more as time goes on. This can be easily seen by the following movie.

The top plot shows the simulation for the Stochastic Traffic Cellular Automaton model. The plot below it shows the exact solution to the continuous Lighthill-Whitham PDE Model which I described in my poster presentation. As you can see, the two models behave very similarly.

Lastly, I ran the simulation 10,000 times. After running the simulations, I calculated the probability that a car exists in a particular road segment. This created a traffic density curve similar to the Lighthill-Whitham model which is shown in the following video.

Comments(0) - Discrete Math

Changing PDF Metadata on Ubuntu

January 4th, 2010

I got a Sony E-Reader as a Christmas present from my parents this year. I absolutely love it! It will definitely change how I read books. I no longer have to lug around several books everywhere I go. The best part is that I can put all my journal articles in it. I no longer have to print off hundreds of pages of journal articles. However, there is one thing that I have found really annoying. It is practically impossible to find programs that change the pdf’s metadata. It took me awhile to find one that actually worked.

I ended up finding a Linux Program called PDFTK. I wanted to go over some basic commands (mainly so that I can find them easily).

Installing PDFTK:

sudo apt-get install pdftk

Edit Existing Meta-Data:

pdftk book.pdf dump_data output report.txt

You can then edit the data in report.txt which we can later upload back to the pdf. The text file will contain key value pairs like:

InfoKey: Title
InfoValue: Coders At Work
InfoKey: Author
InfoValue: Peter Seivel
InfoKey: Subject
InfoValue: Programming

After you edit this file, you can update the new meta-data to the pdf.

Update PDF Meta-Data:

pdftk book.pdf update_info report.txt output bookcopy.pdf

Comments(2) - Software

Finite Difference Heat Equation using NumPy

January 2nd, 2010

To expand my experience with NumPy during the Christmas holiday, I decided to solve the heat equation using a finite difference technique. I know this is a trivial problem to solve. However, like I said in my previous post, the goal is to learn NumPy and not to solve complicated mathematical problems. Likewise, the purpose of this post is not to explain how finite difference methods work, but how to implement it using Python and NumPy.

The problem we are solving is the heat equation

$\frac{ \partial u}{\partial t} = \frac{ \partial^2 u}{\partial x^2}$

with Dirichlet Boundary Conditions ( u(0,t) = u(1,t) = 0 ) over the domain $0 \le x \le 1$ with the initial conditions

$u(x,0) = 10~sin( \pi x )$

You can think of the problem as solving for the temperature in a one-dimensional metal rod when the ends of the rod is kept at 0 degrees. Intuitively, you know that the temperature is going to go to zero as time goes to infinite.

To solve this problem using a finite difference method, we need to discretize in space first. I will be using a second-order centered difference to approximate $\frac{ \partial^2 u}{\partial x^2}$ . This will give the following semi-discrete problem:

$\frac{ \partial U}{\partial t} = A~\vec{U}$

where

$A = \frac{1}{\Delta x^2} \left( \begin{array}{cccccccc}-2 & 1 & 0 & 0 & \ldots & 0 & 0 & 0 \\1 & -2 & 1 & 0 & \ldots & 0 & 0 & 0 \\0 & 1 & -2 & 1 & \ldots & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & \ldots & 1 & -2 & 1 \\0 & 0 & 0 & 0 & \ldots & 0 & 1 & -2 \end{array} \right)$

and

$\vec{U}= \left( \begin{array}{c} u(0,t) \\ u(\Delta x,t) \\ u( 2 \Delta x, t) \\ \vdots \\ u(1-\Delta x,t) \\ u(1,t) \end{array} \right)$

The next step is to discretize in time. We can do this by using the Crank-Nicolson method which is

$\frac{\partial y}{\partial t} \approx \frac{ u_{n+1}-u_n }{\Delta t} = \frac{1}{2} ( f_{n+1}(x) + f_n(x) )$

where

$\frac{\partial y}{\partial t} = f(x)$

If we apply this method to the semi-discrete problem, we will get

$\frac{ \vec{U}_{n+1}-\vec{U}_n }{\Delta t} = \frac{1}{2} ( A~\vec{U}_{n+1} + A~\vec{U}_n)$

If we want to solve for $\vec{U}_{n+1}$ , we get the following system of equations

$(I -\frac{\Delta t}{2}~A)~\vec{U}_{n+1} = ( I + \frac{\Delta t}{2}~A)~\vec{U}_n$

We can implement this method using the following python code. This code will then generate the following movie.

Please note that the movie was sped up by a factor of ten.

Comments(0) - Numerical Calculations

Creating Movies with PyPlot/PyLab

January 2nd, 2010

Over the Christmas holidays, I started to dive more into Python. I have been trying to learn more about NumPy, SciPy, and PyLab. So far, I love these libraries. However, there doesn’t seem to be a good way to create movies/animations.

The best way to create a movie using Python is to generate plots using PyLab or PyPlot and save them as pictures. You can then use another program to stitch them together (ie. mencoder or ffmpeg). This process takes substantial amount of time and seems a little bit hacky at best. Please leave comments if you know of a better way.

It took some time to figure out how to get everything configured on MacOS X. My first problem was trying to get a copy of mencoder or ffmpeg. We need one of these programs in order to stitch all the pictures together. I found an interesting post which allowed me to easily install ffmpeg onto my mac. You can download my code and a copy of ffmpeg here.

Once you downloaded the zip file, you will need to

Copy ffmpeg to /usr/local/bin
Execute the following commands from your terminal:
sudo chown root:wheel /usr/local/bin/ffmpeg
sudo chmod 755 /usr/local/bin/ffmpeg

Once you have installed ffmpeg, we can look at how to create a movie. I created a function called CreateMovie which will easily create a movie

def CreateMovie(plotter, numberOfFrames, fps=10):
	import os, sys
	import matplotlib.pyplot as plt
 
	for i in range(numberOfFrames):
		plotter(i)
		fname = '_tmp%05d.png'%i
 
		plt.savefig(fname)
		plt.clf()
 
	os.system("rm movie.mp4")
	os.system("ffmpeg -r "+str(fps)+" -b 1800 -i _tmp%05d.png movie.mp4")
	os.system("rm _tmp*.png")

This function creates a movie called movie.mp4 in the current working directory which has a default frame rate of 10 frames per second. Make sure that you don’t have any files called movie.mp4 and _tmp*.png in the current working directory because they will be deleted.

The function takes the following parameters:

plotter: This parameter is a function with the following form:

def plotter(frame_number)

where frame_number is the current frame that needs to be plotted using the matplotlib.pyplot library.

numberOfFrames: The total number of frames in the movie.

fps: The frames per second. The default is 10.

In order to illustrate the use of this function, check out this code for a moving gaussian.

import CreateMovie as movie
import matplotlib.pyplot as plt
import numpy as np
 
x = np.linspace(-.5,4,500)
 
# Plots a given frame
def plotFunction( frame ):
	plt.plot(x, np.exp( -10*(x - frame/10.0)**2) )
	plt.axis((-.5,4,0,1.1))
 
movie.CreateMovie(plotFunction, 50)

This code will generate the following movie.

Comments(2) - Programming

Poster Presentation on Traffic Flow

December 17th, 2009

Lighthill-Whitham Model

The continuity equation for traffic flow is defined as:

$\frac{\partial \rho}{\partial t} + \frac{\partial q}{\partial x} = 0$

where $\rho$ is the car density, is the car speed, and is the traffic flow ( $q = \rho u$ ). This equation corresponds to a one dimensional fluid-dynamic model where no vehicles will enter or leave the traffic flow except at the boundary.

Let us assume that the car’s velocity depends only on the density ( $u=u(\rho)$ ). The simplest model
would be to consider the following linear relationship:

$u(\rho) = u_{max} ( 1-\frac{\rho}{\rho_{max}})$

where $\rho_{max}$ is the traffic density in which you have bumper-to-bumper traffic and $u_{max}$ is the maximum speed that a vehicle can travel. This usually corresponds to the speed limit.

You will notice that when there are no vehicles on the road an individual vehicle would travel at $u_{max}$ . Likewise, when you have bumper-to-bumper traffic, the vehicles on the road won’t be moving.

The relationship between density and traffic flow can be shown by the Fundamental Diagram of Road Traffic which has a maximum flow at $\frac{\rho_{max}}{2}$ .

If we put our velocity model into the continuity equation, we get the following:

$\frac{\partial \rho}{\partial t} + u_{max} ( 1-2 \frac{\rho}{\rho_{max}} ) \frac{\partial \rho}{\partial x} = 0$

Shock Waves

The Lighthill-Whitham Model is a non-linear wave equation. This equation can be solved using the method of characteristic lines. One of the defining features of this equation is the creation of expansion and shock waves. Shock waves are formed when characteristic lines intersect which causes $\rho (x)$ to become discontinuous. It can be shown that shock waves will have speed

$S(\rho_+,\rho_-) = \frac{q(\rho_+)-q(\rho_-)}{\rho_+-\rho_-}$

where $\rho_+$ and $\rho_-$ is the density when approaching from the right and left of the discontinuity.

Dealing with shock waves numerically has serious difficulties. To avoid shock waves, we can add a small diffusive term to the LW Model. If we define the traffic flow as

$q = u_{max} ( 1-\frac{\rho}{\rho_{max}} ) \rho-\epsilon \frac{\partial \rho}{\partial x}$

we can get a variation of burgers equation

$\frac{\partial \rho}{\partial t} + u_{max} ( 1-2 \frac{\rho}{\rho_{max}} ) \frac{\partial \rho}{\partial x} = \epsilon \frac{\partial^2 \rho}{\partial x^2}$

Fourier-Galerkin Spectral Method

The diffusive LW model can be solved using a Fourier-Galerkin Spectral Method. If we assume the problem has periodic boundary conditions and $x \in [0, 2\pi]$ , we can expand $\rho$ using $e^{inx}$ as a basis. This will lead us to the following system of equations

$\frac{d a_n(t)}{dt} + i~n~u_{max}~a_n(t) + \epsilon~n^2~a_n(t)$ $= \frac{2~i~u_{max}}{\rho_{max}} \displaystyle\sum_{j=-N/2}^{N/2} j~a_j(t)~a_{n-j}(t)$

$a_n(0) = \begin{cases}\frac{1}{2 \pi} \displaystyle\int_{0}^{2\pi} \rho(x,0)~e^{inx}~dx & \text{for $n \in [-\frac{N}{2}, \frac{N}{2}]$} \\ 0 & \text{otherwise} \end{cases}$

where

$\rho(x,t) = \displaystyle\sum_{j=-N/2}^{N/2} a_j(t)~e^{inx}$

After a Traffic Light Turns Green

Let’s analyze the behavior of traffic after a light turns green. We can do this by using the following initial conditions:

$\rho(x,0) = \begin{cases}\rho_{max} & \text{for $0 \leq x \leq \pi$} \\0 & \text{for $\pi < x \leq 2\pi$}\end{cases}$

If we use the Fourier-Galerkin Spectral Method, we get the following results

when $\rho_{max} = 3$ , $\epsilon = .01$ , $u_{max} = .5$ , and N=100 . You will notice the Gibbs phenomenon in our initial conditions. This is because a step-function can’t be exactly represented by our Fourier basis functions. However, the diffusive term will quickly damp any highly oscillating term.

If we look at the exact solution to this problem for the Non-Diffusive LW Model, we will get an expansion wave. The length of expansion can be shown to be $2 u_{max} t$ .

If we look at the Diffusive LW Model, we can see that the expansion length is dependant on $\epsilon$ .

After a Traffic Light Turns Red

Let’s analyze the behavior of traffic after a light turns red. We can do this by using the following initial conditions:

$\rho(x,0) = \begin{cases} \rho_0 & \text{for $0 \leq x \leq \pi$} \\ \rho_{max} & \text{for $\pi < x \leq 2\pi$}\end{cases}$

If we use the Fourier-Galerkin Spectral Method, we get the following results

when $\rho_{max} = 3$ , $\rho_0 = 1$ , $\epsilon = .01$ , $u_{max} = .5$ , and N=100 . After the Gibbs phenomenon dies down from our diffusive term, you will notice a shock that moves at speed -0.1685 .

If we look at the exact solution to this problem for the Non-Diffusive LW Model, we will get a shock that moves at speed

$u = - \frac{\rho_0 ~ u_{max}}{\rho_{max}-\rho_0} ( 1-\frac{\rho_0}{\rho_{max}} )$

If we use the same parameters as the above numerics, we get a shock speed of -0.166667 .

If we look at the Diffusive LW Model, we can see that the shock speed doesn’t vary much for different values of $\epsilon$ .

References:

[1] Kopriva, D. (2009). Implementing Spectral Methods for Partial Differential Equations. Springer Science
[2] Gottlieb, D. & Orszag, S.A. (1977). Numerical analysis of spectral methods: Theory and Applications. Montpelier,Vermont: Capital City Press
[3] Haberman, R. (1977). Mathematical Models: Mechanical Vibrations, Population Dynamics, and Traffic Flow. Englewood Cliffs, NJ: Prentice-Hall
[4] Helbing, D. (2001). Traffic and related self-driven many-particle systems. Rev. Mod. Phys. 73 (4), 1067–1141.

Comments(1) - Uncategorized

Interesting Books

Interesting Blogs

Interesting Articles

Simulating when a Green Light turns Red

Changing PDF Metadata on Ubuntu

Finite Difference Heat Equation using NumPy

Creating Movies with PyPlot/PyLab

Poster Presentation on Traffic Flow

Lighthill-Whitham Model

Shock Waves

Fourier-Galerkin Spectral Method

After a Traffic Light Turns Green

After a Traffic Light Turns Red