%              Stationary.m
%
%  Steady-state probability distribution of a Markov chain
%  
%
% J. Lunze
%
% from: ES  3. Auflage: 11.1.2017
% 23.2.2022
% Version of 16.3.2022
%
%  
%
clear
close all
Gp=[0 0.8 1; 0.3 0.2 0; 0.7 0 0];
p0=[1; 0; 0];                   % Initial state distribution
Tmax=51;
P1=dinitial(Gp, zeros(3,1), eye(3,3), zeros(3,1), p0, Tmax);
figure(1)
subplot(2,1,1)
plot([0:Tmax-1], P1, 'b');
hold on
plot(0, P1(1,1), 'bo', 0, P1(1,2), 'bo', 0, P1(1,3), 'bo');
plot([0 0], [1.1, -0.1], 'k:');
latexxlabel('$$k$$');
xticks([0 10 20 30 40 50]);
latexylabel('$$p_i$$');
rotateY;
yticks([0 1]);
axis([-1 Tmax-1 -0.1 1.1]);
%
%  stationary probability distribution
pbar=inv(eye(3,3)-Gp + ones(3,3))*ones(3,1);
latextext(Tmax+1,0.55, '$$\bar{p}_1$$')
latextext(Tmax+1,0.15, '$$\bar{p}_2$$')
latextext(Tmax+1,0.35, '$$\bar{p}_3$$')
plot(Tmax-1, pbar(1), 'bo');
plot(Tmax-1, pbar(2), 'bo');
plot(Tmax-1, pbar(3), 'bo');
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Stationary distribution of a reducible Markov chain
%
figure(2)
%
alpha=0.45;
beta=0.4;
Gp=[0.5-alpha 0.5 0 0;
    0.5 0.5-beta 0 0;
    alpha 0 1 0;
    0 beta 0 1];
