Introduction

Our purpose is to show how to do a variety of network related simulations involving random variables with Python. In particular we'll show some of the great features of Python, MatPlotLib and SimPy libraries. All code has been tested with Python version 3.6 and current versions of the referenced libraries as of June 2017.

This tutorial will not cover theoretical aspects of probability, statistics, or queueing theory. For those interested in digging deeper or reviewing probability and queueing theory the following free resources are highly recommended and are listed in order of most advanced to most fundamental.

• Zukerman, "Introduction to Queueing Theory and Stochastic Teletraffic Models" http://arxiv.org/abs/1307.2968, July 2013. Advanced (suitable for a whole grad course or two)
• Grinstead & Snell "Introduction to Probability" http://www.clrn.org/search/details.cfm?elrid=8525 Junior/Senior level treatment
• Illowsky & Dean, "Collaborative Statistics" http://cnx.org/content/col10522/latest/ Web based, easy lookups, Freshman/Sophomore level

Random Variables in Python and SciPy

Python makes it extremely easy to generate samples of a random number for simulation or other purposes. The Python standard library includes a module random explicitly for this purpose.

import random

print("Here is a single sample from a uniform random variable")
print(random.random())
print("Here is a list of three samples:")
uniSamples = [random.random(), random.random(), random.random()]
print(uniSamples)
print("Here is a list of three exponential samples:")
expSamples = [random.expovariate(1.0), random.expovariate(1.0), random.expovariate(1.0)]
print(expSamples)

Code Fragment 1. Simple example of Python random number generation.

The previous code should work in any Python 2.7 release that includes the standard library. As a bit of a review, random number generators implemented via computer algorithms are not "truely" random but appear random as judged by a number of different criteria. The underlying random number generator in the Python standard library is based on the very high quality Mersenne Twister algorithm. Which is very suitable for almost all types of simulations. What it is not suitable for is any cryptographic purpose!

Visualizing Distributions: Histograms with Matplotlib

One of the best ways to look at a probability distribution is via a histogram. Matplotlib makes this very easy as the following two code snippets show.

import random
import matplotlib.pyplot as plt

uniSamples = [random.random() for i in range(100000)]
print(uniSamples[0:10])  #Take a look at the first 10
fig, axis = plt.subplots()
axis.hist(uniSamples, bins=100, normed=True)
axis.set_title("Histogram for uniform random number generator")
axis.set_xlabel("x")
axis.set_ylabel("normalized frequency of occurrence")
#fig.savefig("UniformHistogram.png")
plt.show()
    

Code Fragment 1.Generate uniform random variables and plot histogram. histogram1.py

import random
import matplotlib.pyplot as plt

normSamples = [random.normalvariate(9.0, 2.0) for i in range(100000)]
print(normSamples[0:10])  #Take a look at the first 10
fig, axis = plt.subplots()
axis.hist(normSamples, bins=100, normed=True)
axis.set_title(r"Histogram of an Normal RNG $\mu = 9$ and $\sigma = 2$")
axis.set_xlabel("x")
axis.set_ylabel("normalized frequency of occurrence")
plt.show()

Code Fragment 2.Generate Normal random variables and plot histogram. histogram2.py

Distributions and Network Simulations

Although probability distributions can play a wide variety of roles in modeling and simulating networks we will focus on just two uses here. First we will use a probability distribution to model the time between packet arrivals, the inter-arrival time. Second we will use a probability distribution to model the sizes of the packets sent.

A notion closely related to the packet inter-arrival time is the count of the number of packets received by a certain time. In the fancy language of stochastic processes we would call the count of the number of packets received a counting process. A point of potential confusion is that one can talk about packet arrivals either in terms of their inter-arrival time process or the counting process of the (total) number of packet arrivals by a given time. For example to enable analytical analysis of the buffer backlog at the output of a packet switch one utilizes a Poisson process to model the packet arrival count. A Poisson process is a counting process that has exponentially distributed inter-arrival times. In discrete event simulations we must compute the time until the next "event" hence we will always work with inter-arrival time distributions. Hence to generate Poisson process packet arrivals we will use an exponentially distributed random number generator and not a Poisson distributed random variable! Do not confuse the two notions. The difference between these two notions is a typical quiz or exam question.

Network Simulations Directly in Python

In network simulations we can have a number of tasks such as packet generation, transmission, queueing, that appear to run at the same time. In modern computers, operating system mechanisms such as processes and threads providing this illusion. For simplicity of implementation, ease of use, and scalability such mechanisms are not used for discrete event simulation (DES), instead techniques more akin to coroutines are frequently employed.

Although Python's generators provide functionality very similar to coroutines there is still a fair amount of work needed to create a discrete event simulation (DES) system. This is where SimPy, a very nice, open source, DES package comes in. SimPy is a general purpose DES package, not networking specific. I wanted my students to jump right into network simulations without learning too much Python or SimPy so I came up with a set of components, SimComponents.py, to assist them and I in creating network simulations. The key components/classes are:

