% NOTE: This script can easily be adapted to do apply the FES method (for a
% single node) on other closed networks with M/M/1 nodes

% routing probabilities: P(i,j) contains the routing probability from node i to
% node j
P=[
  0.0 1.0 1.0
  0.6 0.0 0.0
  0.4 0.0 0.0; ];


% service rates 
mu=[ 1 0.65 0.4 ];

% number of jobs in the network
K=3;

% number of nodes in the network
N=3;

% sc defines the number of the node, which will be short-circuited
sc=1; % nodes 2 and 3 will be combined to a FES node

% computation of the visit ratios (note that node 1 is chosen to be
% the node with a visit ratio of 1)
p=P-eye(size(P));
b=-p(2:end,1);
p=p(2:end,2:end);

fprintf( '=== visit ratios ===\n one column per node\n' );
e=[ 1.0 (inv(p)*b)' ]


% construction of the short-circuited network
fesMu=mu( find( (1:N) ~= sc ) ); % service of the nodes != shortcut node
fesE=e( find( (1:N) ~= sc ) );   % visit ratios othe nodes != shortcut node
fesN=N-1;                        % number of nodes in the sc-network

% the analysis of the short-circuited network (node 1 short-circuited) is
% carried out with the mean value analysis

% STEP 1: initialization
fesMeanT=[ ];              % mean response time
fesLambda=[ ];             % the throughput of the network
fesMeanK=zeros( 1, fesN ); % mean number of jobs per node (initialized to 0)

% STEP 2:
for k = 1:K                

  % STEP 2.1: calculation of the mean response time of each node
  % when there are k jobs in the system (Eq. 8.38)
  % note that the values are appended to the previous intermediate
  % results
  fesMeanT=[
   fesMeanT;
   ( 1 ./ fesMu ) .* ( 1 + fesMeanK( end, : ) ) ];

  % STEP 2.2: calculation of the throughput (Eq. 8.39)
  fesLambda=[
   fesLambda;
   k ./ ( sum( fesE .* fesMeanT( end, : ) ) ) ];

  % STEP 2.3: calculation of the mean number of jobs (Eq. 8.41)
  fesMeanK=[
   fesMeanK;
   fesE .* fesLambda(end) .* fesMeanT( end, : ) ];
end

fprintf( '==============================================\n' );
fprintf( '===== MVA on the short-circuited network.===== \n' );
fprintf( '==============================================\n' );
fprintf( '\n=== mean response times ===\n' );
fprintf( ' c1: node 2 / c2: node 3\n' )
fesMeanT(end,:)
fprintf( '\n=== mean number of jobs ===\n' );
fprintf( ' c1: node 2 / c2: node 3\n' )
fesMeanK(end,:)
fprintf( '\n=== throughputs ===\n' );
fprintf( ' overall throughput\n' );
fesLambda(end,:)

% construction of the reduced network

% save the old visit ratios and service rates
eOld=e;
muOld=mu;


% calculation of the load-dependent service rates: line 1 contains
% the (constant) service rates of the sc-node, line 2 the load-
% dependent service rates (determined the load-dependent throughputs
% of the short-circuited network) of the FES node

fprintf( '\n=== Load-dependent service rates of the reduced network (l1=short-cut node; l2=FES node) === \n' );
mu=[ ones(1,K) .* mu( sc ); e(sc) .* fesLambda' ]

% visit ratios
e=[ e(sc), e(sc) ];

% Apply the load-dependent MVA

% STEP 1: Initialization
pi1=zeros( 4, 4 ); % pi1(k,j) contains the conditional probability that there are k-1 jobs in node 1 when there
                   % are j-1 jobs in the system
pi2=zeros( 4, 4 ); % likewise for node 2
pi1(1,1)=1;        % probability that there is no job in node 1/2 when there is no job in the system is 1
pi2(1,1)=1;
meanT=[];          % mean response time
lambda=[];         % throughput
G=[ 1 ];           % normalization constant

% STEP 2: MVA iteration
for k = 1:K

  % STEP 2.1: mean response time by Eq. 8.57
  meanT=[
    meanT;
     [ sum( ( (1:k) ./ mu(1,1:k) ) .* pi1(1:k,k)'), sum( ( (1:k) ./ mu(2,1:k) ) .* pi2(1:k,k)') ] ];

  % STEP 2.2: throughput by Eq. 8.58
  lambda=[
    lambda;
    k ./ sum( e .* meanT( end, : ) ) ];

  % normalization constant
  G=[ G; G(end) / lambda(end) ];

  % STEP 2.3: conditional probabilities by Eq. 8.59
  sum1 = 0;
  sum2 = 0;
  for l = k:-1:1
    pi1( l+1, k+1 ) = ( lambda(end) / mu(1,l) ) * pi1( l, k ) * e(1);
    sum1 = sum1 + pi1( l+1, k+1 );
    pi2( l+1, k+1 ) = ( lambda(end) / mu(2,l) ) * pi2( l, k ) * e(2);
    sum2 = sum2 + pi2( l+1, k+1 );
  end
  pi1( 1, k+1 ) = 1 - sum1;
  pi2( 1, k+1 ) = 1 - sum2;

end

% the mean number of jobs can be determined by the mean response times
% and the throughput
meanK=meanT(end,:) .* lambda(end) .* e;

fprintf( '=====================================================\n' );
fprintf( '===== Load-dependent MVA on the reduced network =====\n' );
fprintf( '=====================================================\n' );
fprintf( '\n=== mean response time ===\n' );
fprintf( ' c1: node 1 / c2: FES node\n' );
meanT(end,:)
fprintf( '\n=== throughput ===\n' );
lambda(end,:)
fprintf( '\n=== mean number of jobs ===\n' );
fprintf( ' c1: node 1 / c2: FES node\n' );
meanK(end,:)
fprintf( '\n=== normalization constant ===\n' );
fprintf( ' dependent on the number of jobs 0:K in the system\n' );
G

% recover the visit ratios and service rates of the original network
e=eOld;
mu=muOld;

% the mean number of jobs in the original system by Eq. 8.16
meanK=zeros( 1, N );
for k = 1:K
  meanK = meanK + (e./mu).^k .* ( G( K-k+1 ) / G( K+1 ) );
end

% throughputs
lambdas= e .* ( G(K) / G(K+1) );

% mean response time by Little's Law
meanT= meanK ./ lambdas;

fprintf( '=============================================\n' );
fprintf( '===== FES results for the whole network =====\n' );
fprintf( '=============================================\n' );
fprintf( '\n=== mean number of jobs ===\n' );
fprintf( ' c1: node 1 / c2: node 2 / c3: node 3\n' );
meanK
fprintf( '\n=== throughputs ===\n' );
fprintf( ' c1: node 1 / c2: node 2 / c3: node 3\n' );
lambda'
fprintf( '\n=== mean response time ===\n' );
fprintf( ' c1: node 1 / c2: node 2 / c3: node 3\n' );
meanT
