|
Dr Jon Shiach Senior Lecturer Department of Computing and Mathematics Manchester Metropolitan University Email: j.shiach@mmu.ac.uk |
|
Dr Stephen Lynch Reader Department of Computing and Mathematics Manchester Metropolitan University Email: s.lynch@mmu.ac.uk |
|
Dr Killian O'Brien Senior Lecturer Department of Computing and Mathematics Manchester Metropolitan University Email: k.m.obrien@mmu.ac.uk |
|
Mathematical modelling is a description of a situation or scenario that is described using mathematical concepts. These models can be very useful as they can allow us to identify possible solutions without first implementing them in the real world. For example, NASA used mathematical models to predict the forces experienced by astronauts to see whether they could survive a trip into space. Here we are going to look at a type of mathematical model called an Agent Based Model (ABM) which can be used to model the behaviour of a group of animals or people. Agent based models consist of a number of autonomous agents which are given a set of rules that govern the agents behaviour in the system. For example, an agent may represent a predator which is given a rule that it will chase the nearest prey. This Jupyter notebook will introduce you to agent based models and how can we write Python programs to model complex behaviour. The programs are written for you but you can make changes to the parameters or even add to the programs if you wish.
Jupyter notebooks are documents that combine text and Python code which allow readable documents such as this one to be created that contain executable code used to perform calculations. To run code in a notebook simply click on the code and then on the run button, or alternatively, press the ctrl + enter keys at the same time. Since we will be using commands to perform calculations and produce our models we need to import some commands to help us to do this. Run the code below to import the commands.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation, patches, collections
from IPython.display import HTML
from ipywidgets import IntProgress
matplotlib.rcParams['animation.embed_limit'] = 2**128
Now we can perform some calculations. The Python code below calculates the length of the hypotenuse of a right angled triangle using the lengths of the two smallest sides and prints the result. Note how the equation $c = \sqrt{a^2 + b^2}$ is entered in a similar way to how we write it on a piece of paper. Can you calculate the length of the hypotenuse when $a = 5$ and $b = 12$
a = 3
b = 4
c = np.sqrt(a ** 2 + b ** 2)
print(c)
5.0
An Agent Based Model (ABM) uses individual autonomous agents to which we can give rules that will govern their movement and behaviour. In the models we will use here the agents exist within a finite rectangular region of a given width and height. Each agent is described by its position within the region and velocity. The position is given by a set of co-ordinates $(x,y)$ where $x$ and $y$ are the horizontal and vertical distances from the bottom left-hand corner of the region. The velocity, which governs the direction the agent moves in, is given by the velocity vector (a vector is an arrow with a direction and length) which is represented by $(v_x,v_y)$ where $v_x$ and $v_y$ are the horizontal and vertical speeds of the agent. We can calculate the speed of an agent from the velocity vector using Pythagoras' theorem
$$\textsf{speed} = \sqrt{ v_x ^2 + v_y ^ 2}.$$This is known as the norm of the velocity vector and is equivalent to the length of the vector. The Python code below defines a function called norm() which calculates the norm of the vector (vx, vy).
def norm(vx, vy):
return np.sqrt(vx ** 2 + vy ** 2)
To program an agent based model in Python we start by defining an agent object which contains all of the commands that are applied to an agent. The code below defines an agent and gives it a random position and velocity as well as some other attributes which are explained by the comments following the # symbol.
class boid:
def __init__(self):
self.maxspeed = 10 # maximum speed of the agent
self.minspeed = 5 # minimum speed of the agent
self.x = np.random.uniform(0, width) # horizontal co-ordinate
self.y = np.random.uniform(0, height) # vertical co-ordinate
self.vx = self.maxspeed * np.random.uniform(-1, 1) # horizontal velocity
self.vy = self.maxspeed * np.random.uniform(-1, 1) # vertical velocity
self.neighbours = [] # list of neighbouring agents
self.status = "boid" # status of the agent
self.colour = "blue" # colour of the agent (for plotting)
self.move = True # does the agent move?
def get_agent_list():
np.random.seed(0)
agent_list = [boid() for _ in range(nagents)]
agent_types = ["boid"]
agent_colours = ["blue"]
return agent_list, agent_types, agent_colours
To model the movement of agents we divide the time which we are simulating into lots of individual time steps. If we know the position and velocity of an agent at one time step we can calculate the position at the next time step using
$$\textsf{new position} = \textsf{current position} + \textsf{velocity vector} \times \textsf{time step}.$$This is shown in the diagram below where $\mathbf{p}_1$ is the current position of the agent and $\mathbf{p}_2$ and $\mathbf{p}_3$ are the positions in the next two time steps.
In the diagram above we see that the velocity vector does not change so this agent will continue to move along that same path forever. The velocity vector will change if there is some external force applied to the agent, for example, an agent may be attracted towards a particular target. The change in the velocity is known as acceleration which will be determined by any rules that we impose on our agents.
class boid(boid):
def move_agent(self):
if self.move:
self.x += self.vx * dt
self.y += self.vy * dt
The steering rules that affect the movement of an agent will each generate a steering vector that are added together to give an acceleration vector
$$\textsf{acceleration vector} = \textsf{steering vector 1} + \textsf{steering vector 2} + \textsf{steering vector 3}$$
An example of a rule imposed on an agent may be that it will travel towards a particular target as shown in the diagram below. The steering vector pointing in the direction of the target is calculated using
$$\textsf{steering vector} = \textsf{target} - \textsf{position}.$$
Since we may be combining several steering rules we may wish to prioritise some over the others. We can do this by limiting the steering vector by a steering factor which is a number that we choose the value of in order to replicate the behaviour we want.
$$\textsf{acceleration vector} = \textsf{steering vector} \times \textsf{steering factor}.$$
The diagram below shows how an agent will move towards a target over several time steps. Note that since the acceleration vector is limited by the steering factor the agents does immediately turn in the direction of the target, rather it gradually turns over a number of time steps.
We must ensure that an agent does not pass through the edges of the region and escape. This is shown in the diagram below where the agent is approaching the bottom edge and if no external for is applied it would pass through the edge and out of the region.
To ensure agents remain in the region we apply a steering force if an agent gets too close to an edge which forces it to move away from the edge. A buffer zone is defined that contains the region within some distance of the edge, if an agent is in the buffer zone then a steering force is applied calculated using
\begin{align*} \textsf{steering vector} &= \frac{\textsf{direction vector}}{\textsf{distance from edge}}, \end{align*}where the direction vector is a vector that points towards the interior of the region. By dividing the direction vector by the distance from the edge the steering force is week when the agent is further away from the edge but becomes progressively stronger as the distance from the decreases. This means that the agent will gradually turn away from the edge.
The code below defines the function avoid_edges() which applies this rule for all four edges of the region defined the by width and height of the region that contains our agents and the size of the buffer zone which is defined as 2 here.
class boid(boid):
def avoid_edges(self):
x_steer, y_steer, buffer = 0, 0, 2
if self.x <= buffer:
x_steer += buffer / self.x
if self.x >= width - buffer:
x_steer -= buffer / (width - self.x)
if self.y <= buffer:
y_steer += buffer / self.y
if self.y >= height - buffer:
y_steer -= buffer / (height - self.y)
return x_steer, y_steer
The first model we are going to look at is the modelling of flocking behaviour that can be seen in shoals of fish. Each individual fish in the group reacts to the other fish that surround it and adjusts its movement accordingly. The model Boids which was developed by Craig Reynolds uses three simple rules that are applied to the agents in the model to simulate flocking behaviour, these rules are:
The acceleration forces for all three of these rules is calculated and added to the acceleration vector for an agent.
We want our agents to react to other agents which are the vicinity, for example, if we have an agent that represents a predator it will want to chase after any prey agent and a prey agent will want to flee away from a predator. In order to do this we need to determine which of the other agents can be detected by our agent. The simplest way to do this is to define a detection zone that is a circle surround our agent and all other agents that are within this circle are considered to have been detected by our agent. This could represent a predator using their senses of vision and smell to hunt for their prey. Agents that are within the detection zone are known as neighbouring agents because that are the agents that are nearest to our agent.
The code below defines the function Neighbours() which loops through all of the other agents in the model and calculates the squared distance to a current agent. If this distance is less than the square of the detection radius then it is added to a list of neighbours. The reason why the square of the distance is used is to avoid having to calculate a square root which is slower for a computer to calculate.
class boid(boid):
def find_neighbours(self, agent_list):
self.neighbours = []
for neighbour in agent_list:
if id(self) == id(neighbour):
continue
if norm(self.x - neighbour.x, self.y - neighbour.y) < detection_radius:
self.neighbours.append(neighbour)
The alignment rule states that an agent will align itself to the average velocity vector of its neighbouring agents. The average velocity is calculated using
$$\textsf{average velocity} = \frac{\displaystyle\sum \textsf{velocity of neighbour}}{\textsf{number of neighbours}}.$$The symbol $\sum$ is the Greek character sigma which is equivalent to the Latin character s and is used to denote the sum of a group of numbers. The equation above says that we add up all of the velocity vectors of the neighbouring agents and divide by the number of neighbours to get the average velocity. The steering vector for aligning to the neighbouring agents is then calculated by subtracting the agent velocity from the average velocity of the neighbourng agents
\begin{align*} \textsf{alignment vector} &= \textsf{average velocity} - \textsf{agent velocity}. \end{align*}The alignment vector is then limited to ensure that its norm is not greater than the alignment factor which controls the extent to which the agent aligns to its neighbours. The code below defines the function alignment() applies the alignment rule to the current agent.
class boid(boid):
def alignment(self):
n, sum_vx, sum_vy, x_alignment, y_alignment = 0, 0, 0, 0, 0
for neighbour in self.neighbours:
sum_vx += neighbour.vx
sum_vy += neighbour.vy
n += 1
if n > 0:
x_alignment = sum_vx / n - self.vx
y_alignment = sum_vy / n - self.vy
return limit(x_alignment, y_alignment, 0, alignment_factor)
The function limit() used in the alignment() function is defined below. It uses the norm() function to limit a vector v to ensure that its norm is within min_norm and max_norm.
def limit(vx, vy, min_norm, max_norm):
v_norm, scale = norm(vx, vy), 1
if v_norm > 0 and v_norm < min_norm:
scale = min_norm / (v_norm)
elif v_norm > max_norm:
scale = max_norm / v_norm
return vx * scale, vy * scale
The cohesion rule states that an agent will steer towards the average position of it's neighbouring agents. The steering vector is calculated by subtracting the agent position from the average position of the neigbouring agents
\begin{align*} \textsf{cohesion vector} &= \frac{\displaystyle\sum \textsf{position of neighbour}}{\textsf{number of neighbours}} - \textsf{agent position}. \end{align*}Like with the alignment vector, the cohesion vector is then limited so that it's norm does not exceed the cohesion factor. The code below defines the function Cohesion() applied the cohesion rule to the current agent.
class boid(boid):
def cohesion(self):
n, sum_x, sum_y, x_cohesion, y_cohesion = 0, 0, 0, 0, 0
for neighbour in self.neighbours:
sum_x += neighbour.x
sum_y += neighbour.y
n += 1
if n > 0:
x_cohesion = sum_x / n - self.x
y_cohesion = sum_y / n - self.y
return limit(x_cohesion, y_cohesion, 0, cohesion_factor)
The separation rule states that an agent will steer away from agents that are within some small distance of it. For those neighbouring agents that are too close we calculate the steering vector using
\begin{align*} \textsf{separation vector} &= \sum \dfrac{\textsf{position of agent} - \textsf{position of neighbour}}{\textsf{distance to neighbour}}. \end{align*}By dividing the vector that points away from the neighbouring agent by the distance to that neighbour we scale the separation vector so that the closer a neighbour is, the stronger the seperation force from that neighbour. Once again the separation vector is limited so that it does not exceed the separation factor.
The code below defines a function called separation() applies the separation rule for the current agent.
class boid(boid):
def separation(self):
x_separation, y_separation = 0, 0
for neighbour in self.neighbours:
d = norm(self.x - neighbour.x, self.y - neighbour.y)
if d < separation_radius:
x_separation += (self.x - neighbour.x) / d
y_separation += (self.y - neighbour.y) / d
return limit(x_separation, y_separation, 0, separation_factor)
The code below defines the function agent_rules() which applies the rules that govern the behaviour of the agent. The individual steering vectors for avoiding the edges of the region, alignment, cohesion and separation are calculated and added to the current velocity vector for the agent. The velocity vector is then limited to ensure that it's norm is between the minimum and maximum speed of the agent.
class boid(boid):
def agent_rules(self):
x_edge, y_edge = self.avoid_edges()
x_align, y_align = self.alignment()
x_cohesion, y_cohesion = self.cohesion()
x_separate, y_separate = self.separation()
self.vx += x_edge + x_align + x_cohesion + x_separate
self.vy += y_edge + y_align + y_cohesion + y_separate
self.vx, self.vy = limit(self.vx, self.vy, self.minspeed, self.maxspeed)
Now that we have defined the rules for our model we now need to apply these rules for each time step. For each time step we first calculate the rules for each agent in the model and then calculate the new positions of each agent using the acceleration vector that has been calculated when we applied the rules.
The function abm() defined initialises the positions and velocities of the agents in the model and applies the model through time.
def abm(model):
def update(n):
# Calculate the rules for each agent
for agent in agent_list:
agent.find_neighbours(agent_list)
agent.agent_rules()
# Plot agents
patches = []
colours = []
for agent in agent_list:
# Move agent
agent.move_agent()
# Plot agent as an isocelese triangle
angle = np.arctan2(agent.vy, agent.vx)
cos, sin, s = np.cos(angle), np.sin(angle), min(width, height) / 100
vertices = np.array([[-1, -1, 1], [2, 0, 1], [-1, 1, 1]])
A = np.array([[s * cos, s* sin, 0], [-s * sin, s * cos, 0], [agent.x, agent.y, 1]])
vertices = np.dot(vertices, A)
triangle = plt.Polygon(vertices[:,:2])
patches.append(triangle)
colours.append(agent.colour)
collection.set_paths(patches)
collection.set_color(colours)
progress_bar.value += 1
return collection,
# Generate agent lists
agent_list, agent_types, agent_colours = model.get_agent_list()
# Setup plot
fig, ax = plt.subplots(figsize=(8, 8 * height / width))
plt.axis([0, width, 0, height])
collection = collections.PatchCollection([plt.Polygon(np.zeros((3,2))) for _ in range(nagents)])
ax.add_collection(collection)
progress_bar = IntProgress(min=0, max=int(fps * tmax), description="Running...")
display(progress_bar)
# Create animation
anim = animation.FuncAnimation(fig, update, frames=int(tmax * fps), blit=True)
plt.close()
return anim
The following code defines the various parameters that define the model before running the model and generating the animation. Try experimenting with the model by changing these parameters and then clicking on the 'Run' button.
# Model parameters
nagents = 50 # number of agents
width, height = 40, 30 # width and height of the area
tmax = 20 # max time for the simulation
fps = 30 # frames per second for the animation
dt = 1 / fps # seconds elapsed during 1 frame of the animation
detection_radius = 5 # radius for detecting neighbouring agents
separation_radius = 1 # radius for separating agents
# Steering factors (how much does each behaviour affect the movement of the agent)
alignment_factor = 0.2
cohesion_factor = 0.2
separation_factor = 0.4
# Run simulation
anim = abm(boid)
HTML(anim.to_jshtml(fps=fps))
IntProgress(value=0, description='Running...', max=600)