M/G/c bank example
Mathematica code accompanying handout for the course Stochastic Simulation
Jan-Pieter Dorsman / April 15 2017

(* Function that defines parameters and performs preliminary work *)
preliminaries[] := Module[{},
						(* Model parameters *)						
						c =3; (* Number of bank clerks *)
						λ = 60; (* Arrival rate *)
						lowerParamServiceTime = 0; (* The lower bound on the uniform service time *)
						upperParamServiceTime = (1/15);(* The upper bound on the uniform service time *)
						numberOfHoursInADay =7; (* The number of hours to be simulated *)
						
						(* Definitions *)
						KINDARRIVAL =1;
						KINDDEPARTURE =2;

						(* Statistical Counters *)
						totalNumberOfArrivals = 0;
						totalWaitingTimeOfCustomers =0;
						totalAmountOfTimeClerkIsBusy =0;
						totalSurfaceUnderGraph = 0; (* counter used for the average queue length *)					
						
						(* Other parameters *)
						eventList = {}; (* Events in the list have the form {timeOfOccurence, kind}, e.g. {5.2857, KINDDEPARTURE} *)
						queue = {}; (* Customers in the queue have the form {time of arrival}, e.g. {3.4874} *)
						simulationClock = 0; 

						currentNumberOfClerksBusy = 0;								
							
						SeedRandom[123]; (* Reset the random number generator and set its seed equal 12345 *)
					  ];

(* Function that draws a number from the uniform distribution on [a,b], using a standard uniform sample given by the function RandomReal[] *)
drawUniformNumber[a_, b_]  := a + (b-a) * RandomReal[];

(* Function that draws a number from the exponential ( distribution, using a standard uniform sample given by the function RandomReal[] *)
drawExponentialNumber[rate_] := ((-Log[RandomReal[]])/rate);

(* Function that adds an event to the event list *)
addEventToEventList[event_] := Module[{},
									eventList  = Sort[Append[eventList, event]];
								];

(* Function that removes the next event from the event list and returns it *)
removeAndReturnNextEventFromList[] := Module[{nextEvent},
											nextEvent = First[eventList];
											eventList = Drop[eventList, 1];
											nextEvent
										];

(* Function that removes the next customer from the queue and returns it *)
removeAndReturnNextCustomerFromQueue[] := Module[{nextCustomer},
											nextCustomer = First[queue];
											queue = Drop[queue, 1];
											nextCustomer
										];

(* Function that starts the simulation: schedules the first arrival to start up the event list *)
initialise[] :=Module[{},
					firstArrivalTime = drawExponentialNumber[λ];				
					addEventToEventList[{firstArrivalTime, KINDARRIVAL}]; 
				];

(* Function that processes an arrival event *)
handleArrivalEvent[arrivalTime_] := Module[{nextArrivalTime, serviceTime},
											(* Handle the arrival *)
											totalNumberOfArrivals += 1;				

											If[currentNumberOfClerksBusy < c,
												(* The arriving customer can be taken into service immediately *)
												serviceTime = drawUniformNumber[lowerParamServiceTime, upperParamServiceTime];
												addEventToEventList[{simulationClock + serviceTime, KINDDEPARTURE}];
												currentNumberOfClerksBusy +=1												
											,
												 (* The arriving customer is forced to queue *) 
													queue = Append[queue, arrivalTime];
											];

											(* Now that the arrival is handled, schedule the next arrival *)
											nextArrivalTime = simulationClock + drawExponentialNumber[λ];
											addEventToEventList[{nextArrivalTime, KINDARRIVAL}];
									];

(* Function that processes a departure event *)
handleDepartureEvent[departureTime_] := Module[{nextServiceTime},
											(* Check whether there is a customer queued up.  *)
	
											If[Length[queue] > 0,
												(* There is a customer queued: put this customer into service and schedule departure *)
												arrivalTimeOfCustomerToBePutIntoService = removeAndReturnNextCustomerFromQueue[];
												totalWaitingTimeOfCustomers += (simulationClock
																				-arrivalTimeOfCustomerToBePutIntoService);
												
												nextserviceTime = drawUniformNumber[lowerParamServiceTime, upperParamServiceTime];
												addEventToEventList[{simulationClock + nextserviceTime,KINDDEPARTURE}];
											,
												(* There is no customer queued: make the clerk available again *)
												currentNumberOfClerksBusy -=1;
											];
										];

(* Function that takes the next event from the event list and processes it. *)
handleNextEvent[] :=Module[{nextEvent, eventTime, eventKind},
							nextEvent = removeAndReturnNextEventFromList[];
							eventTime = nextEvent[[1]];
							eventKind = nextEvent[[2]];

							(* Update counters and parameters *)
							totalSurfaceUnderGraph += (eventTime - simulationClock)*Length[queue];
							totalAmountOfTimeClerkIsBusy += (eventTime - simulationClock)*currentNumberOfClerksBusy;

							simulationClock = eventTime;

							(* Checks whether the next event is an arrival or a departure and calls functions accordingly *)
							If[eventKind==KINDARRIVAL , 
								handleArrivalEvent[eventTime]; 
							,
								handleDepartureEvent[eventTime];
							];
						];

(* Function that reports the input parameters and the findings of the simulation program *)
report[] := Module[{},
						Print["========================================"];
						Print["= Simulation run of M/G/c bank example ="];
						Print["========================================"];
						Print["Input parameters:"];
						Print["c equals ", c, "."];
						Print["λ equals ", λ, " per hour."];
						Print["Lower bound on uniform service time is ", lowerParamServiceTime, " hours."];
						Print["Upper bound on uniform service time = ", upperParamServiceTime, " hours."];
						Print["Number of hours in a day equals ", numberOfHoursInADay, "."];
						Print["========================================"];
						Print["Results:"];
						Print["Duration of the simulation is ", simulationClock, " hours"];
						
						Print["The total number of customers that arrived equals ", totalNumberOfArrivals ,"."];
						Print["The average queue length during the day equals ", totalSurfaceUnderGraph/simulationClock ,"."];
						Print["The average waiting time of customers equals ", totalWaitingTimeOfCustomers / totalNumberOfArrivals, " hour(s)."];
						Print["The average fraction of the time a bank clerk has been busy equals ", totalAmountOfTimeClerkIsBusy / simulationClock / c, "."];

			];

(* Main Function *)
runProgram[]  := Module[{},
							preliminaries[];
							initialise[];
				
							(* Handle events when the time of the earliest upcoming event, eventList[[1]][[1]], occurs before the end of the day, marked by numberOfHoursInADay *)				
							While[eventList[[1]][[1]] < numberOfHoursInADay,
									handleNextEvent[];
							];
			
							report[];				
			];

(* Actually run the program: *)

runProgram[]