% Multi-server queue equations with infinite population

% A.A. Trani (Mar 99)
% Revised February 2004 by A. Trani

% S = Number of servers
% Lambda = arrival rate
% Mu = Service rate per server
% Rho = utilization factor
% Po = Idle probability
% L = Expected no of entities in the system
% Lq = Expected no of entities in the queue
% nlast - last probability to be computed

% Initial condition


% Clear screen
clc

disp(' ')
disp('Multi-server Stochastic Queueing Analysis')
disp('Developed by A.A. Trani')
disp('Introduction to Transportation Engineering')
disp(' ')
S = input ('Number of Servers  ');

disp(' ')
disp('Enter the Arrival Rate - Lambda ');
Lambda = input ('(units are vehicles per time)   ');

disp(' ')
disp('Enter the Service Rate - Mu');
Mu = input ('(units are vehicles per time per server)   ');

disp(' ')
disp('Enter the last Probability to be Estimated ');
nlast = input (' ');

% Start computations

% Estimate utilization factor (Rho)
  
Rho=Lambda/(S*Mu); 
  
if Rho >= 1.0
        disp(' ')
        disp([blanks(5),'Queue Never Reaches Steady State Conditions '])
        disp([blanks(5),'The Results of the Analysis are Meaningless '])
        disp(' ')
else

% Find Po
    Po_inverse=0;
    sum_den=0;

  for i=0:S-1           % for the first term in the denominator (den_1) 
    den_1=(Lambda/Mu)^i/factorial(i);
    sum_den=sum_den+den_1;
  end

    den_2=(Lambda/Mu)^S/(factorial(S)*(1-Rho)); % for the second part of denominator (den_2)   		
    Po_inverse=sum_den+den_2;
    Po=1/Po_inverse;
	
% Find probabilities (Pn)

	Pn(1) = Po;		% Initializes the first element of Pn column vector to be Po
	n(1) = 0;		% Vector to keep track of number of entities in system
	
% loop to compute probabilities of n entities in the system

	for j=1:1:nlast		
	   n(j+1) = j;
	   if (j) <= S
	       Pn(j+1) = (Lambda/Mu)^j/factorial(j) * Po;
		else
		   Pn(j+1) = (Lambda/Mu)^j/(factorial(S) * S^(j-S)) * Po;
		end
	end

% Queue measures of effectiveness

  Lq=(Lambda/Mu)^S*Rho*Po/(factorial(S)*(1-Rho)^2);
              
  L=Lq+Lambda/Mu;
       
  Wq=Lq/Lambda; 
  
  W = L/Lambda;  

  disp(' ')
  disp([blanks(5),'Queueing Parameters '])
    disp(' ')
  disp([blanks(5),'System utilization (%) =   ',num2str(Rho * 100)])
  disp(' ')
  disp([blanks(5),'Idle probability (dim) =   ',num2str(Po)])
  disp(' ')  
    disp([blanks(5),'Expected No. of vehicles in queue (Lq)  =   ',num2str(Lq)])
  disp(' ')  
    disp([blanks(5),'Expected No. of vehicles in system (L)   =   ',num2str(L)])
  disp(' ')  
    disp([blanks(5),'Average Waiting Time in Queue (Time units) =   ',num2str(Wq)])
  disp(' ') 
      disp([blanks(5),'Average Waiting Time in System (includes service) (Time units)  =   ',num2str(W)])
  disp(' ') 
  
  
  stem(n,Pn,'*')
  xlabel('Number of entities')
  ylabel('Probability')
  grid
  
end
  