2.4 Simulating a Queueing System By Hand

This section builds on the concepts discussed in the previous section in order to provide insights into how discrete event simulations operate. In order to do this, we will simulate a simple queueing system by hand. That is, we will process each of the events that occur as well as the state changes in order to trace the operation of the system through time.

Consider again the simple banking system described in the previous section. To simplify the situation, we assume that there is only one teller that is available to provide service to the arriving customers. Customers that arrive while the teller is already helping a customer form a single waiting line, which is handled in a first come, first served manner. Let’s assume that the bank opens at 9 am with no customers present and the teller idle. The time of arrival of the first eight customers is provided in the following Table 2.4

Table 2.4: Time of arrivals and service times for 8 customers
Customer Number Time of Arrival Service Time
1 3 4
2 11 4
3 13 4
4 14 3
5 17 2
6 19 4
7 21 3
8 27 2

We are going to process these customers in order to recreate the behavior of this system over from time 0 to 31 minutes. To do this, we need to be able to perform the “bookkeeping” that a computer simulation model would normally perform for us. Thus, we will need to define some variables associated with the system and some attributes associated with the customers. Consider the following system variables.

System Variables

  • Let \(t\) represent the current simulation clock time.

  • Let \(N(t)\) represent the number of customers in the system (bank) at any time \(t\).

  • Let \(Q(t)\) represent the number of customers waiting in line for the teller at any time \(t\).

  • Let \(B(t)\) represent whether or not the teller is busy (1) or idle (0) at any time \(t\).

Because we know the number of tellers available, we know that the following relationship holds between the variables:

\[N\left( t \right) = Q\left( t \right) + B(t)\]

Note also that, knowledge of the value of \(N(t)\) is sufficient to determine the values of \(Q(t)\) and \(B(t)\) at any time \(t.\) For example, if we know that there are 3 customers in the bank, then \(N\left( t \right) = 3\), and since the teller must be serving 1 of those customers, \(B\left( t \right) = 1\) and there must be 2 customers waiting, \(Q\left( t \right) = 2\). In this situation, since \(N\left( t \right)\), is sufficient to determine all the other system variables, we can say that \(N\left( t \right)\) is the system’s state variable. The state variable(s) are the minimal set of system variables that are necessary to represent the system’s behavior over time. We will keep track of the values of all of the system variables during our processing in order to collect statistics across time.

Attributes are properties of things that move through the system. In the parlance of simulation, we call the things that move through the system entity instances or entities. In this situation, the entities are the customers, and we can define a number of attributes for each customer. Here, customer is a type of entity or entity type. If we have other things that flow through the system, we identity the different types (or entity types). The attributes of the customer entity type are as follows. Notice that each attribute is subscripted by \(i\), which indicates the \(i^{th}\) customer instance.

Entity Attributes

  • Let \(\mathrm{ID}_{i}\) be the identity number of the customer. \(\mathrm{ID}_{i}\) is a unique number assigned to each customer that identifies the customer from other customers in the system.

  • Let \(A_{i}\) be the arrival time of the \(i^{\mathrm{th}}\) customer.

  • Let \(S_{i}\) be the time the \(i^{\mathrm{th}}\) customer started service.

  • Let \(D_{i}\) be the departure time of the \(i^{\mathrm{th}}\) customer.

  • Let \(\mathrm{ST}_{i}\) be the service time of the \(i^{\mathrm{th}}\) customer.

  • Let \(T_{i}\) be the total time spent in the system of the \(i^{\mathrm{th}}\) customer.

  • Let \(W_{i}\) be the time spent waiting in the queue for the \(i^{\mathrm{th}}\) customer.

Clearly, there are also relationships between these attributes. We can compute the total time in the system for each customer from their arrival and departure times:

\[T_{i} = D_{i} - A_{i}\]

In addition, we can compute the time spent waiting in the queue for each customer as:

\[W_{i} = T_{i} - \mathrm{ST}_{i} = S_{i} - A_{i}\]