p01=[1; 0; 0; 0];
Tmax=10;
P1=dinitial(Gp, zeros(4,1), eye(4,4), zeros(4,1), p01, Tmax);
subplot(4,1,[1,2])
plot([0:Tmax-1], P1, 'b');
hold on
plot(0, P1(1,1), 'bo', 0, P1(1,2), 'bo', 0, P1(1,3), 'bo', 0, P1(1,4), 'bo');
plot([0 0], [1.1, -0.1], 'k:');
xleer;
xticksempty([0 1 2 3 4 5 6 7 8]);
latexylabel('$$p_i$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
latextext(8,0.8,'$$p_3$$');
%
p02=[0; 1; 0; 0];
P2=dinitial(Gp, zeros(4,1), eye(4,4), zeros(4,1), p02, Tmax);
%
subplot(4,1,[3,4])
plot([0:Tmax-1], P2, 'b');
hold on
plot(0, P2(1,1), 'bo', 0, P2(1,2), 'bo', 0, P2(1,3), 'bo', 0, P2(1,4), 'bo');
plot([0 0], [1.1, -0.1], 'k:');
latexxlabel('$$k$$');
xticks([0 1 2 3 4 5 6 7 8]);
latexylabel('$$p_i$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
latextext(8,0.8,'$$p_4$$');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      Birth-death process
%
figure(3)
%
alpha=0.45;
Tmax=16;
N1=11;
Gp=diag(alpha*ones(N1-1,1),1)+diag((1-alpha)*ones(N1-1,1),-1);
Gp(1,1)=alpha;
Gp(N1,N1)=1-alpha;
%
p0=zeros(N1,1);
p0(1,1)=1;
P=dinitial(Gp, zeros(N1,1), eye(N1,N1), zeros(N1,1), p0, Tmax);

%  stationary probability distribution
pbar=inv(eye(N1,N1)-Gp + ones(N1,N1))*ones(N1,1);
%
subplot(6,1,[1,2]);
stairs([0:Tmax-1], P(:,1), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,1), 'bo');
xleer;
xticksempty([0 5 10 15 20 25 30]);
latexylabel('$$p_0$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(1), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(1), '$$\bar{p}_0$$');
latextitle('$${\tt Stationary:}$$ $$\alpha=$$ 0.45');
latextext(0, 1.3, '$$\alpha=0.45$$');
hold off
%
subplot(6,1,[3,4]);
stairs([0:Tmax-1], P(:,3), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,3), 'bo');
xleer;
xticksempty([0 5 10 15 20 25 30]);
latexylabel('$$p_2$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(3), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(3), '$$\bar{p}_2$$');
hold off
%
subplot(6,1,[5,6]);
stairs([0:Tmax-1], P(:,11), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,11), 'bo');
latexxlabel('$$k$$');
xticks([0 5 10 15 20 25 30]);
latexylabel('$$p_{10}$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(11), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(11), '$$\bar{p}_{10}$$');
hold off

%%
figure(4)
alpha=0.5;
Gp=diag(alpha*ones(N1-1,1),1)+diag((1-alpha)*ones(N1-1,1),-1);
Gp(1,1)=alpha;
Gp(N1,N1)=1-alpha;
P=dinitial(Gp, zeros(N1,1), eye(N1,N1), zeros(N1,1), p0, Tmax);
stairs([0:Tmax-1], P, 'b', 'LineWidth',1);
%  stationary probability distribution
pbar=inv(eye(N1,N1)-Gp + ones(N1,N1))*ones(N1,1);
%
subplot(6,1,[1,2]);
stairs([0:Tmax-1], P(:,1), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,1), 'bo');
xleer;
xticksempty([0 5 10 15 20 25 30]);
latexylabel('$$p_0$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(1), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(1), '$$\bar{p}_0$$');
latextitle('$${\tt Stationary:}$$ $$\alpha=$$ 0.5');
latextext(0, 1.3, '$$\alpha=0.5$$');
hold off
%
subplot(6,1,[3,4]);
stairs([0:Tmax-1], P(:,3), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,3), 'bo');
xleer;
xticksempty([0 5 10 15 20 25 30]);
latexylabel('$$p_2$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(3), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(3), '$$\bar{p}_2$$');
hold off
%
subplot(6,1,[5,6]);
stairs([0:Tmax-1], P(:,11), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,11), 'bo');
latexxlabel('$$k$$');
xticks([0 5 10 15 20 25 30]);
latexylabel('$$p_{10}$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(11), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(11), '$$\bar{p}_{10}$$');
hold off



%%
figure(5)
alpha=0.55;
Gp=diag(alpha*ones(N1-1,1),1)+diag((1-alpha)*ones(N1-1,1),-1);
Gp(1,1)=alpha;
Gp(N1,N1)=1-alpha;
P=dinitial(Gp, zeros(N1,1), eye(N1,N1), zeros(N1,1), p0, Tmax);
stairs([0:Tmax-1], P, 'b', 'LineWidth',1);
%  stationary probability distribution
pbar=inv(eye(N1,N1)-Gp + ones(N1,N1))*ones(N1,1);
%
subplot(6,1,[1,2]);
stairs([0:Tmax-1], P(:,1), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,1), 'bo');
xleer;
xticksempty([0 5 10 15 20 25 30]);
latexylabel('$$p_0$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(1), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(1), '$$\bar{p}_0$$');
latextitle('$${\tt Stationary:}$$ $$\alpha=$$ 0.55');
latextext(0, 1.3, '$$\alpha=0.55$$');
hold off
%
subplot(6,1,[3,4]);
stairs([0:Tmax-1], P(:,3), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,3), 'bo');
xleer;
xticksempty([0 5 10 15 20 25 30]);
latexylabel('$$p_2$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(3), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(3), '$$\bar{p}_2$$');
hold off
%
subplot(6,1,[5,6]);
stairs([0:Tmax-1], P(:,11), 'b', 'LineWidth',1);
hold on
plot([0:Tmax-1], P(:,11), 'bo');
latexxlabel('$$k$$');
xticks([0 5 10 15 20 25 30]);
latexylabel('$$p_{10}$$');
rotateY;
yticks([0 1]);
axis([-0.5 Tmax-1 -0.1 1.1]);
plot(Tmax-1, pbar(11), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(Tmax, pbar(11), '$$\bar{p}_{10}$$');
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Periodic and aperiodic Markov chains
%
%figure(6)
%
G1=[0 0 1 0; 1 0 1 0; 0 1 0 1; 0 1 0 0];
G2=[0 0 1 0 0 0; 0 0 1 0 0 1; 0 1 0 0 1 0; 0 1 0 0 0 0; 0 0 0 1 0 0; 1 0 0 0 0 0];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Controlled Markov chain
%
figure(7)
pm=0.6;
pw=0.54;
pm1=0.3;
pm2=pm-pm1;
pw1=0.27;
pw2=pw-pw1;
%   state-transition probabilities of states with 6 successor states
%   edges are enumerated clockwise
pi1=pw1*(1-pm);
pi2=pw2*(1-pm);
pi3=pw2*pm1;
pi4=(1-pm)*pm1;
pi5=(1-pw)*pm2;
pi6=pw1*pm2;
%
%
Gp=zeros(25,25);
% Gp(1,1)=1-pw;
Gp(1,2)=(1-pw)*pm1;
Gp(1,3)=(1-pw)*pm2;
Gp(2,1)=pw1;
% Gp(2,2)=1-Gp(1,2)-Gp(4,2)-Gp(5,2)-Gp(3,2);
Gp(2,3)=pw1*pm2;
Gp(2,4)=(1-pw)*pm1;
Gp(2,5)=(1-pw)*pm2;
Gp(3,1)=pw2;
Gp(3,2)=pw2*pm1;
% Gp(3,3)=1-Gp(1,3)-Gp(2,3)-Gp(5,3)-Gp(6,3);
Gp(3,5)=(1-pw)*pm1;
Gp(3,6)=(1-pw)*pm2;
Gp(4,2)=pw1*(1-pm);
% Gp(4,4)=1-Gp(2,4)-Gp(5,4)-Gp(7,4)-Gp(8,4);
Gp(4,5)=pw1*pm2;
Gp(4,7)=(1-pw)*pm1;
Gp(4,8)=(1-pw)*pm2;
Gp(5,2)=pw2*(1-pm);
Gp(5,3)=pw1*(1-pm);
Gp(5,4)=pw2*pm1;
% Gp(5,5)=1-Gp(2,5)-Gp(3,5)-Gp(4,5)-Gp(6,5)-Gp(8,5)-Gp(9,5);
Gp(5,6)=pw1*pm2;
Gp(5,8)=(1-pw)*pm1;
Gp(5,9)=(1-pw)*pm2;
Gp(6,3)=pw2*(1-pm);
Gp(6,5)=pw2*pm1;
% Gp(6,6)=1-Gp(3,6)-Gp(5,6)-Gp(9,6)-Gp(10,6);
Gp(6,9)=(1-pw)*pm1;
Gp(6,10)=(1-pw)*pm2;
Gp(7,4)=pw1*(1-pm);
% Gp(7,7)=1-Gp(4,7)-Gp(8,7)-Gp(11,7)-Gp(12,7);
Gp(7,8)=pw1*pm2;
Gp(7,11)=(1-pw)*pm1;
Gp(7,12)=(1-pw)*pm2;
Gp(8,4)=pw2*(1-pm);
Gp(8,5)=pw1*(1-pm);
Gp(8,7)=pw2*pm1;
% Gp(8,8)=1-Gp(4,8)-Gp(5,8)-Gp(7,8)-Gp(9,8)-Gp(12,8)-Gp(13,8);
Gp(8,9)=pw1*pm2;
Gp(8,12)=(1-pw)*pm1;
Gp(8,13)=(1-pw)*pm2;
Gp(9,5)=pw2*(1-pm);
Gp(9,6)=pw1*(1-pm);
Gp(9,8)=pw2*pm1;
% Gp(9,9)=1-Gp(5,9)-Gp(6,9)-Gp(8,9)-Gp(10,9)-Gp(13,9)-Gp(14,9);
Gp(9,10)=pw1*pm2;
Gp(9,13)=(1-pw)*pm1;
Gp(9,14)=(1-pw)*pm2;
Gp(10,6)=pw2*(1-pm);
Gp(10,9)=pw2*pm1;
% Gp(10,10)=1-Gp(6,10)-Gp(9,10)-Gp(14,10)-Gp(15,10);
Gp(10,14)=(1-pw)*pm1;
Gp(10,15)=(1-pw)*pm2;
Gp(11,7)=pw1*(1-pm);
% Gp(11,11)=1-Gp(7,11)-Gp(12,11)-Gp(16,11);
Gp(11,12)=pw1*pm2;
Gp(11,16)=(1-pw)*pm2;
Gp(12,7)=pw2*(1-pm);
Gp(12,8)=pw1*(1-pm);
Gp(12,11)=pw2*pm1;
% Gp(12,12)=1-Gp(7,12)-Gp(8,12)-Gp(11,12)-Gp(13,12)-Gp(16,12)-Gp(17,12);
Gp(12,13)=pw1*pm2;
Gp(12,16)=(1-pw)*pm1;
Gp(12,17)=(1-pw)*pm2;
Gp(13,8)=pw2*(1-pm);
Gp(13,9)=pw1*(1-pm);
Gp(13,12)=pw2*pm1;
% Gp(13,13)=1-Gp(8,13)-Gp(9,13)-Gp(12,13)-Gp(14,13)-Gp(17,13)-Gp(18,13);
Gp(13,14)=pw1*pm2;
Gp(13,17)=(1-pw)*pm1;
Gp(13,18)=(1-pw)*pm2;
Gp(14,9)=pw2*(1-pm);
Gp(14,10)=pw1*(1-pm);
Gp(14,13)=pw2*pm1;
% Gp(14,14)=1-Gp(9,14)-Gp(10,14)-Gp(13,14)-Gp(15,14)-Gp(18,14)-Gp(19,14);
Gp(14,15)=pw1*pm2;
Gp(14,18)=(1-pw)*pm1;
Gp(14,19)=(1-pw)*pm2;
Gp(15,10)=pw2*(1-pm);
Gp(15,14)=pw2*pm1;
% Gp(15,15)=1-Gp(10,15)-Gp(14,15)-Gp(19,15);
Gp(15,19)=(1-pw)*pm1;
Gp(16,11)=pw2*(1-pm);
Gp(16,12)=pw1*(1-pm);
% Gp(16,16)=1-Gp(11,16)-Gp(12,16)-Gp(17,16)-Gp(20,16);
Gp(16,17)=pw1*pm2;
Gp(16,20)=(1-pw)*pm2;
Gp(17,12)=pw2*(1-pm);
Gp(17,13)=pw1*(1-pm);
Gp(17,16)=pw2*pm1;
% Gp(17,17)=1-Gp(12,17)-Gp(13,17)-Gp(16,17)-Gp(18,17)-Gp(20,17)-Gp(21,17);
Gp(17,18)=pw1*pm2;
Gp(17,20)=(1-pw)*pm1;
Gp(17,21)=(1-pw)*pm2;
Gp(18,13)=pw2*(1-pm);
Gp(18,14)=pw1*(1-pm);
Gp(18,17)=pw2*pm1;
% Gp(18,18)=1-Gp(13,18)-Gp(14,18)-Gp(17,18)-Gp(19,18)-Gp(21,18)-Gp(22,18);
Gp(18,19)=pw1*pm2;
Gp(18,21)=(1-pw)*pm1;
Gp(18,22)=(1-pw)*pm2;
Gp(19,14)=pw2*(1-pm);
Gp(19,15)=pw1*(1-pm);
Gp(19,18)=pw2*pm1;
% Gp(19,19)=1-Gp(14,19)-Gp(15,19)-Gp(18,19)-Gp(22,19);
Gp(19,22)=(1-pw)*pm1;
Gp(20,16)=pw2*(1-pm);
Gp(20,17)=pw1*(1-pm);
% Gp(20,20)=1-Gp(16,20)-Gp(17,20)-Gp(21,20)-Gp(23,20);
Gp(20,21)=pw1*pm2;
Gp(20,23)=(1-pw)*pm2;
Gp(21,17)=pw2*(1-pm);
Gp(21,18)=pw1*(1-pm);
Gp(21,20)=pw2*pm1;
% Gp(21,21)=1-Gp(17,21)-Gp(18,21)-Gp(20,21)-Gp(22,21)-Gp(23,21)-Gp(24,21);
Gp(21,22)=pw1*pm2;
Gp(21,23)=(1-pw)*pm1;
Gp(21,24)=(1-pw)*pm2;
Gp(22,18)=pw2*(1-pm);
Gp(22,19)=pw1*(1-pm);
Gp(22,21)=pw2*pm1;
% Gp(22,22)=1-Gp(18,22)-Gp(19,22)-Gp(21,22)-Gp(24,22);
Gp(22,24)=(1-pw)*pm1;
Gp(23,20)=pw2*(1-pm);
Gp(23,21)=pw1*(1-pm);
% Gp(23,23)=1-Gp(20,23)-Gp(21,23)-Gp(24,23)-Gp(25,23);
Gp(23,24)=pw1*pm2;
Gp(23,25)=(1-pw)*pm2;
Gp(24,21)=pw2*(1-pm);
Gp(24,22)=pw1*(1-pm);
Gp(24,23)=pw2*pm1;
% Gp(24,24)=1-Gp(21,24)-Gp(22,24)-Gp(23,24)-Gp(25,24);
Gp(24,25)=(1-pw)*pm1;
Gp(25,23)=pw2*(1-pm);
Gp(25,24)=pw1*(1-pm);
% Gp(25,25)=1-Gp(23,25)-Gp(24,25);
%
Gdiag=sum(Gp)';
Gp=Gp+diag(ones(25,1)-Gdiag);
Gp0=Gp;            %  Graph without control
%  stationary probability distribution
pbar=inv(eye(25,25)-Gp + ones(25,25))*ones(25,1);
%
%  Evaluation of the result
%  according to the number of workpieces in the system
P=[];
P(1)=pbar(1);
P(2)=pbar(2)+pbar(3);
P(3)=pbar(4)+pbar(5)+pbar(6);
P(4)=pbar(7)+pbar(8)+pbar(9)+pbar(10);
P(5)=pbar(11)+pbar(12)+pbar(13)+pbar(14)+pbar(15);
P(6)=pbar(16)+pbar(17)+pbar(18)+pbar(19);
P(7)=pbar(20)+pbar(21)+pbar(22);
P(8)=pbar(23)+pbar(24);
P(9)=pbar(25);
subplot(2,1,1);
plot([0:8],P,'b--');
hold on
plot([0:8],P,'bo');
latextitle('$${\tt Stationary:}$$ Controlled machine tool');
latexylabel('$$p_{\rm s}$$');
rotateY;
latexxlabel('$$z_{\rm s}$$');
xticks([0 1 2 3 4 5 6 7 8]);
yticks([0 0.1 0.2]);
axis([0 8 0 0.2]);
%
%  Evaluation of the result
%  according to the difference of workpieces in both waiting areas
D=[];
D(9)=pbar(11);                                   % difference=4
D(8)=pbar(7)+pbar(16);
D(7)=pbar(4)+pbar(12)+pbar(20);
D(6)=pbar(2)+pbar(8)+pbar(17)+pbar(24);
D(5)=pbar(1)+pbar(5)+pbar(13)+pbar(21)+pbar(25); % difference=0
D(4)=pbar(3)+pbar(9)+pbar(18)+pbar(24);
D(3)=pbar(6)+pbar(14)+pbar(22);
D(2)=pbar(10)+pbar(19);
D(1)=pbar(15);                                   % difference=-4
subplot(2,1,2);
plot([-4:4],D,'b--');
hold on
plot([-4:4],D,'bo');
latexylabel('$$p_{\rm \Delta}$$');
rotateY;
latexxlabel('$$z_1-z_2$$');
xticks([-4 -3 -2 -1 0 1 2 3 4]);
yticks([0 0.1 0.2]);
axis([-4 4 0 0.25]);

%%   Control: change the direction of workpieces if the waiting areas are (nearly) full
%
pm=0.6;
pw=0.54;
pm1=0.3;
pm2=pm-pm1;
pw1=0.27;
pw2=pw-pw1;
alpha=0.7;
Gp=Gp0;
%   Reduce the workpiece flow to service station 1 if there are already 
%   3 workpieces in the waiting area
%   states: 7, 12, 17
pw1=0.27;
pw1=pw1*alpha;
pw2=pw-pw1;
Gp(8,7)=pw2*pm1;
Gp(11,7)=pw1*(1-pm);
Gp(12,7)=pw2*(1-pm);
Gp(11,12)=pw1*pm2;
Gp(13,12)=pw2*pm1;
Gp(16,12)=pw1*(1-pm);
Gp(17,12)=pw2*(1-pm);
Gp(16,17)=pw1*pm2;
Gp(18,17)=pw2*pm1;
Gp(20,17)=pw1*(1-pm);
Gp(21,17)=pw2*(1-pm);
%   Reduce the workpiece flow to service station 2 if there are already 
%   3 workpieces in the waiting area
%   states: 10, 14, 18
pw1=0.27;
pw2=pw-pw1;
pw2=pw2*alpha;
pw1=pw-pw2;
Gp(9,10)=pw1*pm2;
Gp(14,10)=pw1*(1-pm);
Gp(15,10)=pw2*(1-pm);
Gp(13,14)=pw1*pm2;
Gp(15,14)=pw2*pm1;
Gp(18,14)=pw1*(1-pm);
Gp(19,14)=pw2*(1-pm);
Gp(17,18)=pw1*pm2;
Gp(19,18)=pw2*pm1;
Gp(21,18)=pw1*(1-pm);
Gp(22,18)=pw2*(1-pm);
%   Stop the workpiece flow to service station 1
%   states: 11, 16, 20, 23
pw1=0;
pw2=pw;
Gp(12,11)=pw2*pm1;
Gp(16,11)=pw2*(1-pm);
Gp(17,16)=pw2*pm1;
Gp(20,16)=pw2*(1-pm);
Gp(21,20)=pw2*pm1;
Gp(23,20)=pw2*(1-pm);
Gp(24,23)=pw2*pm1;
Gp(25,23)=pw2*(1-pm);
%   Stop the workpiece flow to service station 2
%   states: 15, 19. 22. 24
pw2=0;
pw1=pw;
Gp(14,15)=pw1*pm2;
Gp(19,15)=pw1*(1-pm);
Gp(18,19)=pw1*pm2;
Gp(22,19)=pw1*(1-pm);
Gp(21,22)=pw1*pm2;
Gp(24,22)=pw1*(1-pm);
Gp(23,24)=pw1*pm2;
Gp(25,24)=pw1*(1-pm);
%
Gdiag=sum(Gp)';
Gp=Gp+diag(ones(25,1)-Gdiag);
%  stationary probability distribution
pbar=inv(eye(25,25)-Gp + ones(25,25))*ones(25,1);
%
%  Evaluation of the result
%  according to the number of workpieces in the system
P=[];
P(1)=pbar(1);
P(2)=pbar(2)+pbar(3);
P(3)=pbar(4)+pbar(5)+pbar(6);
P(4)=pbar(7)+pbar(8)+pbar(9)+pbar(10);
P(5)=pbar(11)+pbar(12)+pbar(13)+pbar(14)+pbar(15);
P(6)=pbar(16)+pbar(17)+pbar(18)+pbar(19);
P(7)=pbar(20)+pbar(21)+pbar(22);
P(8)=pbar(23)+pbar(24);
P(9)=pbar(25);
subplot(2,1,1);
plot([0:8],P,'r')
plot([0:8],P,'ro')
hold off
%
%  Evaluation of the result
%  according to the difference of workpieces in both waiting areas
D=[];
D(1)=pbar(11);                                   % difference=4
D(2)=pbar(7)+pbar(16);
D(3)=pbar(4)+pbar(12)+pbar(20);
D(4)=pbar(2)+pbar(8)+pbar(17)+pbar(24);
D(5)=pbar(1)+pbar(5)+pbar(13)+pbar(21)+pbar(25); % difference=0
D(6)=pbar(3)+pbar(9)+pbar(18)+pbar(24);
D(7)=pbar(6)+pbar(14)+pbar(22);
D(8)=pbar(10)+pbar(19);
D(9)=pbar(15);                                   % difference=-4
subplot(2,1,2);
plot([4:-1:-4],D,'r');
plot([4:-1:-4],D,'ro');
hold off





%%
epsfigc('Stationary');
