
mu=[ 1/1 1/30 1/50 1/50 ];
P=[ 0.000 0.500 0.250 0.250      % the routing probability matrix
    0.950 0.050 0.000 0.000  
    1.000 0.000 0.000 0.000  
    1.000 0.000 0.000 0.000 ];

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

K=8;
N=4;
epsilon = 0.02;

% === INITIALIZATION ===
meanK = ones(1,N) * K/N;
meanKminus1 = [];
lambda=[];
lambdas=[];
meanT=[];
err=[];

while 1
 
  % === estimate values (Eq.(9.1)) ===
  meanKminus1 = [ meanKminus1; (K-1)/K * meanK(end,:) ];

  % === one MVA iteration ===

  % mean number of jobs
  meanT = [ meanT; 1./mu .* ( meanKminus1(end,:) + 1 ) ];
  meanT(end,2) = 1/mu(2);

  % throughput
  lambda = [ lambda; K/(meanT(end,:)*e') ];
  lambdas = [ lambdas; e * lambda(end) ];

  % mean number of jobs
  meanK = [ meanK; lambdas(end,:) .* meanT(end,:) ];
   
  err = [ err; max( abs( meanK(end,:) - meanK(end-1,:) ) ) ];
  if err(end) <= epsilon
    break
  end
end

fprintf( '=== Estimates for meanK(K-1) ===\n' );
fprintf( ' rows = iterations / columns = nodes (1 - 4)\n ' );
meanKminus1
fprintf( '=== Mean response times ===\n' );
fprintf( ' rows = iterations / columns = nodes (1 - 4)\n ' );
meanT
fprintf( '=== Throughputs ===\n' );
fprintf( ' rows = iterations / columns =  nodes (1 - 4) and last column is the overall throughput\n ' );
[ lambdas, lambda ]
fprintf( '=== Mean number of jobs ===\n' );
fprintf( ' rows = iterations / columns = nodes (1 - 4)\n ' );
meanK
fprintf( '=== error ===\n' );
fprintf( ' rows = iterations\n' );
err