As in the previous section, there are two types of events: arrivals and departures. Let \(E(t)\) be (A) for an arrival event at time \(t\) and be (D) for a departure event at time \(t\). As we process each event, we will keep track of when the event occurs and the type of event. We will also keep track of the state of the system after each event has been processed. To make this easier, we will keep track of the system variables and the entity attributes within a table as follows.

Table 2.5: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA

In the table, we see that the initial state of the system is empty and idle. Since there are no customers in the bank, there are no tabulated attribute values within the table. Reviewing the provided information, we see that customer 1 arrives at \(t = 3\) and has a service time of 4 minutes. Thus,\(\ \text{ID}_{1} = 1\), \(A_{1} = 3\), and \(\text{ST}_{1} = 4\). We can enter this information directly into our bookkeeping table. In addition, because the bank was empty and idle when the customer arrived, we can note that the time that customer 1 starts service is the same as the time of their arrival and that they did not spend any time waiting in the queue. Thus, \(S_{1} = 3\) and \(W_{1} = 0\). The table has been updated as follows.

Table 2.6: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)

Because customer 1 arrived to an empty system, they immediately started service at time 3. Since we know the service time of the customer, \(\text{ST}_{1} = 4\), and the current time, \(t = 3\), we can determine that customer 1, will depart from the system at time 7 (\(t = 3 + 4 = 7\)). That means we will have a departure event for customer 1 at time 7. According to the provided data, the next customer, customer 2, will arrive at time 11. Thus, we have two pending events, a departure of customer 1 at time 7 and the arrival of customer 2 at time 11. This fact is noted in the column labeled “Pending E(t)” for pending events. Here E(7) = D(1), E(11) = A(2) indicates that customer 1 with depart,\(\ D_{1},\) at time 7, \(E\left( 7 \right)\) and customer 2 will arrive, \(A_{2}\), at the event at time 11, \(E(11)\). Clearly, the next event time will be the minimum of 7 and 11, which will be the departure of the first customer. Thus, our bookkeeping table can be updated as follows.

Table 2.7: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)

Since there are no customers in the bank and only the one pending event, the next event will be the arrival of customer 2 at time 11. The table can be updated as follows.

Table 2.8: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)
11 A 1 1 0 2 11 4 11 NA NA 0 E(13) = A(3), E(15) = D(2)

Since the pending event set is E(13) = A(3), E(15) = D(2) the next event will be the arrival of the third customer at time 13 before the departure of the second customer at time 15. We will now have a queue form. Updating our bookkeeping table, yields:

Table 2.9: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)
11 A 1 1 0 2 11 4 11 NA NA 0 E(13) = A(3), E(15) = D(2)
13 A 2 1 1 3 13 4 15 NA NA 2 E(14) = A(4), E(15) = D(2)

Customer 3 will start service, when customer 2 departs. When customer 2 departs at time 15, we can go back and update row 5 (t = 13) to indicate that customer 3’s value for \(S(i)\) is now 15. With \(S(i) = 15\), we can compute the waiting time \(W(i)\) for customer 3 as \(S(i) - A(i) = 15 - 13 = 2\) for row 5. Reviewing the pending event set, we see that the next event will be the arrival of customer 4 at time 14, which yields the following bookkeeping table.

Table 2.10: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)
11 A 1 1 0 2 11 4 11 NA NA 0 E(13) = A(3), E(15) = D(2)
13 A 2 1 1 3 13 4 15 NA NA 2 E(14) = A(4), E(15) = D(2)
14 A 3 1 2 4 14 3 19 NA NA 5 E(15) = D(2), E(17) = A(5)

Notice that we now have 3 customers in the system, 1 in service and 2 waiting in the queue. Reviewing the pending event set, we see that customer 2 will finally complete service and depart at time 15.

Table 2.11: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)
11 A 1 1 0 2 11 4 11 NA NA 0 E(13) = A(3), E(15) = D(2)
13 A 2 1 1 3 13 4 15 NA NA 2 E(14) = A(4), E(15) = D(2)
14 A 3 1 2 4 14 3 19 NA NA 5 E(15) = D(2), E(17) = A(5)
15 D 2 1 1 2 11 4 11 15 4 0 E(17) = A(5), E(19) = D(3)

