Excursions in Programming:
Traffic Engineering with Cellular Automata

Mathematica in Education and Research
Vol. 4 No. 3
Copyright 1995 TELOS/Springer-Verlag

by
Richard J. Gaylord and Kazume Nishidate

Introduction

Traffic engineering is concerned with the design of streets and highways and with the control of traffic on these roads. In an urban area, the routing of traffic needs to be considered when detours are set up so that road repairs can be made. When roadwork is done on a highway, the effect of temporary lane closures on traffic flow needs to be taken into account. That these matters are in fact, quite seriously regarded by the people in the Department of Transportation, is evidenced by the fact that roadwork is invariably carried out at times and locations that will mazimize both inconvenience and time delay (and flagmen are stationed to ensure that no one proceeds unimpeded through a work area). Cellular automata can be used to study both routing and traffic flow on streets and highways.

Solving a Maze

One aspect of routing is determining a course by which one can go from one location to another. This is closely related to the problem of finding a path through a labyrinth. We will show how this can be done using the following simple maze on a rectangular grid.

labyrinth = 
               {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, 
                {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1}, 
                {1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1}, 
                {e, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1}, 
                {1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1},  
                {1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}, 
                {1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1}, 
                {1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, e}, 
                {1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1}, 
                {1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1}, 
                {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},  
                {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};

In labyrinth,  1 represents a wall, 0 represents a path, and e indicates an entrance and an exit. The labyrinth is drawn below using red for walls, yellow for paths, and green for the entrance and exit.

Show[Graphics[RasterArray[Reverse[labyrinth] /. {e -> RGBColor[0, 1, 0],  
       1 -> RGBColor[1, 0, 0], 0 -> RGBColor[1, 1, 0]}] ], AspectRatio -> Automatic];

The usual procedure for solving a maze involves employing a 'search' strategy while the maze is being traversed (eg., start at the entrance and proceed along a path such that your left hand always remains in contact with a wall). While a search strategy is useful, indeed necessary, when a person is placed in a three-dimensional maze, there is a more efficient method for solving a maze map (a maze on a sheet of paper). Rather than tracing a path through the maze with a pencil (or a pen for the over-confident), it is possible to eliminate all of the cul-de-sac's (paths that culminate in a dead-end) in the maze so that all that remains is the path through the maze. This strategy can be easily implemented using a cellular automaton (CA).

Maze-solving cellular automaton

The maze is a rectangular lattice consisting of sites with value 1 (a wall site), 0 (a path site), or e (an entance or exit site). All  turns in the paths and the walls are at ninety degrees.

The CA rule for eliminating the site at the end of a dead-end path in the maze is:

A path site (a site with value 0) that has 3 or 4 nearest neighbor wall sites (sites with value 1) becomes a wall.

The CA rule for all of the other sites in the maze is: 

Any other path site (including the exit and entrance) and any wall site remains unchanged.

These rules are implemented in the following six rewrite rules:

mazeSolve[0, 1, 1, 1, 0]   := 1; 
mazeSolve[0, 1, 1, 0, 1]   := 1;
mazeSolve[0, 1, 0, 1, 1]   := 1;
mazeSolve[0, 0, 1, 1, 1]   := 1;
mazeSolve[0, 1, 1, 1, 1]   := 1;
mazeSolve[x_, _, _, _, _] := x;

The five arguments of the mazeSolve rules are the values of: a site, the north nearest neighbor to the site, the east nearest neighbor to the site, the south nearest neighbor to the site, and the west nearest neighbor to the site. 

The VonNeumannValues function defined below takes a matrix, lat, and a function, func, and returns a matrix whose elements are func applied to the value of the corresponding site in lat and its N, E, S and W nearest neighbor sites (ie., func[site, N, E, S, W]) [see Gaylord and Nishidate, 1995]

VonNeumannValues[func_, lat_] := 
              MapThread[func, {#, 
				RotateRight[#, {0, 1}], RotateRight[#, {1, 0}],
				RotateRight[#, {0, -1}], RotateRight[#, {-1, 0}]},
			2]&[lat]

The VonNeumannValues function is applied repeatedly to the mazeSolve rule and the maze lattice using the FixedPoint function until all of the dead-end paths have been removed.

FixedPoint[VonNeumannValues[mazeSolve, #]&, maze]

The program for the maze-solving cellular automaton is written as:

PathToEnlightenment[ maze_] := 
   Module[ {mazeSolve, VonNeumannValues},

 mazeSolve[0, 1, 1, 1, 0]   := 1; 
 mazeSolve[0, 1, 1, 0, 1]   := 1;
 mazeSolve[0, 1, 0, 1, 1]   := 1;
 mazeSolve[0, 0, 1, 1, 1]   := 1;
 mazeSolve[0, 1, 1, 1, 1]   := 1;
 mazeSolve[x_, _, _, _, _] := x;

VonNeumannValues[func_, lat_] := 
              MapThread[func, {#, 
				RotateRight[#, {0, 1}], RotateRight[#, {1, 0}],
				RotateRight[#, {0, -1}], RotateRight[#, {-1, 0}]},
                                  2]&[lat];
                                  
FixedPoint[VonNeumannValues[mazeSolve, #]&, maze]
   ]

We can use the PathToEnlightenment program and a maze input to create a graphic showing side-by-side, the maze and the path through the maze. This is done below for the labyrinth maze.

Show[GraphicsArray[
   Map[ Show[Graphics[
                    RasterArray[Reverse[#] /. {e -> RGBColor[0, 1, 0], 
                        1 ->RGBColor[1, 0, 0], 0 -> RGBColor[1, 1, 0]}] ], 
                     AspectRatio -> Automatic, DisplayFunction -> Identity]&, 
           {#, PathToEnlightenment[#]}&[labyrinth]] ]];

Traffic Flow

We'll look at the traffic flow on one side of a divided highway. The system consists of two lanes (with no entrances or exits along the lanes) containing cars and empty spaces (which are the same size as a car). This system is represented by a '2 x n' lattice whose sites have the value 1 (indicating a space occupied by a car) or 0 (indicating an empty space). Initially, each site is occupied, with probability p, by a car (for a very long road, the density of cars on the road will approach p).

road = Table[Floor[p + Random[]],{2},{n}]

Traffic flows in the same direction along both lanes of the road (ie., along both rows of the lattice). There is no distribution of car velocities in this simple model; at each step, a car either moves forward one space, moves sideways, or remains in place. When a car reaches one end of the road and then moves forward, it disappears from that end and re-appears at the other end (this behavior, which is implemented in the CA by using periodic boundary conditions, corresponds to having an infinitely long road). 

The behavior of a car is determined by presense of other cars in its vicinity.

=======================
                                              | 6 | 4 | 5 |       ->               
                              ----------------------------------------
                                              | 3 | 1 | 2 |       ->                
                              =======================

The rules for the behavior of a site (which is either empty or occupied by a car) during one time step in the highway CA depend on the values of the following sites (the enclosed numbers below correspond to the numbering of sites in the figure above): the site [1], the site ahead [2], the site behind [3], the site adjacent [4], the site ahead of the adjacent site [5], and the site behind the adjacent site [6]. The values of these sites appear (in the same order as they are numbered above) as arguments to the CA update rules which can be easily stated (in the rule descriptions that follow, the sites in the above figure that correspond to the sites in the rules are indicated in square brackets).

rule # 1 
an occupied site [1] becomes empty when the site ahead of the site [2] is empty (ie., a car that has an empty space ahead of it moves into the space).

drive[1, 0, _, _, _, _] :=  0

rule # 2 
an occupied site [1] becomes empty when the site ahead of the site [2] is occupied and the sites adjacent to the site [4] and behind the adjacent site [6] are both empty (ie., a car that is blocked by a car ahead of it, switches lanes when there is an adjacent empty space and an empty space behind that empty space).

drive[1, 1, _, 0, _, 0] :=  0

rule # 3 
an occupied site [1] remains occupied when the sites ahead of the site [2] and adjacent to the site [4] are both occupied (ie., a car that is blocked by cars ahead and next to it doesn't move).

drive[1, 1, _, 1, _, _] :=  1

rule # 4 
an empty site [1] that has an occupied site behind it [3] becomes occupied (an empty space ahead of a car becomes occupied by the car).

drive[0, _, 1, _, _, _] :=  1

rule # 5 
an empty site [1] that has an empty space behind it [3], an occupied site adjacent to it [4] and an occupied site ahead of the adjacent site [5] becomes occupied (ie., an empty space that is followed by an empty space and is next to a car that is blocked by a car ahead of it becomes occupied as the car switches lanes)

drive[0, _, 0, 1, 1, _] :=  1

rule # 6 
in situations other than those given above, empty sites remain empty and occupied sites remain occupied.

drive[x_, _, _, _, _, _] := x

The following anonymous function is used to apply the drive rules to the road lattice.

MapThread[drive, {#, RotateRight[#, {0, -1}], 
                                RotateRight[#, {0, 1}], RotateRight[#, {1, 0}], 
                                RotateRight[#, {1, -1}], RotateRight[#, {1, 1}]},
                                  2]&

To illustrate how this anonymous function works, we can use the following small roadTest lattice (the numbering of sites in the lattice corresponds to the numbering used in the figure above).

(roadTest= {{a, 6, 4, 5, b}, {c, 3, 1, 2, d}})//TableForm

Looking at what happens to the site numbered as 1 in roadTest when the anonymous function above is applied to roadTest and an undefined function, rule, we see that the sites that are numbered as 2 through 6 are taken, in order, as arguments to rule.

MapThread[rule, {#, RotateRight[#, {0, -1}], 
                                RotateRight[#, {0, 1}], RotateRight[#, {1, 0}], 
                                RotateRight[#, {1, -1}], RotateRight[#, {1, 1}]},
                                  2]&[roadTest][[2, 3]]

This anonymous function is used with the NestList function to repeatedly apply the drive rules to the road lattice t times.

NestList[MapThread[drive, {#, RotateRight[#, {0, -1}], RotateRight[#, {0, 1}],
                       RotateRight[#, {1, 0}],  RotateRight[#, {1, -1}], RotateRight[#, {1, 1}]},
                                  2]&, road, t]

We can write the program for the evolution of this CA as follows:

keepOnTrucking[n_, p_, t_] := 
   Module[{road, drive},

             road = Table[Floor[p + Random[]],{2},{n}]  ;

             drive[1, 0, _, _, _, _] =  0;
             drive[1, 1, _, 0, _, 0] =  0;
             drive[1, 1, _, 1, _, _] =  1;
             drive[0, _, 1, _, _, _] =  1;
             drive[0, _, 0, 1, 1, _] =  1;
             drive[x_, _, _, _, _, _] := x;
    
             
             NestList[MapThread[drive, {#, RotateRight[#, {0, -1}], 
                                RotateRight[#, {0, 1}], RotateRight[#, {1, 0}], 
                                RotateRight[#, {1, -1}], RotateRight[#, {1, 1}]},
                                  2]&, road, t]
                  ]

The dynamic behavior of traffic on one side of the two-lane divided highway can be seen in an animation. This is shown below for 10 successive time steps.

Map[Show[Graphics[RasterArray[# /. 
     {0 -> RGBColor[0.7, 0.7, 0.7], 1 -> RGBColor[0, 1, 0]}]],
  AspectRatio -> Automatic]&, keepOnTrucking[50, 0.7, 10]];

Road Obstacles

The effect of roadwork or traffic accidents on the flow of traffic in the two-lane divided highway model can be computed by making some modifications in the indy500 program. 

Note the ease with which these changes can be incorporated into the model; this is one of the main advantages of a rule-based programming approach. 

Using c to represent an obstacle, we can place a single obstacle at a random location along the road with

roadWithObstacle = ReplacePart[Table[Floor[p + Random[]],{2}, {n}],  c, 
                                                    {Random[Integer,{1, 2}], Random[Integer,{1, n}]}]

The following six rules are used to describe the movement of cars on the road when obstacles are present:

A car that has an empty space ahead of it moves ahead.

driveWithObstacle[1, 0, _, _, _, _] =  0

A car that is blocked by either a car or an obstacle ahead of it, switches lanes when there is an adjacent empty space and an empty space or an obstacle behind that empty space (note that, as in the previous model, a car does not switch lanes if there is a car in the empty space behind the empty adjacent space).

driveWithObstacle[1, 1 | c, _, 0, _, 0 | c] =  0

A car that is blocked by cars or obstacles ahead and next to it doesn't move.

driveWithObstacle[1, 1 | c, _, 1 | c, _, _] =  1

An empty space ahead of a car becomes occupied by the car.

driveWithObstacle[0, _, 1, _, _, _] =  1

An empty space that is followed by an empty space or an obstacle and is next to a car that is blocked by a car or an obstacle ahead of it becomes occupied as the car switches lanes.

driveWithObstacle[0, _, 0 | c, 1, 1 | c, _] =  1

Obstacles remain in place and in situations other than those given above cars don't move.

driveWithObstacle[x_, _, _, _, _, _] := x

The overall program for traffic flow in the presense of obstacles on the road is given by:

indy500[n_, p_, t_] := 
   Module[{roadWithObstacle, driveWithObstacle},

            roadWithObstacle = ReplacePart[Table[Floor[p + Random[]],{2}, {n}],  c, 
                                                    {Random[Integer,{1, 2}], Random[Integer,{1, n}]}] ;

            driveWithObstacle[1, 0, _, _, _, _] =  0;
            driveWithObstacle[1, 1 | c, _, 0, _, 0 | c] =  0;
            driveWithObstacle[1, 1 | c, _, 1 | c, _, _] =  1;
            driveWithObstacle[0, _, 1, _, _, _] =  1;
            driveWithObstacle[0, _, 0 | c, 1, 1 | c, _] =  1;
            driveWithObstacle[x_, _, _, _, _, _] := x;
    
             
             NestList[MapThread[driveWithObstacle, {#, RotateRight[#, {0, -1}], 
                                              RotateRight[#, {0, 1}], RotateRight[#, {1, 0}], 
                                              RotateRight[#, {1, -1}], RotateRight[#, {1, 1}]},
                                              2]&, roadWithObstacle, t]
                  ]

Computer Simulation Projects

1. Modify the keepOnTrucking program to calculate the mean velocity of the system - the number of cars moving ahead in a time step divided by the total number of cars [see Gaylord and Wellin 1995].

2. Use the program from the previous project to create the fundamental diagram of traffic flow - a plot of the flow or mean velocity of cars per time step versus the car density.

General References

Nayfeh, Basem A. Cellular automata for solving mazes. Dr. Dobb's Journal, February, p. 32 (1993).

Gaylord, Richard J. and Nishidate, Kazume. A simple rule-based programming approach to modeling complex random walk systems with cellular automata. Mathematica In Education and Research, 4:2  pp. 37-44 (1995).

Gaylord, Richard J. and Wellin, Paul R. Computer Simulations with Mathematica: Explorations in Complex Physical and Biological Systems. pp. 140-3, TELOS/Springer-Verlag (1995).