Discrete Event Simulation Guidance

I recognize the complexity of homework 8, especially for students without much programming experience. Therefore, I am going to write here a step by step guidance on how to approach this problem. Remember, our goal is to simulate a \(M/E2/1/\infty\) queue, that is, Poisson arrival - Two-stage Erlang service time - one server - infinite queue capacity.

Feel free to ask questions following this post. By providing this guidance, I am in the hope that all students, even those without programming experience, could solve this problem, at least partially. I will provide snippets of code to clarify my points, those are likely C code, as people in this department (ECE) should all know it. Having said that, I would recommend you to use a scripting language like Python or Matlab if you know them. They will save you from a lot of coding effort. This is an on-going document and I will update it as I find time to.

Random Deviates

Uniform Random Deviates Generation

The very first step is to have a way to generate independent random variates from uniformly distributed random variable [0, 1]. It turns out to be a very difficult task in a computer where everything is deterministic. Instead, we can use some complex algorithms to generate pseudorandom numbers that look random, but actually not. It is not truly random because the generated sequence is predetermined by an initial setting, called the seed. We do not expect you to implement one yourself. Luckily, almost all languages have standard libraries that provide decent random number generators. In C, for example, the rand function will generate a random integer between [0, RAND_MAX]. Then, random deviates of [0, 1] uniform a random variable can be generated by

#include <stdlib.h>

double uniform()
{
    return (double)rand() / (double)RAND_MAX ;
}

If you are using Python, then simply

import random
def uniform():
    return random.uniform(0,1)

I am sure Matlab has similar functions. Just search it.

Exponential Random Deviates Generation

A uniform random variate can be transformed into an exponential one using a method called inverse transform sampling, which you should investigate yourself. Your goal is to write a function exponential that generates exponential random deviates. For C,

double exponential_trans(double uniform_variate, double rate)
{
    //inverse transform
}
double exponential(double rate)
{
    return exponential_trans(uniform(), rate)
}

If you are using Python, do not use the built-in exponential generator.

Two-stage Erlang Random Deviates Generation

Not all types of variates can be easily generated through inverse transform sampling, especially when CDF is very complex. However, if you notice the relation between Erlang and exponential, you will find out how to generate Erlang random variates from exponential ones.

Poisson Arrival Process

For Poisson arrival processes \(N(t)\), you need to generate a sequence of arrival times, at which \(N(t)\) is increased by 1. Again, Poisson arrival process is heavily related to exponential distribution. With the ability to generate exponential random variates, you can easily generate a sequence of Poisson arrival times.

Discrete Event Simulator

Discrete event simulation is most suitable for systems with discrete state, for example, our queuing system. A system with discrete state only changes its state at discrete time points, which means nothing is happening between these time points. Therefore, during simulation, our simulator can jump directly from one of these time points to the next. These state changes are caused by events (job arrival, job departure, …). An event can be described by its type, the time point at which it happens, and optionally some extra information. In C, we can use a struct to represent events:

enum EventType
{
    JobArrival,
    JobDeparture
}

struct Event
{
    EventType type;
    double time;
    //possibly more fields
}

The state of a queuing system is not a very interesting one: the simulation time and the number of jobs/customers in the queue:

struct State
{
    double t;
    int num_jobs;
};

Now, try to simulate this system with a pencil and a piece of paper. Keep track of how state changes due to events and also think about how to compute different output measures (waiting time, queue length, etc..).

While you simulate it on paper, you will find that you need a data structure to store future events. You will also need to retrieve the nearest event from that data structure (and possibly remove events from that data structure). That data structure is called priority queue. In most cases, heap is a good choice for implementing a priority queue. However, if you do not know what a heap is, you can also use a sorted array. Be warned that a sorted array will be much slower than a heap but will nonetheless work for this problem.

The main loop of your simulator should look like this:

1. take the nearest event (called E) from the event queue.
2. advance the simulator time to the time point of E.
3. change system state depending on the type of E.
4. possibly generate new future events and put them into the event queue.