
mu = [ 1/3 1/1 1/1 1/1 1/1 1/20 ];           % the visit ratios
P=[ 0.000 0.125 0.125 0.125 0.125 0.500;     % the routing probability matrix
    1.000 0.000 0.000 0.000 0.000 0.000; 
    1.000 0.000 0.000 0.000 0.000 0.000; 
    1.000 0.000 0.000 0.000 0.000 0.000; 
    1.000 0.000 0.000 0.000 0.000 0.000; 
    1.000 0.000 0.000 0.000 0.000 0.000 ];

b=[ 1 0 0 0 0 0 ];      % linear equations system to calcualte the visit ratios
A=P - eye( size( P ) );
A(:,1)=[ 1 0 0 0 0 0 ];
e=b*inv(A);             % calculation of the visit ratios

K=5; % five iterations are needed

% INITIALIZATION
lambda=[];
lambdas=[];
meanT=[];
meanK=[ 0 0 0 0 0 0 ];

% ITERATION
for k = (1:K)

   % Step 2.1: Mean response times
   meanT = [ meanT; 1./mu .* ( meanK(end,:) + 1 ) ];
   % Node 6 is a IS node
   meanT(end,6) = 1/mu(6);

   % Step 2.2: Throughput
   lambda = [ lambda; k/(e * meanT(end,:)') ];
   lambdas = [ lambdas; lambda(end,:) * e ];

   % Step 2.3: Mean number of jobs
   meanK = [ meanK; lambdas(end,:) .* meanT(end,:) ];

end

% Calculation of mean value of T_S
meanKS=K-meanK(end,6);
lambdaS = lambdas(end,6);
meanTS = meanKS / lambdaS;

% OUTPUT
fprintf( '=== Mean response times ===\n' );
fprintf( ' rows = number of jobs in the network (0 - 5) / columns = nodes (1 - 6)\n ' );
meanT
fprintf( '=== Mean number of jobs ===\n' );
fprintf( ' rows = number of jobs in the network (0 - 5) / columns = nodes (1 - 6)\n ' );
meanK
fprintf( '=== Throughputs ===\n' );
fprintf( ' rows = number of jobs in the network (0 - 5) / columns =  nodes (1 - 6) and last column is the overall throughput\n ' );
[ lambdas, lambda ]

fprintf( '=== Other performance measures ===\n' );
meanKS
lambdaS
meanTS
