plot_steady_states

% Read the matrices P0 and C

0001 %% Read the matrices P0 and C
0002 [P0,Txt,Raw]=xlsread('P_matrix'); % Original probability transition matrix
0003 [C,Txt,Raw]=xlsread('C'); % Optimal minimal energy perturbation matrix
0004 P = P0+C; 
0005 
0006 %% Compute pi_0 as the steady-state of P0
0007 [V,D] = eig(P0');
0008 % find the eigenvector corresponding to eigenvalue 1
0009 for k = 1:128
0010     if(abs(D(k,k)- 1) < 10^(-6))
0011         pi0 = V(:,k)/sum(V(:,k)); % pi0 is probability vector
0012         break;
0013     end
0014 end
0015 
0016 %% Define pi_d
0017 pid1 = 0.015525 * ones(64,1);
0018 pid2 = 10^(-4) * ones(64,1);
0019 pid = [pid1;pid2];
0020 
0021 %% Compute the steady-state of the controlled matrix P
0022 [V,D] = eig(P');
0023 % find the eigenvector corresponding to eigenvalue 1
0024 for k = 1:128
0025     if(abs(D(k,k)- 1) < 10^(-5))
0026         pic = V(:,k)/sum(V(:,k)); % pic is probability vector
0027         break;
0028     end
0029 end
0030 
0031 %% Plot
0032 states = [1:128];
0033 figure; plot(states, pi0, 'r', states, pid, 'x-', states, pic, 'g+');
0034 legend('Initial steady-state distribution', 'Desired steady-state distribution', 'Controlled steady-state distribution')
0035 
0036 
0037