% Multi-server queue equations with infinite population

% A.A. Trani (Mar 99)
%
% Revisions:

% 1) Changed function call to FCT.m to avoid use of FCT. Instead use of the Matlab
% factorial function (FACTORIAL) - October 23, 2003 (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 conditions

  S=2;
  Lambda=60/25;
  Mu = 60/45;
  nlast = 15;				
  
  Rho=Lambda/(S*Mu);   	

% 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  =   ',num2str(Po)])
  disp(' ')  
    disp([blanks(5),'Expected No. of customers in queue (Lq)  =   ',num2str(Lq)])
  disp(' ')  
    disp([blanks(5),'Expected No. of customers in system (L)   =   ',num2str(L)])
  disp(' ')  
    disp([blanks(5),'Average Waiting Time in Queue  =   ',num2str(Wq)])
  disp(' ') 
      disp([blanks(5),'Average Waiting Time in System (includes service)  =   ',num2str(W)])
  disp(' ') 

  % Fraction of time servers are idle
  
  FTI=0;
  for i=0:S-1
    Pi=(Rho*S)^i/factorial(i)*Po; 
    FTI=FTI+1/S*(S-i)*Pi; 
  end   
  
  stem(n,Pn,'*')
  xlabel('Number of entities')
  ylabel('Probability')
  grid
  
  % Save results to a file
  % Create a vector r with results [ Lambda Mu_S Po Rho Lq L W Wq] and append values in row format
  
  r(1,1) = Lambda; r(2,1) = Mu; r(3,1) = Po; r(4,1) = Rho; r(5,1) = Lq; r(6,1) = L; r(7,1) = W; r(8,1) = Wq; r(9,1) = S;
  
  fid = fopen ('queue_results.txt','a');
  fprintf(fid,'%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n',r);
  fclose(fid);