Packet

This is a simple data object that represents a packet. Generally created by PacketGenerators. Key fields include: generation time, size, flow_id, packet id, source, and destination. We do not model upper layer protocols, i.e., packets don't contain a payload. The size (in bytes) field is used to determine transmission time.

PacketGenerator

A PacketGenerator simulates the sending of packets with a specified inter-arrival time distribution and a packet size distribution. One can set an initial delay and a finish time for packet generation. In addition one can set the source id and flow ids for the packets generated. The packet generators out member variable is used to connect the generator to any component with a put() member function.

PacketSink

Packet sinks record arrival time information from packets. This can take the form of raw arrival times or inter-arrival times. In addition the packet sink can record packet waiting times. Supports the put() operation.

SwitchPort

Models a first-in, first-out (FIFO) queued output port on a packet switch/router. You can set the rate of the output port and a queue size limit (in bytes). Keeps track of packets received and packets dropped.

PortMonitor

A PortMonitor is used if you want to monitor queue size over time for a SwitchPort. You need to specify a sampling distribution, i.e., a distribution that gives the time between successive samples.

More detailed documentation can be found in the SimComponents.py file.

Example 1: Two Packet Generators and a Sink

In code fragment 3 below we show basic usage of the PacketGenerator and PacketSink classes. In particular:

• On lines 3-12 we define functions that return either inter-arrival times or packet size distributions.
• On line 14 we create a SimPy Environment object that is used when creating our packet generators and sink. In addition we use this Environment object to run the simulation on line 23.
• On line 16 we create a packet sink and enable debugging which prints a summary for every packet received.
• On lines 17 and 18 we create our packet generators and on lines 20-21 we "wire" our packet generators to our packet sink. Note that we can wire as many outputs to an input as we like but to hook a single output to multiple inputs requires more work (packet duplication) which we have a separate component.

from random import expovariate
import simpy
from SimComponents import PacketGenerator, PacketSink

def constArrival():  # Constant arrival distribution for generator 1
return 1.5

def constArrival2():
return 2.0

def distSize():
return expovariate(0.01)

env = simpy.Environment()  # Create the SimPy environment
# Create the packet generators and sink
ps = PacketSink(env, debug=True)  # debugging enable for simple output
pg = PacketGenerator(env, "EE283", constArrival, distSize)
pg2 = PacketGenerator(env, "SJSU", constArrival2, distSize)
# Wire packet generators and sink together
pg.out = ps
pg2.out = ps
env.run(until=20)

Code Fragment 3.Simple packet generation and debugging. CompExample1.py

Partial example output:

id: 1, src: EE283, time: 1.5, size: 67.4340963648
id: 1, src: SJSU, time: 2.0, size: 92.0114018208
id: 2, src: EE283, time: 3.0, size: 66.6183260991
id: 2, src: SJSU, time: 4.0, size: 5.07004835759
id: 3, src: EE283, time: 4.5, size: 256.491861122
id: 3, src: SJSU, time: 6.0, size: 302.18680785
id: 4, src: EE283, time: 6.0, size: 45.2649887809
id: 5, src: EE283, time: 7.5, size: 142.46239579

Code Fragment 4.Example output from above code.

Example 2: Overloaded Switch Port

In this example, Code Fragment 5, we have a packet generator connected to a switch port which is then connected to a packet sink. We see from lines 4-5 that packets will be generated with a constant interarrival time of 1.5 seconds and from lines 7-8 that the packets will have a constant size of 100 bytes. This gives an average arrival rate of about 533.3 bps. On line 14 we create our switch port but it only has a line rate of 200 bps and a queue limit of 300 bytes. Hence should be pushed into dropping packets quickly.

We see from the output, Code Fragment 6, that only four packets are received at the sink by the time the simulation is stopped (the simulation can be resumed with just another call to the run() method). In addition to the debugging output for each packet received at the sink we have the information from the packet sink and switch port as printed in lines 20-21. Note how the first 100 byte packet incurs a delay of 4 seconds (due to the 200 bps line rate) and also the received and dropped packet counts.

import simpy
from SimComponents import PacketGenerator, PacketSink, SwitchPort

def constArrival():
return 1.5    # time interval

def constSize():
return 100.0  # bytes

env = simpy.Environment()  # Create the SimPy environment
ps = PacketSink(env, debug=True) # debug: every packet arrival is printed
pg = PacketGenerator(env, "SJSU", constArrival, constSize)
switch_port = SwitchPort(env, rate=200.0, qlimit=300)
# Wire packet generators and sinks together
pg.out = switch_port
switch_port.out = ps
env.run(until=20)
print("waits: {}".format(ps.waits))
print("received: {}, dropped {}, sent {}".format(ps.packets_rec,
switch_port.packets_drop, pg.packets_sent))

