plot_prop5

% Read the matrices P0 and C

0001 
0002 %% Read the matrices P0 and C
0003 [P0,Txt,Raw]=xlsread('P_matrix'); % Original probability transition matrix
0004 [C,Txt,Raw]=xlsread('C'); % Optimal minimal energy perturbation matrix
0005 P = P0+C; 
0006 OneVec = ones(128,1);
0007 
0008 %% Define pi_0
0009 [V,D] = eig(P0');
0010 % find the eigenvector corresponding to eigenvalue 1
0011 for k = 1:128
0012     if(abs(D(k,k)- 1) < 10^(-6))
0013         pi0 = V(:,k)/sum(V(:,k)); % pi0 is probability vector
0014         break;
0015     end
0016 end
0017 
0018 %% Define pi_d
0019 pid1 = 0.015525 * ones(64,1);
0020 pid2 = 10^(-4) * ones(64,1);
0021 pid = [pid1;pid2];
0022 
0023 %% Find \pi_d(s): the steady-state of Q(s)
0024 NormVectorQs = [];
0025 UpperBound = [];
0026 i = 1;
0027 for s = 0:10^(-3):1
0028     Qs = (1-s)*P0 + s * OneVec*pid'; 
0029     [Vs,Ds] = eig(Qs');
0030     % find the eigenvector corresponding to eigenvalue 1
0031     for k = 1:128
0032         if(abs(Ds(k,k)- 1) < 10^(-6))
0033             pids = Vs(:,k)/sum(Vs(:,k)); % pids is probability vector
0034             break;
0035         end
0036     end   
0037     NormVectorQs(i) = norm(pids-pid); 
0038     UpperBound(i) = (2*(1-s)* norm(pi0 - pid))/(2-s);
0039     i = i +1;
0040 end
0041 figure; plot([0:10^(-3):1],NormVectorQs); hold on; plot([0:10^(-3):1],UpperBound, '--r');
0042 xlabel('s'); legend('||pid(s) - pid||', 'Upper bound');
0043 
0044 
0045 
0046 
0047 
0048 
0049 
0050 
0051 
0052