Because customer 2 completes service at time 15 and customer 3 is waiting in the line, we see that customer 3’s attributes for \(S_{i}\) and \(W_{i}\) were updated within the table. Since customer 3 has started service and we know their service time of 4 minutes, we know that they will depart at time 19. The pending events set has been updated accordingly and indicates that the next event will be the arrival of customer 5 at time 17.

Table 2.12: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)
11 A 1 1 0 2 11 4 11 NA NA 0 E(13) = A(3), E(15) = D(2)
13 A 2 1 1 3 13 4 15 NA NA 2 E(14) = A(4), E(15) = D(2)
14 A 3 1 2 4 14 3 19 NA NA 5 E(15) = D(2), E(17) = A(5)
15 D 2 1 1 2 11 4 11 15 4 0 E(17) = A(5), E(19) = D(3)
17 A 3 1 2 5 17 2 22 NA NA 5 E(19) = D(3), E(19) = A(6)

Now, we have a very interesting situation with the pending event set. We have two events that are scheduled to occur at the same time, the departure of customer 3 at time 19 and the arrival of customer 6 at time 19. In general, the order in which events are processed that occur at the same time may affect how future events are processed. That is, the order of event processing may change the behavior of the system over time and thus the order of processing is, in general, important. However, in this particular simple system, the order of processing will not change what happens in the future. We will simply have an update of the state variables at the same time. In this instance, we will process the departure event first and then the arrival event. If you are not convinced that this will not make a difference, then I encourage you to change the ordering and continue the processing. In more complex system simulations, a priority indicator is attached to the events so that the events can be processed in a consistent manner. Rather than continuing the step-by-step processing of the events through time 31, we will present the completed table. We leave it as an exercise for the reader to continue the processing of the customers. The completed bookkeeping table at time 31 is as follows.

Table 2.13: Bank by hand simulation bookkeeping table.
System Variables
Entity Attributes
t E(t) N(t) B(t) Q(t) ID(i) A(i) ST(i) S(i) D(i) T(i) W(i) Pending E(t)
0 NA 0 0 0 NA NA NA NA NA NA NA NA
3 A 1 1 0 1 3 4 3 NA NA 0 E(7) = D(1), E(11) = A(2)
7 D 0 0 0 1 3 4 3 7 4 0 E(11) = A(2)
11 A 1 1 0 2 11 4 11 NA NA 0 E(13) = A(3), E(15) = D(2)
13 A 2 1 1 3 13 4 15 NA NA 2 E(14) = A(4), E(15) = D(2)
14 A 3 1 2 4 14 3 19 NA NA 5 E(15) = D(2), E(17) = A(5)
15 D 2 1 1 2 11 4 11 15 4 0 E(17) = A(5), E(19) = D(3)
17 A 3 1 2 5 17 2 22 NA NA 5 E(19) = D(3), E(19) = A(6)
19 D 2 1 1 3 13 4 15 19 6 2 E(19) = A(6)
19 A 3 1 2 6 19 4 24 NA NA 5 E(21) = A(7), E(22) = D(4)
21 A 4 1 3 7 21 3 28 NA NA 7 E(22) = D(4), E(24) = D(5)
22 D 3 1 2 4 14 3 19 22 8 5 E(24) = D(5), E(27) = A(8)
24 D 2 1 1 5 17 2 22 24 7 5 E(27) = A(8), E(28) = D(6)
27 A 3 1 2 8 27 2 31 NA NA 4 E(28) = D(6), E(31) = D(7)
28 D 2 1 1 6 19 4 24 28 9 5 E(31) = D(7)
31 D 1 1 0 7 21 3 28 31 10 7 E(33) = D(8)

Given that we have simulated the bank over the time frame from 0 to 31 minutes, we can now compute performance statistics for the bank and measure how well it is operating. Two statistics that we can easily compute are the average time spent waiting in line and the average time spent in the system.

