Simulation of classical queueing systems with Mathematica

This paper describes a simple and useful tool for simulating every kind of Markovian arrival process and Markovian service process queues. With this program, the user can study and observe the behavior of Marakovian systems and modify the parameters for various explorations.

by  C?. D'Apice, B?. D'Auria, G?. Gargiulo and S?. Salerno

A queue or queuing system is a stochastic system whose structure is as follows: an input, or incoming flow (the arriving stream of customers of some sort), one or more servers, a waiting line and some kind of output or exit. In other words, a queue is a dynamical system where some kind of commodity flows in a non-deterministic way. The incoming flow is made up by the users requiring service, which may be in a finite or (potentially) infinite number. The system's output is made up by the users who received service. The servers or service stations (clerks, computers…) may in turn be finite or infinite, and in typical examples serve one customer at a time. The waiting line is formed by users finding all servers to be busy. Again, we may have a finite upper bound for its length (finite storage capacity), or no bound (infinite storage capacity). When the maximum length is finite, customers are refused until a server is idle. A queuing system may be specified in various ways. Thus, the arriving stream of the population of potential users may be described through the probability distribution of the interarrival times (time elapsed between two consecutive arrivals).

More precisely, the interarrival time (Subscript[τ, i])= time interval between the i-th and the (i+1)-th service request within the observation period.

The time sequence i is the realisation of a discrete parameter stochastic process, called the input process of the queuing system. A typical assumption on this process (at least when studying steady states) is that it can be modelled as a stationary renewal process, i.e. the interarrival times are identically distributed independent random variables (i.d.i.r.v.). In particular, the exponential distribution is among the most used models. By the same token, the service process is given by the service times i = service time for the i-th customer (Subscript[C, i]). Again, it is generally assumed that i are i.d.i.r.v. and the exponential distribution is among the most used models. The average service frequency is by definition 1/average-servicetime. The behaviour of a queuing system may then be analysed once we are given the waiting line's length and the queuing discipline, at least when the system's input and output processes are independent. Among the most widely used are: first-come-first-serve (FCFS), last-come-first-serve (LCFS) and random order of service (RO) of which FCFS is the most common. In order to describe the stationary (long run) behaviour of the queuing system we shall use the following parameters:

k = number of customers in the system;
δ = time spent by a generic user in the system;
l = number of the customers in the waiting line;
η = time spent by a generic user in the waiting line or waiting time;
h = number of customers being served or number of busy servers;

A shorthand notation for queuing systems requires a string of five symbols, where the first two refer to the input and output (arrival-departure) processes, and denote respectively the probability distribution of interarrival and service times as follows:

({{Symbol, Distribution}, {M, exponential}, {(Subscript[ϵ, r]), order r Erlang}, {(Subscript[H, r]), order r hyperexponential}, {G, general}, {D, deterministic}})

The third, fourth and fifth symbols denote respectively the number m of servers, the lenght L of the waiting line and the size n of the population of potential users. The five symbols are separated by /, as in M / M / 1 / 150 / 200.

Package description

The package is loaded with the standard Mathematica command.

<<  Simulation`SimulationQueue.m

Then, the kind of queue to be simulated has to be specified
The Queue command accepts a four-parameter input:

\(TextForm\`The\ first\ parameter\ represents\ the\ interarrival\ times\ distribution . 
The\ package\ \ manages\ the\ Markov - 
        type\ distributions\ with\ a\ parameter\ depending\ on\ the\ number\ of\ customers\ in\ the\ system\ at\ a\ given\ epoch . 

            Three\ predefined\ functions\ are\ \(available : 
\(Poisson[λ]\  : \ \ λ\_k\)\) = λ\)
Bernoulli[λ]  :  (Subscript[λ, k]) = λ*(n - 
            k)
\(TextForm\`Pascal[α, β]\  : \ λ\_k = α + β\ k\)
\(TextForm\`where\ k\ is\ the\ \ number\ of\ customers\ in\ the\ system\ and\ n\ is\ the\ number\ of\ customers\ in\ the\ whole\ population\ \((\ 
        in\ the\ Bernoullian\ case, \ 
        n\ must\ be\  < ∞)\)\)

The second parameter represents the service times distribution. The package  manages the Markov-type distributions with a constant parameter μ
The third parameter represents the number of servers (service stations) in the system. The fourth parameter represents the length of the waiting line.
The fifth parameter represents the number of customers in the whole population.
In general:

q := Queue[M[f[#]&], M[μ], m, K, M]

where f is a function whose output is λ>0 once k ( number of customers in the system) is given . Once f has been defined, both static and dynamical simulation are possible. For static simulation, the Simulation function may be used:

ss = Simulation[q, N]

where N is the number of users that will go through the system. Simulation 's output is a 21st order list as reported in Table 1. For dynamic simulation, the DynamicSimulation function may be used:

ds = DynamicSimulation[q, N]

where N is the number of users that will go through the system. DynamicSimulation 's output is a 21+26 order list. The first 21 outputs are equal to the 21 Simulation 's outputs. The 22nd output ( ds[[22]] )is the number of events generated by the system's simulation. The outputs from ds[[23]][n] to ds[[46]][n] are functions that report the interest value vs. the number of the event ( 1≤n≤ds[22]=events). The list of DynamicSimulation 's outputs is reported in Table 2.

Notebook description

The "Queue" notebook makes usage of the "Simulation`QueueSimulation"  package which is more simple and immediate. Once the "Code.nb" file has been loaded the notebook is displayed as follows.

