% NOTE: This script can easily be adapted to do apply the FES method (for 
% multiple nodes) 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 0.0 1.0 1.0;
  1.0 0.0 0.0 0.0;
  0.0 0.5 0.0 0.0;
  0.0 0.5 0.0 0.0  ];


% The visit ratios are 4,3,2,1
mu=4:-1:1;


% number of jobs in the network
K=3;

% number of nodes in the network
N=4;

% STEP 1 of ext. FES

% sn1 defines the nodes of subnetwork 1, sn2 respectively
sn1=[ 1 2 ];
sn2=[ 3 4 ];

% 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);
e=[ 1.0 (inv(p)*b)' ];

fprintf( '\n=== visit ratios ===\n' );
fprintf( ' one column per node\n' );
e

% STEP 2 of ext. FES

% construction of the short-circuited network: 
% subnetwork 1 is short-circuited
fesMu=mu( sn2 );
fesE=e( sn2 );
fesN=N-length( sn2 );

% the analysis of the short-circuited network (subnetwork 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)
G=[ 1 ];                    % normalization constant

% 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, : ) ];
  % normalization constant
  G=[ G; G(end)/fesLambda(end)];

end

% store the results
fesLambdaSN2=fesLambda;
GSN2=G;

fprintf( '==============================================\n' );
fprintf( '===== MVA on the short-circuited network.===== \n' );
fprintf( '       subnetwork 1 is short-circuited         \n' );
fprintf( '==============================================\n' );
fprintf( '\n=== mean response times ===\n' );
fprintf( ' c1: node 3 / c2: node 4\n' )
fesMeanT(end,:)
fprintf( '\n=== mean number of jobs ===\n' );
fprintf( ' c1: node 3 / c2: node 4\n' )
fesMeanK(end,:)
fprintf( '\n=== throughputs ===\n' );
fprintf( ' load-dependent overall throughput\n' );
fesLambdaSN2
fprintf( '\n=== normalization constant ===\n' );
fprintf( ' one line per number of jobs in the network (0:K)\n' )
GSN2

% construction of the short-circuited network: 
% subnetwork 2 is short-circuited
fesMu=mu( sn1 );
fesE=e( sn1 );
fesN=N-length( sn1 );

% the analysis of the short-circuited network (subnetwork 2 short-circuited) is
% carried out with the mean value analysis
fesMeanT=[ ];              % mean response time
fesLambda=[ ];             % the throughput of the network
fesMeanK=zeros( 1, fesN ); % mean number of jobs per node (initialized to 0)
G=[1];                     % normalization constant

% 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, : ) ];

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

end

% store the values
fesLambdaSN1=fesLambda;
GSN1=G;

fprintf( '==============================================\n' );
fprintf( '===== MVA on the short-circuited network.===== \n' );
fprintf( '       subnetwork 2 is short-circuited         \n' );
fprintf( '==============================================\n' );
fprintf( '\n=== mean response times ===\n' );
fprintf( ' c1: node 1 / c2: node 2\n' )
fesMeanT(end,:)
fprintf( '\n=== mean number of jobs ===\n' );
fprintf( ' c1: node 1 / c2: node 2\n' )
fesMeanK(end,:)
fprintf( '\n=== throughputs ===\n' );
fprintf( ' load-dependent overall throughput\n' );
fesLambdaSN1
fprintf( '\n=== normalization constant ===\n' );
fprintf( ' one line per number of jobs in the network (0:K)\n' )
GSN1

% ALTERNATIVE: Analysis with the convolution method.

% Calculation of the function Fi by Eq. 7.61
%% F1= ( fesE(1) ./ fesMu(1) ) .* ones( K+1, 1 );
%% F2= ( fesE(2) ./ fesMu(2) ) .* ones( K+1, 1 );
%% F=[ F1 .^ ((0:K)'), F2 .^ ((0:K)') ]
%% 
%% % Initialization of the vector G (Eq. 8.5)
%% G=F(:,1);
%% 
%% % Convolution
%% G=conv(G,F(:,2));
%% G=G(1:K+1);
%% 
%% % Calculation of the throughput
%% lambda= G(1:end-1) ./ G(2:end);

% STEP 3 of ext. FES

% save the old values
eOld=e;
muOld=mu;


fprintf( '==========================\n' );
fprintf( '===== Reduced network ====\n' );
fprintf( '==========================\n' );

fprintf( '=== Load-dependent service rates ===\n' );
fprintf( ' l1: node 1 / l2: node 2 / one columns per number of jobs in the node\n' );
[ fesLambdaSN1'; fesLambdaSN2' ]


% STEP 4 of ext. FES

G=conv( GSN1, GSN2 );

fprintf( '=== Normalization constant of the reduced network ===\n' );
fprintf( ' one entry per number of jobs in the network (0:K)\n' );
G=G( 1:K+1 )

% STEP 4 of ext. FES

% Mean number of jobs at each node: 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 (Eq. 8.14)
lambda= e .* ( G(K) / G(K+1) );

% mean response time: Little's Law
meanT= meanK ./ lambda;

fprintf( '=========================================\n' );
fprintf( '===== Results for the whole network ===== \n' );
fprintf( '=========================================\n' );
fprintf( '\n=== mean response times ===\n' );
fprintf( ' one column per node\n' )
meanT
fprintf( '\n=== mean number of jobs ===\n' );
fprintf( ' one column per node\n' )
meanK
fprintf( '\n=== throughputs ===\n' );
fprintf( ' one column per node\n' )
lambda