Code Fragment 5.Packet generator, switch port and packet sink. CompExample2.py

Example output:

id: 1, src: SJSU, time: 1.5, size: 100.0
id: 2, src: SJSU, time: 3.0, size: 100.0
id: 3, src: SJSU, time: 4.5, size: 100.0
id: 4, src: SJSU, time: 6.0, size: 100.0
waits: [4.0, 6.5, 9.0, 11.5]
received: 4, dropped 6, sent 13

Code Fragment 6.Output from above code.

Example 3: An M/M/1 Queueing System

Now let's see how to simulate an switching port with exponential packet inter-arrival times, and exponentially distributed packet sizes. With a FIFO queueing discipline and no limit on queue sizes such a system is equivalent to the well known M/M/1 queueing system which has nice analytical expressions for average queue size and waiting times.

The code to simulate an M/M/1 queue, Code Fragment 7, is very similar to that used in the overloaded switch port example, Code Fragment 5, however we have made the following changes:

• On lines 7-9 we have used the standard Python functools module to slightly ease the definition of functions returning a random sample with a given parameter.
• On line 14 where we create the packet sink we have turned off debugging since we will be processing a great many packets.
• On line 18 we create a Port Monitor to track the queue size.
• On line 23 we run the simulation significantly longer than before.

import random
import functools
import simpy
import matplotlib.pyplot as plt
from SimComponents import PacketGenerator, PacketSink, SwitchPort, PortMonitor

adist = functools.partial(random.expovariate, 0.5)
sdist = functools.partial(random.expovariate, 0.01)  # mean size 100 bytes
samp_dist = functools.partial(random.expovariate, 1.0)
port_rate = 1000.0

env = simpy.Environment()  # Create the SimPy environment
# Create the packet generators and sink
ps = PacketSink(env, debug=False, rec_arrivals=True)
pg = PacketGenerator(env, "Greg", adist, sdist)
switch_port = SwitchPort(env, port_rate, qlimit=10000)
# Using a PortMonitor to track queue sizes over time
pm = PortMonitor(env, switch_port, samp_dist)
# Wire packet generators, switch ports, and sinks together
pg.out = switch_port
switch_port.out = ps
# Run it
env.run(until=8000)
print("Last 10 waits: "  + ", ".join(["{:.3f}".format(x) for x in ps.waits[-10:]]))
print("Last 10 queue sizes: {}".format(pm.sizes[-10:]))
print("Last 10 sink arrival times: " + ", ".join(["{:.3f}".format(x) for x in ps.arrivals[-10:]]))
print("average wait = {:.3f}".format(sum(ps.waits)/len(ps.waits)))
print("received: {}, dropped {}, sent {}".format(switch_port.packets_rec, switch_port.packets_drop, pg.packets_sent))
print("loss rate: {}".format(float(switch_port.packets_drop)/switch_port.packets_rec))
print("average system occupancy: {:.3f}".format(float(sum(pm.sizes))/len(pm.sizes)))

Code Fragment 7. Exponential inter-arrival times and exponential packet sizes, to a fixed rate FIFO switch output port. CompMM1Example.py.

In Code Fragment 8 we show example output from a simulation run. Note since we are using random numbers your output should be different. The arrival times λ = 0.5 packets per second, and the packet service rate μ = 1.25 packets per second which gives a utilization ρ = 0.4. The theoretical mean queue size = 2/3 and the mean waiting time is = 4/3.

In addition to computing statistics we can use the measurements from the packet sink to create histograms as shown in Figure 2 below. Note the roughly exponential character of plot. This is demonstration of Burke's theorem. Which is useful in the analysis of networks of queues.

Last 10 waits: 4.337, 4.609, 4.328, 3.991, 3.454, 3.623, 2.084, 0.429, 0.817, 0.445
Last 10 queue sizes: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Last 10 sink arrival times: 0.992, 0.702, 0.318, 0.333, 0.216, 0.321, 0.228, 0.236, 6.733, 0.289
average wait = 1.365
received: 4141, dropped 0, sent 4141
loss rate: 0.0
average system occupancy: 0.698

Code Fragment 8. Example output from the previous program.

Simulating & Analyzing Networks of Queues

Explanation to be written... But we have example code to simulate the network of queues shown inf Figure 3 here: QueueNet2.py

Traffic Shaping

Explanation to be written... But we have implemented a Token bucket shaper in the SimComponents.py module and the example code that uses it to produce Figure 4, shown below, is given here: ShaperExample1.py

Queueing Disciplines

Explanation to be written... But we have implemented Virtual Clock and Weighted Fair Queueing disciplines in the SimComponents.py module and have example code that uses the WFQ component to produce Figure 5, shown below, here: WFQueueExample.py