Let \(n\) be the number of customers observed to depart during the simulation. In this simulation, we had \(n = 7\) customers depart during the time frame of 0 to 31 minutes. The time frame over which we analyze the performance of the system is call the simulation time horizon. In this case, the simulation time horizon is fixed and known in advance. When we perform a computer simulation experiment, the time horizon is often referred to as the simulation run length (or simulation replication length). To estimate the average time spent waiting in line and the average time spent in the system, we can simply compute the sample averages ( \(\overline{T}\) and \(\bar{W})\) of the observed quantities (\(T_{i}\) and \(W_{i}\))for each departed customer.

\[\overline{T} = \frac{1}{7}\sum_{i = 1}^{7}T_{i} = \frac{4 + 4 + 6 + 8 + 7 + 9 + 10}{7} = \frac{48}{7} \cong 6.8571\]

\[\overline{W} = \frac{1}{7}\sum_{i = 1}^{7}W_{i} = \frac{0 + 0 + 2 + 5 + 5 + 5 + 7}{7} = \frac{24}{7} \cong 3.4286\]

Besides the average time spent in the system and the average time spent waiting, we might also want to compute the average number of customers in the system, the average number of customers in the queue, and the average number of busy tellers. You might be tempted to simply average the values in the \(N(t)\), \(B(t)\), and \(Q(t)\) columns of the bookkeeping table. Unfortunately, that will result in an incorrect estimate of these values because a simple average will not take into account how long each of the variables persisted with its values over time. In this situation, we are really interested in computed a time average. This is because the variables, \(N(t)\), \(B(t)\), and \(Q(t)\) are time-based variables. This type of variable is always encountered in discrete event simulations.

To make the discussion concrete, let’s focus on the number of customers in the queue, \(Q(t)\). Notice the number of customers in the queue, \(Q(t)\) takes on constant values during intervals of time corresponding to when the queue has a certain number of customers. \(Q(t) = \{ 0,1,2,3,\ldots\}\). The values of \(Q(t)\) form a step function in this particular case. The following figure illustrates the step function nature of this type of variable. A realization of the values of variable is called a sample path.

Number of customers waiting in queue for the bank simulation

Figure 2.17: Number of customers waiting in queue for the bank simulation

That is, for a given (realized) sample path, \(Q(t)\) is a function that returns the number of customers in the queue at time \(t\). The mean value theorem of calculus for integrals states that given a function, \(f( \bullet )\), continuous on an interval \((a, b)\), there exists a constant, c, such that

\[\int_{a}^{b}{f\left( x \right)\text{dx}} = f(c)(b - a)\]

\[f\left( c \right) = \frac{\int_{a}^{b}{f\left( x \right)\text{dx}}}{(b - a)}\]

The value, \(f(c)\), is called the mean value of the function. A similar function can be defined for \(Q(t)\) This function is called the time-average (and represents the mean value of the \(Q(t)\) function):

\[\overline{Q}\left( n \right) = \frac{\int_{t_{0}}^{t_{n}}{Q\left( t \right)\text{dt}}}{t_{n} - t_{0}}\]

This function represents the average with respect to time of the given state variable. This type of statistical variable is called time-based because \(Q(t)\) is a function of time.

In the particular case where \(Q(t)\) represents the number of customers in the queue, \(Q(t)\) will take on constant values during intervals of time corresponding to when the queue has a certain number of customers. Let \(Q\left( t \right) = \ q_{k}\ \)for\(\ t_{k - 1} \leq t \leq t_{k}\). Then, the time-average can be rewritten as follows:

\[\overline{Q}\left( t \right) = \frac{\int_{t_{0}}^{t_{n}}{Q\left( t \right)\text{dt}}}{t_{n} - t_{0}} = \sum_{k = 1}^{n}\frac{q_{k}(t_{k} - t_{k - 1})}{t_{n} - t_{0}}\] Note that \(q_{k}(t_{k} - t_{k - 1})\) is the area under the curve, \(Q\left( t \right)\) over the interval \(t_{k - 1} \leq t \leq t_{k}\) and because \[t_{n} - t_{0} = \left( t_{1} - t_{0} \right) + \left( t_{2} - t_{1} \right) + \left( t_{3} - t_{2} \right) + \ \cdots + \left( t_{n - 1} - t_{n - 2} \right) + \ \left( t_{n} - t_{n - 1} \right)\] we can write, \[t_{n} - t_{0} = \sum_{k = 1}^{n}{t_{k} - t_{k - 1}}\]

