Designing and prototyping your own simulation with Python3.

Written by Derek Groen (Derek.Groen@brunel.ac.uk), with help from Subam Bista, Diana Suleimenova and Gebremariam Assres.

In this tutorial, you will learn how to write a basic agent-based simulation application. The example we use is a simulation that attempts to predict the movement of persons escaping dangerous areas and seeking safety. The principles we cover in this tutorial have been used for a wide range of simulations, for instance to model the escape of people from armed conflicts (https://www.nature.com/articles/s41598-017-13828-9).

The underlying technique we introduce here is more widely known as agent-based modelling, or ABM.

In this tutorial we will cover the following aspects:

  • What is agent-based modelling in general.

  • Creating a simple agent-based model with Python3 and using several of the key concepts that we used in the Flee agent-based modelling code.

  • Observing uncertainty in your simulation results.

  • Basic steps towards designing and prototyping your own model.

Requirements

To do this tutorial, you need a working Python3 installation. Numpy and pyplot are useful, but optional.

Introduction: What is agent-based modelling?

agent-based model (ABM)

a computational technique to model the actions and interactions of autonomous agents, with a view to assessing their effects on the system as a whole. Agents may represent individuals, groups, or abstract entities.

In this tutorial we assume that we are modelling of people moving from one place to another. However, ABM can also be used to model other movements, e.g. of objects such as cars, of animals, of cells, or even transactions or e-mails. The elements may differ on the type of model you wish to create.

In the case of our people movement simulations, we work with three basic elements:

  • The persons themselves (Person class).

  • The locations where the persons reside (Location class).

  • And possibly the paths (or routes) that interconnect the locations (Link class).

In its simplest form, this agent-based model features people that reside at a given location, and that move from one location to another as the time in the simulation progresses.

Network-based versus geographically pixelated

In general there are two widespread basic approaches to ABM. One is network-based, where each location is an agent, and the location agents are interlinked using path agents. A second approach is geographically pixelated, where a region is subdivided into square areas, and the location of agents is indicated by the respective coordinates of the corresponding square areas. It is possible to use a network-based model to approximate a geographically pixelated model, by placing all locations on Carthesian coordinates.

The simulation code

What follows is a step-by-step explanation how you can use Python3 to build an agent-based simulation code. The code works as is, but as part of this tutorial you will be given the opportunity to change some of its features, so that you can make it suit your personal needs.

To start off, please make a new Python file named my-abm.py. Unless specified otherwise, please place any Python fragments following in this tutorial at the _bottom_ of this file.

To start off, we will need to import one library, so that we can use the randomizer.

import random

In this tutorial we use very few dependencies for the main simulation code, but the random library is an essential one, as agent-based simulations strongly rely on randomizers. We do require a few more dependencies for visualizing results, namely matplotlib, pandas, and ffmpeg for storing animations.

Defining the Location Graph

Now any Persons in the simulation will reside at a given place, or Location. To define these places in a networked model, we create a Location object for each place:

class Location:
  def __init__(self, name, x=0.0, y=0.0, movechance=0.001):
    self.name = name
    self.x = x
    self.y = y
    self.movechance = movechance
    self.links = []
    self.numAgents = 0

The Location class, has a number of simple parameters. These represent essential characteristics for individual locations:

  • name: the name of the Location.

  • x: GPS x-coordinate, useful for placing on a map and for calculating distances as the bird flies.

  • y: GPS y-coordinate.

  • movechance: An indicator denoting the safety level of this location. Are people certain to stay put (0.0), certain to move out immediately (1.0) or will there be a mixture (0.0<`movechance`<1.0).

  • links: An array containing routes/links/paths to other Locations.

  • numAgents: A tracking variable that keeps count as to how many people are present at this Location.

Rules for movement and state changes

Great! You’ve established the Locations in your simulation, as well as the Links interconnecting them, as well as Persons that can reside either in Locations or Links. Now each Person will have to make decisions at different moments to determine her/his behavior. In this code we model two types of decisions:

  1. Whether the Person wishes to move from its current location to another one.

  2. If 1 is the case: which route the Person will choose from a set of routes.

We will start with decision 2, which is at the lowest level, and create a simple function that picks a favourite route amongst a list of routes. To do this, we created a simple weighted choice algorithm:

def selectRoute(self):
  total_score = 0.0
  for i in range(0,len(self.location.links)):
    total_score += 40000.0 / (10.0 + self.location.links[i].distance)

  selected_value = random.random() * total_score

  checked_score = 0.0
  for i in range(0,len(self.location.links)):
    checked_score += 40000.0 / (10.0 + self.location.links[i].distance)
    if selected_value < checked_score:
      return i

Here, each option has a weight equal to 40000 (the approximate circumference of the planet in km) divided by (10 + [distance to the endpoint of the route in km]).

Because the function is rather simple, I included a full implementation. However, the exact same functionality can also be accomplished using numpy.random.choice(), if you have access to the numpy library.

selectRoute() is embedded in a more general function (evolve()), which evolves the position of a Person over a single timestep in the simulation. This function essentially captures the mechanics in making decision 1, and relies on the aforementioned selectRoute() to resolve decision 2 when necessary:

def evolve(self):

  if not self.travelling:
    movechance = self.location.movechance
    outcome = random.random()

    if outcome < movechance:
      # determine here which route to take?
      chosenRoute = self.selectRoute()

      # update location to link endpoint
      self.location.numAgents -= 1
      self.location = self.location.links[chosenRoute]
      self.location.numAgents += 1
      self.travelling = True

Here the chance of a Person moving at all at a given time step is given by the movechance. This movechance is a static number for each Location, allowing us to set a high movechance for unsafe locations, and a lower movechance for safer locations.

evolve() places Persons on the Links. To ensure that these Persons reach there destination we create one more function, namely finish_travel()

def finish_travel(self):
  # if the person resides on a link between locations, it is "travelling"
  if self.travelling:

    # increment the distance covered by 10 kilometers.
    self.distance_travelled_on_link += 10

    # get the length of the current route (link).
    link_length = self.location.distance

    # If the distance travelled is longer than the length of the link, we arrive at our destination.
    if self.distance_travelled_on_link > link_length:
      self.location.numAgents -= 1
      self.location = self.location.endpoint
      self.location.numAgents += 1
      self.travelling = False
      self.distance_travelled_on_link = 0

  # Update the X and Y coordinates of each agent
  if self.travelling:
      self.x = self.location.calc_x(self.distance_travelled_on_link)
      self.y = self.location.calc_y(self.distance_travelled_on_link)
  else:
    self.x = self.location.x
    self.y = self.location.y

This function allows us to track agents who are on links, and have them progress gradually.

From state to simulation

We now have people, locations, and links that represent connections between these locations. These are essential components for an agent-based model in this context. It’s easy to think up many other possible components (e.g., conflict events, other types of agents, more parameters regarding age, religion etc.), but most of these are not essential for the simulation in its most basic form. However, what is essential is to be able to model a period of time, i.e. turning out frozen state into a simulation.

To accomplish this, we create an Ecosystem class, which stores the full state (Locations, Links and Persons), and which is able to evolve them in time. We define the class as follows:

class Ecosystem:
  def __init__(self):
    self.locations = []
    self.locationNames = []
    self.agents = []
    self.time = 0

The Ecosystem class has the following attributes:

  • locations: Contains all the locations in our system.

  • locationNames: A shorthand list of the names of the respective locations in our system, to make it easier to write diagnostic information.

  • agents: A list of all the agents in our system.

  • time: Basically a clock, which contains the number of time steps that have been taken.

Next, we need a member function that adds locations to the Ecosystem:

def addLocation(self, name, x="0.0", y="0.0", movechance=0.1):
  l = Location(name, x, y, movechance)
  self.locations.append(l)
  self.locationNames.append(l.name)
  return l

…a function that adds Agents to the Ecosystem:

def addAgent(self, location):
  self.agents.append(Person(location))

…and a function that adds Links to the Ecosystem:

def linkUp(self, endpoint1, endpoint2, distance="1.0"):
  """ Creates a link between two endpoint locations
  """
  endpoint1_index = 0
  endpoint2_index = 0
  for i in range(0, len(self.locationNames)):
    if(self.locationNames[i] == endpoint1):
      endpoint1_index = i
    if(self.locationNames[i] == endpoint2):
      endpoint2_index = i

  self.locations[endpoint1_index].links.append( Link(self.locations[endpoint1_index], self.locations[endpoint2_index], distance) )
  self.locations[endpoint2_index].links.append( Link(self.locations[endpoint2_index], self.locations[endpoint1_index], distance) )

Crucially, we want to evolve the system in time. This is actually done using the following function:

def doTimeStep(self):
  #update agent locations
  for a in self.agents:
    a.evolve()

  #update agent travel on links
  for a in self.agents:
    a.finish_travel()

  self.time += 1

Lastly, we add two functions to aid us in writing out some results.

def numAgents(self):
  return len(self.agents)

def printLocationInfo(self):
  my_file = open("locations.csv", "w")
  my_file.write("#name,x,y\n")
  for l in self.locations:
    my_file.write("%s,%s,%s\n" % (l.name, l.x, l.y))
  my_file.close()

def printInfo(self):
  print("Time: ", self.time, ", # of agents: ", len(self.agents))
  for l in self.locations:
    print(l.name, l.numAgents)

  my_file = open("agents.%s.csv" % (self.time), "w")

  my_file.write("#id,x,y\n")
  for id,a in enumerate(self.agents):
    my_file.write("%s,%s,%s\n" % (id, a.x, a.y))
  my_file.close()

Creating and running an Agent-based Simulation

We have now created all the essential classes to perform an agent-based simulation. Here we describe how you can construct and run a simple ABM simulation. We start off by creating an Ecosystem, and creating a location graph with six locations in it. The location graph will roughly look like this:

_images/locations.png

And the source code required to add the locations for this involves:

if __name__ == "__main__":
  print("A first ABM implementation")

  e = Ecosystem()

  l1 = e.addLocation("Source1",x=200,y=0)
  l2 = e.addLocation("Source2",x=100,y=100)
  l3 = e.addLocation("Transit1",x=100,y=0)
  l4 = e.addLocation("Transit2",x=200,y=100)
  l5 = e.addLocation("Sink1",x=300,y=0)
  l6 = e.addLocation("Sink2",x=0,y=100)

Next, we establish two paths, each of which connects the source location to one of the two sink locations. As a test, we specify one of the paths to have a length of 10 kilometers, and one to have a length of 5 kilometers:

e.linkUp("Source1","Transit1","100.0")
e.linkUp("Source1","Transit2","50.0")
e.linkUp("Source2","Transit1","100.0")
e.linkUp("Source2","Transit2","50.0")
e.linkUp("Transit1","Sink1","200.0")
e.linkUp("Transit2","Sink2","200.0")

With the location and links in place, we can now insert a hundred agents in the source location l1. To do that, we use the addAgent() function a hundred times.

for i in range(0,100):
  e.addAgent(location=l1)

With all the agents in place, we can now proceed to run the simulation. We first print all the locations to a CSV file for later reference. Next, we run the simulation for a duration of 10 time steps, and we print basic diagnostic information after each time step:

e.printLocationInfo()

duration=10
for t in range(0,duration):
  e.doTimeStep()
  e.printInfo()

…and with that all in place, you have just established your first working ABM model!

You can run your simulation using: python3 <name_of_the_python_script_in_which_you_stored_your_code>

If it runs successfully, it will create 11 CSV files in the directory that you launch it from. These files include locations.csv as well as 10 agent log files, named agents.1.csv all the way to agents.10.csv.

Optional: Visualizing your results

In this section we will try to create a visualization of your simulation, so you can explore how the agents move around the location graph.

installing dependencies

For this section, you will need the Python3 matplotlib and pandas packages. In addition, to save the visualizations, you will need to install imagemagick. One way to install these on Linux platforms is by using the following commands:

pip3 install pandas, pip3 install matplotlib and sudo apt install imagemagick.

Main visualization

To show an animation of your results, you can paste the following code into a file named make_animation.py.

import numpy as np
import glob
import matplotlib.pyplot as plt
import sys
import pandas as pd
from matplotlib.animation import FuncAnimation

data_path = "."

def plot_location():
  # sample data in data directory
  location_df = pd.read_csv('%s/locations.csv' % data_path, index_col="#name")
  city_names = location_df.index.tolist()
  x = location_df.x
  y = location_df.y
  plt.scatter(x, y, s=300, alpha=0.5)
  # label the points with the city names
  for i, txt in enumerate(city_names):
    plt.annotate(txt, (x[i], y[i]), xytext=(5, 5), textcoords='offset points', fontsize='12')


def read_csv_to_df():
  # Reads data from data directory
  df_list = []

  num_files = len(glob.glob('%s/agents.*.csv' % data_path))
  for i in range(1,num_files+1):
    file_path = '%s/agents.%s.csv' % (data_path, i)
    print(file_path)
    dataframe = pd.read_csv(file_path, index_col='#id')
    dataframe.apply(pd.to_numeric)
    df_list.append(dataframe)
  return df_list


def animate(i, df, scat):
  scat.set_offsets(np.c_[df[i].x, df[i].y])
  return scat


def save_animation(anim):
  """
  Requires your host system to have the "ffmpeg" pacakage installed
  For mac use home brew: brew install ffmpeg
  this will install a lot of other dependenies required as well
  """
  # Assumes output directory exists
  anim.save('%s/agent_location.gif' % data_path, writer='imagemagick')
  print('Animation saved in output directory')


def main():

  if len(sys.argv)>1:
    global data_path
    data_path = sys.argv[1]

  fig, ax = plt.subplots(figsize=(5, 3))
  #ax.set(xlim=(-10, 110), ylim=(-10, 110))

  num_files = len(glob.glob('%s/agents*.csv' % data_path))
  scat = ax.scatter([], [])
  plot_location()

  print("# of frames: ",num_files)

  dataframe_list = read_csv_to_df()
  # time between frames can be changed by adjusting the interval param which is in milliseconds
  anim = FuncAnimation(
    fig, animate, interval=1000, frames=range(num_files), fargs=(dataframe_list, scat))

  plt.draw()
  # shows the output on screen
  plt.show()
  # uncomment line below to save as mp4 video file
  # save_animation(anim)


if __name__ == "__main__":
  main()

Once you have done so, you can then create an animation on your screen using the command: python3 make_animation.py <name_of_directory_with_output_files>