Simulation

Load the package         Functions      Help

Clicking on "Load the package"  automatically loads the  "Simulation`QueueSimulation" package. "Help" opens an help notebook. "Functions" shows the notebook's functions, viz

Simulate     Examples      Theory

"Examples" opens a notebook with evaluated simulations together with the corresponding graphs. "Theory" opens a notebook with basic news about "Queuing Theory." "Simulate" allows the choice the kind of simulation:

SSimulation		DSimulation		DSimulation3D

"SSimulation" shows the homonymous function:

SSimulation[arrtimfun, sertimfun, numser, queuel, popol, numit]

SSimulation cyclically calls  the Simulation function generating a set of simulations of indexed queues with index i variable from 0.25 to 14.5 with  0.25 as step. arrtimfun is a function of the cycle index i and returns the (Subscript[λ, k]) function that characterizes the i-th queue's interrarrival times' process.
Example Poisson[7]& to have every queue with Poisson process of the interarrival times (Subscript[λ, k]) = λ
      Pascal[#,Sqrt[#]]& to have the i-th queue with Pascal process of the interrarrival times (Subscript[λ, k]) = i + k (Sqrt[i])
sertimfun is a function of cycle index i which returns μ  that characterizes the i-th queue's service times' Poisson process.
Example 8& to have every queue with μ=8
       Abs[Ln[#]]& to have the i-th queue with μ=|Ln(i)|
numser is the maximum number of servers that are present in the queue system.  In fact the simulation with index i will be repeated numser times in order to have graphs vs. i and parameterized by m  (1≤m≤numser is the number of server that are currently present in the queue system)
queuel is the length of the waiting queue; it is the same for all the simulated queues. popol is the number of the users in the population; ( same for all  queues). numit is the number of users that will go through the system, ( same for all  queues).

After making a selection, this panel will be shown:

SSimulation		DSimulation		DSimulation3D		Graph[■]		Help		Theory

"Graph" calls  the homonymous function Graph[n] allows drawing a graph of any chosen output n (1≤n≤21) vs. the cycle index i parameterized by m (number of servers in the system 1≤m≤numser)

"DSimulation" calls  the homonymous function:

DSimulation[arrtim, sertim, numser, queuel, popol, numit]

DSimulation calls the inner function DynamicSimulation 
arrtim is the (Subscript[λ, k]) function.
Example: Poisson[7] to have a queue system with a Poissonian interarrival times process (Subscript[λ, k]) = λ. Pascal[10,Sqrt[5]] to have a queue system with a Pascal interarrival times process (Subscript[λ, k]) = 10 + k (Sqrt[5.])  sertimfun is the parameter μ that characterizes the Poisson service times process. 
Example: 8 to have μ=8 .Abs[Ln[6]] to have  μ=|Ln(6)| numser is the number of servers that are present in the queue system.  queuel is the length of the waiting queue;  (the same for all  the simulated queues). popol is the number of the users in the population;  ( same for all  queues). numit is the number of users that will go through the system; ( same for all  queues).

After simulating this choice will be shown:

SSimulation		DSimulation		DSimulation3D		GraphD[■]		Help		Theory

"GraphD" calls  the homonymous function GraphD[n]  allows drawing a graph of any chosen output n  (1≤n≤24) vs. time.

"DSimulation3D" calls  the homonymous function:

DSimulation3D[arrtimfun, sertimfun, numser, queuel, popol, numit]

DSimulation3D calls cyclically the inner DynamicSimulation function generating a set of simulations of indexed queues with index i variable from 0.25 to 14.5 with  0.25 as step. arrtimfun is a function of cycle index i and returns the (Subscript[λ, k]) function that characterizes the i-th queue's interrarrival times' process.
Example Poisson[7]& to have every queue with Poisson process of the interarrival times (Subscript[λ, k]) = λ. Pascal[#,Sqrt[#]]& to have the i-th queue with Pascal process of the interrarrival times (Subscript[λ, k]) = i + k (Sqrt[i])
sertimfun is a function of cycle index i and returns μ  that characterizes the i-th queue's service times' Poisson process.
Example 8& to have every queue with μ=8. Abs[Ln[#]]& to have the i-th queue with μ=|Ln(i)| numser is the number of servers that are present in the queue system;  (the same for all  simulated queues).  queuel is the length of the waiting queue;  ( same for all  queues). popol is the number of the users in the population;  ( same for all  queues). numit is the number of users that will go through the system;  ( same for all  queues).

After simulating this choice will be shown:

Simulation		DSimulation		DSimulation3D		GraphD3D[■]		Help		Theory

"GraphD3D" calls  the homonymous function GraphD3D[n] allows drawing a graph of any chosen output n  (1≤n≤21) vs. the cycle index i and the number of event k (1≤k≤events).

Output explanation

Here there are some examples of output generated by the "Code" notebook.

SSimulation[Sqrt, (1/Sqrt[#])&, 60, 10, 4000]

In this case we plot two output values for a statistical simulation.
The first represents the average serviced users rate (i.e. the output flow of users that get service) as a function of the cycle index i, which has been chosen in such a way as to equal Ao (Average Amplitude of Arrival Traffic). We remark that as the incoming traffic increases, the Average Serviced Users Rate decreases, i.e. the Average Service Time increases, as it does when the number of servers is increased.

The second represents the average lost users rate (i.e.the flow of users refused by the system)  as a function of the cycle index i, which has again been chosen in such a way as to equal Ao (Average Amplitude of Arrival Traffic). We remark that as the incoming traffic increases, the refusal rate increases.

The results agree with the theoretical ones.

Graph[5];

Graph[8];

DSimulation[3, 1/3, 10, 3, 200]

In this case we plot three output values for a 2D dynamical simulation.

GraphD[1];

GraphD[11];

GraphD[14];

GraphD[16];

DSimulation3D[Sqrt, (1/Sqrt[#])&, 60, 5, 100]

In this case we plot three output values for a 3D dynamical simulation.

The first represents the number of users in the waiting queue in the time epochs (events) when a variation of the system's state occurs (arrival, departure) as a function of the cycle index i, which has been as usual chosen in such a way as to equal Ao (Average Amplitude of Arrival Traffic). We observe that, as the incoming traffic increases, the queue is more and more often saturated - in the present example, 5 users in the waiting line.

The second represents the average amplitude of serviced traffic (i.e. the flow of customers that leave the servers after being served), again when a variation occurs and as a function of i as above. The values given by outputs ds[[1]] - ds[[21]] converge as time goes by - i.e. as the events occur - to the corresponding outputs of the statical simulation (namely, ss[[1]] - ss[[21]). This means that the statical symulation yields the section of this graph  along the plane x = events (ds[[22]]). We remark that, as the incoming traffic increases, the amplitude of serviced traffic decreases i.e the service time increases.

The third represents the blocked service probability (i.e. the probability that all servers are busy), as usual, when a variation occurs and as a function of the usual cycle index i. We remark that, as the incoming traffic increases, both the  blocked service probability and the convergence rate increase.

The results agree with the theoretical ones.

GraphD3D[3];

GraphD3D[13];

GraphD3D[19];

A simple practical example

Let us suppose a bank has to buy a bancomat system; two choices are available, either two slow machines or a  single fast one. The choice is to be done according to the minimum waiting time (the Average Arrival Users Rate being expected to be the same in both cases). Let us assume that the Average Arrival Users Rate has a frequency of 1/30 s and let us evaluate the two alternatives according to the resulting average.

Let us first simulate with two slow servers (Average Serviced User Rate 1/120s) and  Average Arrival Users Rate 1/30s

SSimulation[Poisson[1/30]&, 1/120&, 2, ∞, ∞, 1000]

gf2=Graph[20, 2, PlotRange -> All, DisplayFunction -> Identity];

Show[gf2, DisplayFunction -> $DisplayFunction];

MeanMM2 = Mean[Table[Trial[2][[k, 20]], {k, 1, 57, 1}]]

MeanMM2 represents the average waiting time over all simulations.
Let us now simulate with one fast server (Average Serviced User Rate 2/120 s)  and the same Average Arrival Users Rate 1/30 s

SSimulation[Poisson[1/30]&, 2/120&, 1, ∞, ∞, 1000]

gf1=Graph[20, PlotRange -> All, DisplayFunction -> Identity];

Show[gf1, DisplayFunction -> $DisplayFunction];

MeanMM1 = Mean[Table[Trial[1][[k, 20]], {k, 1, 57, 1}]]

MeanMM1 again represents the average waiting time over all simulations.

Means = 
      Plot[{MeanMM1, MeanMM2}, {k, 0.50, 14.5}, 
        PlotStyle -> {{RGBColor[0, 0, 1], Dashing[{0.15, 0.05}]}, {
              RGBColor[0, 1, 0], Dashing[{0.05, 0.02}]}}, 
        DisplayFunction -> Identity];

Show[gf1, gf2, Means, DisplayFunction -> $DisplayFunction];

MeanMM1 < MeanMM2

Since MeanMM2<MeanMM1 the two-server system will be chosen.

Conclusion

As you can see with this package you can study the static behavior of M/M systems, confirming the theorical results and additionally you can utilize it to obtain information on its dinamical behavior. For the latter no general closed formulas are known.

Appendix

Table 1(Simulation' s Outputs)

ss[[1]] = no :   Number of users arrived to the systemss[[2]] = Lambdao :   Average arrival users ratess[[3]] = 
      Ao :  Average Amplitude of Arrival Traffic 
          (λ/μ)ss[[4]] = ns :  Number of users serviced by the systemss[[5]] = Lamdas :  Average serviced users ratess[[6]] = As : Average Amplitude of Serviced Traffic ss[[7]] = np :   Number of users lost by the systemss[[8]] = Lambdap :   Average lost users ratess[[9]] = Ap :  Average Amplitude of Lost Traffic ss[[10]] = 
      Sp :  Blocked System Probability (Waiting Queue is Full)ss[[11]] = 
      Pip :  Lost Probability 
          (Arrival Probability conditioned on Blocked System)ss[[12]] = 
      Sr :  Blocked Service Probability (All Servers are Busy)ss[[13]] = 
      Pir :  Waiting Probability 
          (Arrival Probability conditioned on Blocked Service)ss[[14]] = 
      time :  Time spent by the System during the Simulation 
          (in the same unit of λ and μ)ss[[15]] = 
      Omegatime :  
        Time spent by the System in state of Blocked System during 
          the Simulationss[[16]] = 
      blocktime :  
        Time spent by the System in state of Blocked Service during
           the Simulationss[[17]] = avk :  Average number of users in the Systemss[[18]] = 
      avh :  Average number of contemporary busy serversss[[19]] = 
      avl :  Average number of users in the waiting queuess[[20]] = aveta :  Average waiting timess[[21]] = 
      avdelta :  
        Average system time 
          (Average waiting time  +  Average service time)

\(\(Table\ 2
DynamicSimulation' s\ Outputs)\)\)

\(\(ds[
      \([22]\) = 
        \(events : 
            number\ of\ events\ processed\ by\ the\ simulation

              \(ds[\([23]\)]\)[n] = 
          \(ftime[n] : 
              Time\ of\ event\ occurence
01\  - \ \(ds[\([24]\)]\)[n] = 
            \(fk[n] : \ 
                Number\ of\ users\ in\ the\ System
02\  - \ 
                  \(ds[\([25]\)]\)[n] = 
              \(fl[n] : \ 
                  Number\ of\ contemporary\ busy\ servers
03\  - \ 
                    \(ds[\([26]\)]\)[n] = 
                \(fh[n] : \ 
                    Number\ of\ users\ in\ the\ waiting\ queue
04\  - \ 
                      \(ds[\([27]\)]\)[n] = 
                  \(favk[n] : \ 
                      Average\ number\ of\ users\ in\ the\ System
05\  - \ 
                        \(ds[\([28]\)]\)[n] = 
                    \(favh[n] : \ 
                        Average\ number\ of\ contemporary\ busy\ servers
06
                            \  - \ \(ds[\([29]\)]\)[n] = 
                      \(favl[n] : \ 
                          Average\ number\ of\ users\ in\ the\ waiting\ queue
                              
07\  - \ \(ds[\([30]\)]\)[n] = 
                        \(fno[n] : \ 
                            Number\ of\ users\ arrived\ to\ the\ system
08
                                \  - \ \(ds[\([31]\)]\)[n] = 
                          \(fnp[n] : \ 
                              Number\ of\ users\ lost\ by\ the\ system
09
                                  \  - \ \(ds[\([32]\)]\)[n] = 
                            \(fns[n] : \ 
                                Number\ of\ users\ serviced\ by\ the\ system

                                    10\  - \ \(ds[\([33]\)]\)[n] = 
                              \(fnr[n] : \ 
                                  Number\ of\ users\ that\ waited\ in\ the\ 
                                      system
11\  - \ \(ds[\([34]\)]\)[n] = 
                                \(fAo[n] : \ 
                                    Average\ Amplitude\ of\ Arrival\ Traffic\ 
                                        \((λ/μ)\)
12\  - \ \ 
                                      \(ds[\([35]\)]\)[n] = 
                                  \(fAp[n] : \ 
                                      Average\ Amplitude\ of\ Lost\ Traffic\ 
                                        
13\  - \ \ \(ds[\([36]\)]\)[n] = 
                                    \(fAs[n] : \ 
                                        Average\ Amplitude\ of\ Serviced\ 
                                        Traffic\ 
14\  - \ 
                                        \(ds[\([37]\)]\)[n] = 
                                      \(fLambdao[n] : \ 
                                        Average\ arrival\ users\ rate
15\  - 
                                        \ \(ds[\([38]\)]\)[n] = 
                                        \(fLambdap[n] : \ 
                                        Average\ lost\ users\ rate
16\  - \ 
                                        \(ds[\([39]\)]\)[n] = 
                                        \(fLambdas[n] : \ 
                                        Average\ serviced\ users\ rate
17
                                        \  - \ \(ds[\([40]\)]\)[n] = 
                                        \(fSp[n] : \ 
                                        Blocked\ System\ Probability\ 
                                        \((Waiting\ Queue\ is\ Full)\)
18
                                        \  - \ \(ds[\([41]\)]\)[n] = 
                                        \(fPip[n] : \ 
                                        Lost\ Probability\ 
                                        \((Arrival\ Probability\ conditioned\ 
                                        on\ Blocked\ System)\)
19\  - \ 
                                        \(ds[\([42]\)]\)[n] = 
                                        \(fSr[n] : \ 
                                        Blocked\ Service\ Probability\ 
                                        \((All\ Servers\ are\ Busy)\)
20\  - 
                                        \ \(ds[\([43]\)]\)[n] = 
                                        \(fPir[n] : \ Waiting\ Probability\ 
                                        \((Arrival\ Probability\ conditioned\ 
                                        on\ Blocked\ Service)\)
21\  - \ 
                                        \(ds[\([44]\)]\)[n] = 
                                        \(fOmegatime[n] : \ \ 
                                        Time\ spent\ by\ the\ System\ in\ 
                                        state\ of\ Blocked\ System
22\  - \ 
                                        \(ds[\([45]\)]\)[n] = 
                                        \(fblocktime[n] : \ 
                                        Time\ spent\ by\ the\ System\ in\ 
                                        state\ of\ Blocked\ Service
23\  - \ 
                                        \(ds[\([46]\)]\)[n] = 
                                        \(faveta[n] : \ Average\ waiting\ time
                                        
24\  - \ \(ds[\([47]\)]\)[n] = 
                                        favdelta[n] : \ Average\ system\ time
                                        \)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)

REFERENCES

Kleinrock Queueing Systems: Computer Applications 1976 Volume II Wiley-Interscience Publication
Donald Gross Carl M.Harris: Fundamentals of  Queuing Theory  Third Edition Wiley-Interscience Publication
G. Albano, A. Cavallone, C. D'Apice and S. Salerno: Mathematica and didactical innovation   IMS 99 Vienna
G. Albano, C. D'Apice, G. Gargiulo: Visualisation of four dimensional curves ICTMT4

ABOUT THE AUTHORS
(*Need first names, not just initials*)
C. D'Apice is a researcher in  Mathematical Analysis in the Department of Information Engineering and Applied Mathematics at the University of Salerno. His current scientific interests are calculus of variations, queueing systems, visualization in Mathematics.
B. D'Auria is a student in Electronic Engineering at University of Salerno.
G. Gargiulo is a researcher in  Mathematical Analysis in the Department of Information Engineering and Applied Mathematics at the University of Salerno. His current scientific interests are calculus of variations, queueing systems, visualization in Mathematics.
S. Salerno is a Full Professor in Mathematical Analysis in the Department of Information Engineering and Applied Mathematics at the University of Salerno. His current scientific interests are calculus of variations, number theory, queueing systems, visualization in Mathematics.
Address
Department of Information Engineering and Applied Mathematics at the University of Salerno
Via Ponte don Melillo
84084 Fisciano (SA)
Italy
e-mail   dapice@diima.unisa.it
             gargiulo@diima.unisa.it
             salerno@ponza.dia.unisa.it

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file classicCuing.nb containing Mathematica code for the material described in this article.