  # Simulating the Nagel-Schreckenberg model

In my Modeling, Simulation, and Decision Making (CS166) class, I learned about the Nagel-Schreckenberg model and how was used as a theoretical model to simulate the behavior of freeway traffic.

## Building the TrafficSimulation Class

First, let's create a `Traffic Simulation` class to model traffic behavior according to the Nagel-Schrekenberg model. This class will contain the following methods: - The initialization method `__init__()` to create the object - The `changeSpeed()` method, which is to simulate the behavior of the cars within the system - if there's no cars in front of them, they'll speed up; if there's a car in front of them, they'll slow down. - The `display()` method, which is just to show the model.

We will be running this model for 30 timesteps to observe the behavior of the Nagel–Schreckenberg model.

```import numpy as np
import random
```
```class TrafficSimulation():

def __init__(self, length = 100, density = 0.2, max_velocity = 5, prob = 0.5):
# Here, we define the model parameters (road length, traffic density, maximum velocity,
# probability of slowing down)
self.length = length # The length of the road (number of cells)
self.density = density # The density of traffic
self.max_velocity = max_velocity # The maximum speed of the cars
self.slow_prob = prob # The probability of slowing down for no apparent reason
self.state = -np.ones(self.length, dtype=int) # "-1" represents an empty cell
filled_cells = np.random.choice(range(self.length), size=int(round(density * self.length)), replace=False)
self.state[filled_cells] = np.random.randint(0, self.max_velocity + 1, size=len(filled_cells))

​
def changeSpeed(self, display = True):
for i in range(self.length):
if self.state[i] != -1: # We only change the speed if there's a car in the cell
dist = 0 # Algorithm that scans the grids in front of the car until it finds the next car
while self.state[(i + (dist + 1)) % self.length] == -1:
dist += 1
# Here we will code for the 3 rules for speed change in the Nagel-Schreckenberg model
if dist > self.state[i] and self.state[i] < self.max_velocity:
self.state[i] += 1
if dist < self.state[i]:
self.state[i] = dist
if self.state[i] > 0 and np.random.random() < self.slow_prob:
self.state[i] -= 1

if display == True: # Display the model
self.display()

# Here we are coding for the last rule in the Nagel-Schrekenberg model - the movement rule
next_state = -np.ones(self.length, dtype=int)
for i in range(self.length):
if self.state[i] != -1:
next_state[(i + self.state[i]) % self.length] = self.state[i]
self.state = next_state

def display(self):
print(''.join('.' if x == -1 else str(x) for x in self.state))

sim = TrafficSimulation()
for i in range(30):
sim.changeSpeed()
```

Output: ## Measuring the Flow Rate

Traffic flow is defined as the average number of cars passing per unit time between two particular cells. Although we can simply measure the flow between any two cells, it was advised that we consider the periodic boundary (the `self.length` parameter in the class above) as the two points. To do that, we could simply count the number of cars that exit the road on the right end of the road and re-enter on the left end in each time step.

The flow rate, on the other hand, is the average of all these flow counts over time. To measure flow rate, we will be adding extra parameters to the `__init__()` method as well as incorporating time steps into our model.

```class TrafficSimulation():

def __init__(self, length = 100, density = 0.2, max_velocity = 5, prob = 0.5):
# Here, we define the model parameters (road length, traffic density, maximum velocity,
# probability of slowing down)
self.length = length # The length of the road (number of cells)
self.density = density # The density of traffic
self.max_velocity = max_velocity # The maximum speed of the cars
self.slow_prob = prob # The probability of slowing down for no apparent reason
self.state = -np.ones(self.length, dtype=int) # "-1" represents an empty cell
filled_cells = np.random.choice(range(self.length), size=int(round(density * self.length)), replace=False)
self.state[filled_cells] = np.random.randint(0, self.max_velocity + 1, size=len(filled_cells))
self.time_step = 0 # Adding new parameters to measure traffic flow
self.total_traffic_flow = 0

def changeSpeed(self, display = False): # Changing this to False because we will be running this multiple times
# and we don't want the model to be printed every time
for i in range(self.length):
if self.state[i] != -1: # We only change the speed if there's a car in the cell
dist = 0 # Algorithm that scans the grids in front of the car until it finds the next car
while self.state[(i + (dist + 1)) % self.length] == -1:
dist += 1
# Here we will code for the 3 rules for speed change in the Nagel-Schreckenberg model
if dist > self.state[i] and self.state[i] < self.max_velocity:
self.state[i] += 1
if dist < self.state[i]:
self.state[i] = dist
if self.state[i] > 0 and np.random.random() < self.slow_prob:
self.state[i] -= 1

if display == True: # Display the model
self.display()

# Here we are coding for the last rule in the Nagel-Schrekenberg model - the movement rule
next_state = -np.ones(self.length, dtype=int)
for i in range(self.length):
if self.state[i] != -1:
next_state[(i + self.state[i]) % self.length] = self.state[i]
self.state = next_state

# For every time step, we add the traffic flow measured in that time step to total_traffic_flow
self.time_step += 1
for i in range(self.max_velocity):
if self.state[i] > i:
self.total_traffic_flow += 1

def display(self):
print(''.join('.' if x == -1 else str(x) for x in self.state))

sim = TrafficSimulation()
for i in range(50):
sim.changeSpeed()
```

## Plotting a graph of the traffic flow rates for different car densities

To do this, we will be running the simulation above 200 time steps, get the flow rate, repeat that for 20 trials, and then calculate the average flow rate. Then, we will repeat the whole procedure with varying densities from 0.01 to 1.00.

Sounds complicated? No worries, this will take Python just a couple of minutes to run...

```density_list = []
traffic_flow_list = []
for density in range(1, 101):
sim = TrafficSimulation(density=density/100)
# print('Traffic density:', sim.density)
density_list.append(sim.density)

trial_flow = []
for trial in range(20): # Repeating each time step by 20 times to get the average value
for i in range(200): # Run the simulation for 200 timesteps to get the average flow rate
sim.changeSpeed()
trial_flow.append(sim.total_traffic_flow / sim.time_step)
traffic_flow_list.append(sum(trial_flow)/len(trial_flow))
# print('Average traffic flow:', sum(trial_flow)/len(trial_flow))
```
```import matplotlib.pyplot as plt

plt.figure(num=None, figsize=(12, 9), dpi=80, facecolor='w', edgecolor='k')
plt.plot(density_list, traffic_flow_list)
plt.xlabel("Density (cars per site)")
plt.ylabel("Flow rate (cars per time step)")
plt.title("Graph of Traffic Flow Rate for Different Car Densities based on the Nagel-Schrekenberg Model")
plt.show()
``` ...and, voilà! As you can see in the graph above, the flow rate increases sharply all the way till the 0.12 (est.) mark, where it hits the optimum capacity of the system, before starting to decrease sharply and then gradually all the way down to 0 - which makes sense, because the cars on a fully packed road (`self.density = 1.0`) will not be able to move, and hence the flow rate will be 0.