The quantity \(t_{n} - t_{0}\) represents the total time over which the variable is observed. Thus, the time average is simply the area under the curve divided by the amount of time over which the curve is observed. From this equation, it should be noted that each value of \(q_{k}\) is weighted by the length of time that the variable has the value. If we define, \(w_{k} = (t_{k} - t_{k - 1})\), then we can re-write the time average as:

\[\overline{Q}\left( t \right) = \frac{\int_{t_{0}}^{t_{n}}{Q\left( t \right)\text{dt}}}{t_{n} - t_{0}} = \sum_{k = 1}^{n}\frac{q_{k}(t_{k} - t_{k - 1})}{t_{n} - t_{0}} = \frac{\sum_{k = 1}^{n}{q_{k}w_{k}}}{\sum_{i = 1}^{n}w_{k}}\]

This form of the equation clearly shows that each value of \(q_{k}\) is weighted by:

\[\frac{w_{k}}{\sum_{i = 1}^{n}w_{k}} = \frac{w_{k}}{t_{n} - t_{0}} = \frac{(t_{k} - t_{k - 1})}{t_{n} - t_{0}}\]

This is why the time average is often called the time-weighted average. If \(w_{k} = 1\), then the time-weighted average is the same as the sample average.

Now we can compute the time average for \(Q\left( t \right),\ N\left( t \right)\) and \(B(t)\). Using the following formula and noting that \(t_{n} - t_{0} = 31\)

\[\overline{Q}\left( t \right) = \sum_{k = 1}^{n}\frac{q_{k}(t_{k} - t_{k - 1})}{t_{n} - t_{0}}\]

We have that the numerator computes as follows:

\[\sum_{k = 1}^{n}{q_{k}\left( t_{k} - t_{k - 1} \right)} = 0\left( 13 - 0 \right) + \ 1\left( 14 - 13 \right) + \ 2\left( 15 - 14 \right) + \ 1\left( 17 - 15 \right) + \ 2\left( 19 - 17 \right)\]

\[\ + 1\left( 19 - 19 \right) + \ 2\left( 21 - 19 \right) + \ 3\left( 22 - 21 \right) + \ 2\left( 24 - 22 \right) + \ \]

\[1\left( 27 - 24 \right) + \ 2\left( 28 - 27 \right) + \ 1\left( 31 - 28 \right) = 28\]

And, the final time-weighted average number in the queue ss:

\[\overline{Q}\left( t \right) = \frac{28}{31} \cong 0.903\]

The average number in the system and the average number of busy tellers can also be computed in a similar manner, resulting in:

\[\overline{N}\left( t \right) = \frac{52}{31} \cong 1.677\]

\[\overline{B}\left( t \right) = \frac{24}{31} \cong 0.7742\]

The value of \(\overline{B}\left( t \right)\) is most interesting for this situation. Because there is only 1 teller, the fraction of the tellers that are busy is 0.7742. This quantity represents the utilization of the teller. The utilization of a resource represents the proportion of time that the resource is busy. Let c represent the number of units of a resource that are available. Then the utilization of the resource is defined as:

\[\overline{U} = \frac{\overline{B}\left( t \right)}{c} = \frac{\int_{t_{0}}^{t_{n}}{B\left( t \right)\text{dt}}}{c(t_{n} - t_{0})}\]

Notice that the numerator of this equation is simply the total time that the resource is busy. So, we are computing the total time that the resource is busy divided by the total time that the resource could be busy, \(c(t_{n} - t_{0})\), which is considered the utilization.

In the following section, we will explore the basic elements of a process-oriented simulation.