function [d,u,nq] = queue2(n,k,lambda,mu)
% Simulate delays in queue for queuing model (M/M/k)
% ------------------------------------
% n the number of customers
% k the number of servers
% lambda the arrival rate of customers
% mu the service rate
% ------------------------------------
% d the delay times for each customer
% u the utilization of servers
% nq the average queue length
% ------------------------------------
t = 0; % Simulation clock
ta = exprnd(lambda); % The time of next customer arrival
td = inf; % Departure time of customers (set to a dummy value)
served = []; % Customers being served
queue = []; % Queued customers
n_in = 0; % Number of customer arrivals
n_out = 0; % Number of customer departures
d = []; % Service delays
u = 0; % Utilization of server
% (time interval * fraction of servers having
% customer in that time interval)
nq = 0; % Number in queue statistic
% (time interval * length of queue)
% The main simulation loop
while n_out < n
% Next event is arrival
if ta < td
% Update statistics
% Number of customers entered the system
n_in = n_in + 1;
% Utilization statistic
u = u + (ta-t) * length(served) / k;
% Number of customers in queue statistic
nq = nq + (ta-t) * length(queue);
% Add either to server or queue
if length(served) < k
% Generate departure time, add to served vector and sort
% so that earliest departure is the first element
served = sort([served ta+exprnd(mu)]);
% Customer went straight to server => no delay
d(end+1)= 0;
else
% Add the arrival time to the queue vector
queue = [queue ta];
end
% Update simulation clock
t = ta;
% Generate next arrival time
ta = ta + exprnd(lambda);
% Set next departure time as the earliest departure time, i.e.,
% the first element of the served vector
td = served(1);
else
% Next event is departure
% Update statistics
n_out = n_out + 1;
u = u + (td-t) * length(served) / k;
nq = nq + (td-t) * length(queue);
% Departure, remove first element of served vector
served = served(2:end);
% Take next customer from queue, if queue is not empty
if ~isempty(queue)
% Add new customer to served list and generate departure time
served = sort([served td+exprnd(mu)]);
% Record delay: Now - time the customer entered the queue
d(end+1)= td - queue(1);
% Remove the customer, i.e., the first element from the queue
queue = queue(2:end);
end
% Update simulation clock
t = td;
% If servers are empty
if isempty(served)
% Set next departure time high so that the event cannot happen
td = inf;
else
% If there is customers in the servers, set the next departure
% time as the earliest departure time in served
td = served(1);
end
end
end
% Calculate utilization: In u, we have the added time intervals
% between the events times the fraction of servers having customer in the
% intervals. To get utilization we simply divide u with the simulation time
% in the end of the simulation
u = u / t;
% Calculate average queue length: Similar to utilization
nq = nq / t;