• Article

March 2016

Volume 31 Number 3

[C#]

Discrete Event Simulation: A Population Growth Example

By Arnaldo Perez

Throughout history, the ability to simulate has aided the development of multiple sciences. Medical models simulating the human body enhance the study of human anatomy. Computer simulation games such as “World of Warcraft” recreate an entire fantasy world and “Flight Simulator” helps train pilots on the ground. Various simulation programs explore responses to terrorist attacks, pandemic diseases and other possible crises. Even the simulated dinosaurs in the film “Jurassic Park” hint at the broad application of simulation and its potential.

Simulation is a technique in which a real-life system or process is emulated by a designed model. The model encapsulates all of the system’s features and behaviors; the simulation is the execution of this system over time. There are several stages to designing a simulation:

  • Defining the system to be modeled, which involves studying the problem at hand, identifying the properties of the environment and specifying the goals to reach.
  • Formulating the model, which includes defining all of its variables and their logical relations and creating the necessary flow diagrams.
  • Defining the data the model will require to produce the desired outcome.
  • Producing a computerized implementation of the model.
  • Verifying whether the implemented model satisfies the design.
  • Validating through comparison that the simulator actually represents the real system being simulated.
  • Experimenting to generate desired data using the simulator.
  • Analyzing and interpreting results from the simulator and making decisions based on these results.
  • Documenting the model created and the simulator as a tool.

Simulations generally comprise either a continuous process or discrete events. To simulate a weather system, for example, the tracking occurs continuously as all elements are constantly changing. Hence, the temperature variable placed against the time variable would be represented by a continuous curve. In contrast, airplane takeoffs or landings occur as points in time and, therefore, a simulation can consider only those precise moments or events and discard everything else. This type of simulation is known as discrete event simulation (DES), and it’s what I’ll discuss in this article.

Discrete Event Simulation

DES models a system or process as an ordered sequence of individual events over time, that is, from the time of one event to the time of the next event. Hence, in a DES simulation, time is usually much shorter than real time.

When developing a DES, there are six main elements to consider:

Objects represent elements of the real system. They have properties, they relate to events, they consume resources, and they enter and leave queues over time. In the airplane takeoff and landing scenario mentioned earlier, the objects would be airplanes. In a health care system, objects might be patients or organs. In a warehouse system, the objects would be the products in stock. Objects are supposed to interact with each other or with the system and they can be created at any time during a simulation.

Properties are features particular to each object (size, landing time, sex, price and so on) that are stored in order to determine responses to various scenarios that might occur in the simulation; such values can be modified.

Events are things that can occur in the system, especially to objects, such as the landing of an airplane, the arrival of a product at a warehouse, the appearance of a particular disease and so forth. Events can occur and reoccur in any order.

Resources are elements that provide services to objects (for example, a landing strip in an airport, storage cells in a warehouse and doctors in a clinic). When a resource is occupied and an object needs it, the object must queue and wait until the resource is available.

Queues are the conduits in which objects are organized to await the release of a resource that’s currently occupied. Queues can have a maximum capacity and they can have different calling approaches: first-in-first-out (FIFO), last-in-first-out (LIFO), or based on some criteria or priority (disease progression, fuel consumption and the like).

Time (as it happens in real life) is essential in simulation. To measure time, a clock is started at the beginning of a simulation and can then be used to track particular periods of time (departure or arrival time, transportation time, time spent with certain symptoms, and so on). Such tracking is fundamental because it allows you to know when the next event should occur.

Because simulation programming can be complicated, there have been many attempts to create languages that embody all the requirements of the simulation paradigm in order to ease development. One such language is SIMULA, invented in the 1960s by Ole Johan and Kristen Nygaard and the first to introduce the concept of object-oriented programming (OOP), today’s leading programming paradigm. Nowadays, the focus is more on creating packages, frameworks or libraries that incorporate what programmers need when creating a simulation. These libraries are meant to be called from ordinary languages like C#, C++, Java or Python.

In “A History of Discrete Event Simulation Programming Languages,” Nance proposed six minimum requirements any DES programming language should satisfy (bit.ly/1ZnvkVn):

  • Random number generation.
  • Process transformers, to allow variables other than the uniform random variables.
  • List processing, to facilitate the creation, manipulation and deletion of objects.
  • Statistical analysis, to provide a descriptive summary of model behavior.
  • Report generation, to assist in the presentation of large sets of data and facilitate decision making.
  • A time-flow mechanism.

These requirements can all be satisfied in C#, so in this article I’ll present a discrete event simulation for population growth developed in C#.

DES for Population Growth

Population growth is one of many aspects considered in the study of how populations (animals and plants) change over time and space and interact with their environment. Populations are groups of organisms of the same species living at the same time, sometimes in the same space or further delimited by particular characteristics.

Why is it important to study population growth? A better understanding of how populations grow or shrink gives scientists the possibility of making better predictions about changes in biodiversity conservation, resource use, climate change, pollution, health care, transportation needs, and so on. It also provides insight into how organisms interact with each other and the environment, a critical aspect when considering whether a population might prosper or decline.

In this article I’ll present a DES for the growth of a population. The goal is to observe how the population evolves over time and get some statistics by analyzing the end results (population size, old people, young people and so on). The population will start with m male and n female individuals, each having an associated age. Clearly, m and n must be greater than zero or the simulation would be pointless. The objects in this model are individuals.

An individual is capable of starting a relationship with another individual after reaching an age that distributes by a Poisson function with parameter λ = 18 years. (You don’t need to fully understand Poisson, normal or exponential distributions right now; they’ll be explained in the next section.)

There’s a less than a 50 percent probability of single and capable opposite-sex individuals engaging with one another, and even that occurs only when the age difference is no greater than 5 years. Figure 1 shows an example of the probability of a relationship between two individuals ending.

Figure 1 Probability of a Relationship Ending

Average Age Probability
14-20 0.7
21-28 0.5
29+ 0.2

Individuals involved in a relationship can have a child after a time that distributes by an exponential function with parameter λ= 8 years.

A woman can get pregnant if she has an age that follows a normal (bell-shaped) distribution function with parameters µ = 28 and σ2 = 8 years. Every woman has a number of children that distributes by a normal function with parameters µ = 2and σ2 = 6 years. (The parameter µ represents the average age while σ2 is the measure of age variability.)

Every individual has a life expectancy that distributes by a Poisson function with parameter λ= 70 years, where λ represents the average life expectancy.

In the previous description it’s possible to identify several events:

  • Starting a relationship.
  • Ending a relationship.
  • Getting pregnant.
  • Having a child.
  • Dying.

Every event is accompanied by a discrete random variable that determines the moment in which the event will occur.

Probabilistic Distributions

A discrete random variable is one whose set of values is finite or countably infinite. That is, the values can be listed as a finite or infinite sequence of values 1, 2, 3 …. For a discrete random variable, its probability distribution is any graph, table or formula that assigns a probability to each possible value. The sum of all probabilities must be 1, and each individual probability must be between 0 and 1. For instance, when throwing a fair die (all sides equally probable), the discrete random variable X representing the possible outcomes will have the probabilistic distribution X(1) = 1/6, X(2) = 1/6, …, X(6) = 1/6. All sides are equally probable, so the assigned probability of each is 1/6.

Parameters l and µ indicate the mean (expected value) in their corresponding distributions. The mean represents the value that the random variable takes on average. In other words, it’s the sum E=[(each possible outcome) × (probability of that outcome)], where E denotes the mean. In the case of the die, the mean would be E = 1/6 + 2/6 + 3/6 + 4/6 + 5/6 + 6/6 = 3.5. Note that the result 3.5 is actually halfway between all possible values the die can take; it’s the expected value when the die is rolled a large number of times.

Parameter σ2 indicates the variance of the distribution. Variance represents the dispersion of possible values of the random variable, and it’s always non-negative. Small variances (close to 0) tend to indicate that values are close to each other and the mean; high variances (close to 1) indicate great distance among values and from the mean.

Poisson is a discrete distribution that expresses probabilities regarding the number of events per time unit (see Figure 2). It’s usually applied when the probability of an event is small and the number of opportunities for it to occur is large. The number of misprints in a book, customers arriving at a business center, cars arriving at traffic lights and deaths per year in a given age group are all examples of applications of the Poisson distribution.

Poisson Distribution, Parameter λ = 18
Figure 2 Poisson Distribution, Parameter λ = 18

An exponential distribution expresses time between events in a Poisson process (see Figure 3). For example, if you’re dealing with a Poisson process describing the number of customers arriving at a business center during a certain time, you may be interested in a random variable that would indicate how much time passed before the first customer arrived. An exponential distribution can serve this purpose. It could also be applied to physics processes, for example to represent the lifetime of particles, where λ would indicates the rate at which the particle ages.

Exponential Distribution, Parameter λ = 8
Figure 3 Exponential Distribution, Parameter λ = 8

The normal distribution describes probabilities that tend to fall around a central value, no bias left or right, as shown in Figure 4. Normal distributions are symmetric and possess bell-shaped density curves with a single peak at the mean. Fifty percent of the distribution lies to the left of the mean and 50 percent to the right. The standard deviation indicates the spread or girth of the bell curve; the smaller the standard deviation the more concentrated the data. Both the mean and the standard deviation must be indicated as parameters to the normal distribution. Many natural phenomena closely follow a normal distribution: blood pressure, people’s height, errors in measurements and many more.

NormalDistribution, Parameters µ = 28 and σ2 = 8 Years
Figure 4 NormalDistribution, Parameters µ = 28 and σ2 = 8 Years

Now I’ll show you how to implement the proposed DES in a popular, sophisticated and elegant language such as C#.

 

Implementation

To develop this simulation I’ll exploit all the benefits of the OOP paradigm. The idea is to obtain the most readable code possible. Scientific applications tend to be complex and difficult to understand, and it’s essential to try to make them as clear as possible so that others can comprehend your code. I’ll develop the application as a console application in C#, with the structure shown in Figure 5.

The Structure of the Simulation Application
Figure 5 The Structure of the Simulation Application

A logical path is essential; note how the namespace structure maintains its own common sense: Simulation.DES.PopulationGrowth.

The Events folder contains an Event.cs file, which defines an enum type representing all possible events in the simulation:

namespace DES.PopulationGrowth.Events
{
public enum Event
  {
Capable Engaging,
Birth EngageDisengage,
GetPregnant,
ChildrenCount,
TimeChildren,
    Die
  }
}

The Objects folder contains every class related to individuals; these are the pieces of code that will benefit most from OOP. There are three classes related to individuals: Individual, Male and Female. The first is an abstract class and the others inherit from it; that is, males and females are individuals. Most of the coding happens in the Individual class. Figure 6 shows its properties and methods.

Properties and Methods of the Individual Class
Figure 6 Properties and Methods of the Individual Class

Here’s a description of every property:

  • Age: Individual age.
  • Couple: Null if individual is single, otherwise his/her couple, another Individual.
  • Engaged: True if individual is involved in a relation, false otherwise.
  • LifeTime: Age at which the individual dies.
  • RelationAge: Age at which the individual can start a relationship.
  • TimeChildren: Time in the simulation (not age) at which the individual can have children.

And here’s a description of their methods:

  • Disengage: Ends the relation between two individuals.
  • EndRelation: Determines whether the relation should end, according to the probability table in Figure 1.
  • FindPartner: Finds an available partner for an individual.
  • SuitablePartner: Determines whether an individual is suitable for starting a relationship (age difference and opposite sex).
  • SuitableRelation: Determines whether some individual can start a relationship; that is, is single and at an age for starting a relation.
  • ToString: Override to represent individual information as a string.

The class starts by defining all properties and later the constructor, which merely sets the age:

public abstract class Individual
  {
public int Age { get; set; }
public int RelationAge{ get; set; }
public int LifeTime{ get; set; }
public double TimeChildren{ get; set; }
public Individual Couple { get; set; }
protected Individual(int age)
{
  Age = age;
}

The SuitableRelation and SuitablePartner methods and the Engaged property are just logical tests, pretty straightforward:

public bool SuitableRelation()
{
return Age >= RelationAge&& Couple == null;
}
public bool SuitablePartner(Individual individual)
{
return ((individual is Male && this is Female) ||
(individual is Female && this is Male)) &&Math.Abs(individual.Age - Age) <= 5;
}  
public bool Engaged
{
get { return Couple != null; }
}

The Disengage method breaks up a relationship by setting the Couple field to null on the couple and later on the individual. It also sets their time for having children to 0 because they’re no longer engaged.

public void Disengage()
{
Couple.Couple = null;
Couple = null;
TimeChildren = 0;
}

The EndRelation method is basically a translation of the probability table for determining the chance of a relationship ending. It uses a uniform random variable, which produces a random value in the range [0, 1] equivalent to producing an acceptance percentage. The distributions dictionary is created in the simulation class (described shortly) and holds pairs (event, probability distribution), thus associating every event with its distribution:

public bool EndRelation(Dictionary<Event, IDistribution> distributions)
{
var sample =
  ((ContinuousUniform) distributions[Event.BirthEngageDisengage]).Sample();
if (Age >= 14 && Age <= 20 && sample <= 0.7)
return true;
if (Age >= 21 && Age <= 28 && sample <= 0.5)
return true;
if (Age >= 29 && sample <= 0.2)
return true;
return false;
}

The FindPartner method, shown in Figure 7, finds an available partner for the instance individual. It receives as input the population list, the current time of the simulation and the distributions dictionary.

Figure 7 The FindPartner Method

public void FindPartner(IEnumerable<Individual> population, int currentTime,
  Dictionary<Event, IDistribution> distributions)
{
foreach (var candidate in population)
if (SuitablePartner(candidate) &&
candidate.SuitableRelation() &&
((ContinuousUniform) distributions[Event.BirthEngageDisengage]).Sample() <= 0.5)
{
// Relate them
candidate.Couple = this;
Couple = candidate;
// Set time for having child
var childTime = ((Exponential)  distributions[Event.TimeChildren]).Sample()*100;
// They can have children on the simulated year: 'currentTime + childTime'.
candidate.TimeChildren = currentTime + childTime;
TimeChildren = currentTime + childTime;
break;
}
}

Last, the ToString method defines the string representation of an individual:

public override string ToString()
{
Return string.Format("Age: {0} Lifetime {1}", Age, LifeTime);
}

The Male class is very simple:

public class Male: Individual
{
public Male(int age) : base(age)
{
}
public override string ToString()
{
return base.ToString() + " Male";
}
}

The Female class is a bit more complicated because it includes treatment for pregnancy, birth and so forth. The class definition starts by declaring all properties and the constructor:

public class Female :Individual
{
public bool IsPregnant{ get; set; }
public double PregnancyAge{ get; set; }
public double ChildrenCount{ get; set; }
public Female(int age) : base(age)
{
}
}

Here are the properties of the Female class:

  • IsPregnant: Determines whether this woman is pregnant.
  • PregnancyAge: Determines age at which a woman can get pregnant.
  • ChildrenCount: Indicates number of children to give birth to.

And here are the methods it contains:

  • SuitablePregnancy: Determines whether this woman can get pregnant.
  • GiveBirth: Indicates a woman giving birth, adding new individual to the population.
  • ToString:override: Used to represent female information as a string.

SuitablePregnancy is a test method that outputs true when the instance woman satisfies every condition for getting pregnant:

public bool SuitablePregnancy(intcurrentTime)
  {
return Age >= PregnancyAge && currentTime <= TimeChildren && ChildrenCount > 0;
}

The GiveBirth method, shown in Figure 8, is the code that adds and initializes new individuals into the population.

Figure 8 TheGiveBirth Method

public Individual GiveBirth(Dictionary<Event, IDistribution> distributions, int currentTime)
  {
var sample =
  ((ContinuousUniform) distributions[Event.BirthEngageDisengage]).Sample();
var child = sample > 0.5 ? (Individual) newMale(0) : newFemale(0);
// One less child to give birth to
ChildrenCount--;
child.LifeTime = ((Poisson)distributions[Event.Die]).Sample();
child.RelationAge = ((Poisson)distributions[Event.CapableEngaging]).Sample();
if (child is Female)
{
(child as Female).PregnancyAge =
  ((Normal)distributions[Event.GetPregnant]).Sample();
(child as Female).ChildrenCount =
  ((Normal)distributions[Event.ChildrenCount]).Sample();
}
if (Engaged && ChildrenCount> 0)
{
TimeChildren =
  currentTime + ((Exponential)distributions[Event.TimeChildren]).Sample();
Couple.TimeChildren = TimeChildren;
}
else
TimeChildren = 0;
IsPregnant = false;
return child;
  }

A uniform sample is generated to first determine if the new individual will be female or male, each with a probability of 50 percent (0.5). The ChildrenCount value is decremented by 1, indicating that this woman has one less child to give birth to. The rest of the code relates to individual initialization and resetting some variables.

The ToString method changes the manner in which Females are represented as strings:

public override string ToString()
  {
return base.ToString() + " Female";
}

Because all functions relating to individuals were placed in separate classes, the Simulation class is now much simpler. Properties and the constructor are at the beginning of the code:

public class Simulation
  {
public List<Individual> Population { get; set; }
public int Time { get; set; }
private int _currentTime;
private readonly Dictionary<Event, IDistribution> _distributions ;

The properties and variables are self-descriptive and some were previously described. Time represents the time (in years) the simulation will last, and _currentTime represents the current year of the simulation. The constructor in this case, shown in Figure 9, is more complicated because it includes the initialization of random variables for each individual.

Figure 9 The Simulation Class Constructor

public Simulation(IEnumerable<Individual> population, int time)
{
Population = new List<Individual>(population);
Time = time;
_distributions = new Dictionary<Event, IDistribution>
{
{ Event.CapableEngaging, new Poisson(18) },
{ Event.BirthEngageDisengage, new ContinuousUniform() },
{ Event.GetPregnant, new Normal(28, 8) },
{ Event.ChildrenCount, new Normal(2, 6) },
{ Event.TimeChildren, new Exponential(8) },
{ Event.Die, new Poisson(70) },
                                 };
foreach (var individual in Population)
{
// LifeTime
individual.LifeTime = ((Poisson) _distributions[Event.Die]).Sample();
// Ready to start having relations
individual.RelationAge = ((Poisson)_distributions[Event.CapableEngaging]).Sample();
// Pregnancy Age (only women)
if (individual is Female)
                {
(individual as Female).PregnancyAge = ((Normal) _distributions[Event.GetPregnant]).Sample();
(individual as Female).ChildrenCount = ((Normal)_distributions[Event.ChildrenCount]).Sample();
}
}
}

Finally, the Execute method, shown in Figure 10, is where all simulation logic occurs.

Figure 10 The Execute Method

public void Execute()
{
while (_currentTime< Time)
{
// Check what happens to every individual this year
for (vari = 0; i<Population.Count; i++)
{
var individual = Population[i];
// Event -> Birth
if (individual is Female&& (individual as Female).IsPregnant)
Population.Add((individual as Female).GiveBirth(_distributions,
  _currentTime));
// Event -> Check whether someone starts a relationship this year
if (individual.SuitableRelation())
individual.FindPartner(Population, _currentTime, _distributions);
// Events where having an engaged individual represents a prerequisite
if (individual.Engaged)
{
// Event -> Check whether arelationship ends this year
if (individual.EndRelation(_distributions))
individual.Disengage();
// Event -> Check whether a couple can have a child now
if (individual is Female &&
(individual as Female).SuitablePregnancy(_currentTime))
(individual as Female).IsPregnant = true;
}
// Event -> Check whether someone dies this year
if (individual.Age.Equals(individual.LifeTime))
{
// Case: Individual in relationship (break relation)
if (individual.Engaged)
individual.Disengage();
Population.RemoveAt(i);
}
individual.Age++;
_currentTime++;
}
}

Iterations in the outer loop represent years that go by in the simulation. The inner loop goes through events that might occur to those individuals in that particular year.

To see how the population evolves over time, I set up a new console application, shown in Figure 11.

Figure 11 Viewing the Population over Time

static void main()
{
var population = new List<Individual>
{
new Male(2),
new Female(2),
new Male(3),
new Female(4),
new Male(5),
new Female(3)
};
var sim = new Simulation(population, 1000);
sim.Execute();
// Print population after simulation
foreach (var individual in sim.Population)
Console.WriteLine("Individual {0}", individual);
Console.ReadLine();
}

Figure 12 shows the population after 1,000 years.

Population After 1,000 Years
Figure 12 Population After 1,000 Years

In this article I developed a discrete event simulation to see how a population evolves in time. The object-oriented approach proved very useful for obtaining readable, concise code that readers can try and improve if necessary.

Arnaldo Pérez Castaño is a computer scientist based in Cuba. He is the author of a series of programming books—“JavaScript Fácil,” “HTML y CSS Fácil,” and “Python Fácil” (Marcombo S.A). His expertise includes Visual Basic, C#, .NET Framework, and Artificial Intelligence, and offers his services as a freelancer through nubelo.com. Cinema and music are some of his passions. Contact him at <arnaldo.skywalker@gmail.com.>

Thanks to the following Microsoft technical expert who reviewed this article: James McCaffrey