%
%      Queue.M
%
%  Queueing systems as Markov chains
%
% Jan Lunze
% from Wartesyst.m, ES, 3. Auflage: 26.1.2017
% 2.3.3022
% Version of 3.3.3022
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 
clear
close all
echo off
%   Arrival probability (pW) 
%   Service probability (pM)
%
pW=0.55;
pM=0.6;
% 
figure(1)
N=4;
G11=1-pW;
G12=(1-pW)*pM;
G32=pW*(1-pM);
G22=1-pW*(1-pM)-pM*(1-pW);
GNN=1-(1-pW)*pM;
Gp=[G11 G12 0 0;
    pW G22 G12 0;
    0  G32 G22 G12;
    0   0  G32 GNN];
p0=[1; 0; 0; 0];
Tend=7;
Time=[0:Tend-1];
X=dinitial(Gp, zeros(N,1), eye(N,N), zeros(N,1), p0, Tend);
%  stationary probability distribution
pbar=inv(eye(4,4)-Gp + ones(4,4))*ones(4,1);
%
subplot(8,1,[1 2])
stairs(Time, X(1:Tend, 1), 'b', 'LineWidth',1);
hold on
plot([0 0],[1.1 -0.1],'k:');
plot(Time,X(:,1),'bo');
plot(Tend-1, pbar(1), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(6.2, pbar(1), '$$\bar{p}_0$$');
hold off
xleer;
latexylabel('$$p_0\;$$');
rotateY;
yticks([0 1]);
axis([-0.4 Tend-1 -0.2 1.2]);
latextitle('Behaviour of the machine tool');
%
subplot(8,1,[3 4])
stairs(Time, X(1:Tend, 2), 'b', 'LineWidth',1);
hold on
plot([0 0],[1.1 -0.1],'k:');
plot(Time, X(:,2),'bo');
plot(Tend-1, pbar(2), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(6.2, pbar(2), '$$\bar{p}_1$$');
hold off
xleer;
latexylabel('$$p_1\;$$');
rotateY;
yticks([0 1]);
axis([-0.4 Tend-1 -0.2 1.2]);
%
subplot(8,1,[5 6])
stairs(Time, X(1:Tend, 3), 'b', 'LineWidth',1);
hold on
plot([0 0],[1.1 -0.1],'k:');
plot(Time, X(:,3),'bo');
plot(Tend-1, pbar(3), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(6.2, pbar(3), '$$\bar{p}_2$$');
hold off
xleer;
latexylabel('$$p_2\;$$');
rotateY;
yticks([0 1]);
axis([-0.4 Tend-1 -0.2 1.2]);
%
subplot(8,1,[7 8])
stairs(Time, X(1:Tend, 4), 'b', 'LineWidth',1);
hold on
plot([0 0],[1.1 -0.1],'k:');
plot(Time, X(:,4),'bo');
plot(Tend-1, pbar(4), 'bo', 'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latextext(6.2, pbar(4), '$$\bar{p}_3$$');
hold off
latexxlabel('$$k$$');
xticks([0 1 2 3 4 5 6]);
latexylabel('$$p_3\;$$');
rotateY;
yticks([0 1]);
axis([-0.4 Tend-1 -0.2 1.2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    ES: Queue.m für 3. Auflage: 26.1.2017
%    figure(4)
%   Gleichverteilte Ankunftszeiten bei lambda=1, mu=1.2
%%   Uniform distribution of the arriving and departing workpieces
%
figure(2)
Nmax=50;      % number of workpieces considered
Tmax=80;
alpha=1.4;
z0=1;         % initial state: z=0
N0=0;
T0=0;
% TA = arrival times
% TDelta = temporal distance of arriving workpieces
% TB  = service times
%
%Queue1=rng;               % save s to restore the generator setting
%save('Queue1', 'Queue1');
load('Queue1');       % Seed (initialisation) of the number generator
rng(Queue1);

TDelta=1/pW + alpha*rand(1,Nmax);
TB=1/pM + alpha*rand(1,Nmax);
TA=[];
TA(1)=0.5;
TW=[];
TW(1)=0;
for k1=1:Nmax-1
    TW(k1+1)=max(TW(k1)+TB(k1)-TDelta(k1), 0);
    TA(k1+1)=TA(k1)+TDelta(k1);
end
%
TD=TA+TW+TB;
subplot(6,1,[1 2])
queueplotblue(TA, TD);
hold on
plot([0 Tmax],[0 0], 'k:')
axis([0 Tmax -0.3 0.3]);
latextitle('$${\tt Queue:}$$ Machine tool with random arrival times');
xticksempty([0 20 40 60 80]);
xleer;
hold off
%
%  Bestimmung der Anzahl N von Kunden im Wartesystem
kA=0;
kD=0;
kAmax=length(TA);
kDmax=length(TD);
kmax=kAmax+kDmax;
N(1)=N0;
T(1)=T0;
k1=1;
while (kD<kDmax) & (kA<kAmax)
    k1=k1+1;
    if TA(kA+1) < TD(kD+1)
        N(k1)=N(k1-1)+1;
        T(k1)=TA(kA+1);
        kA=kA+1;
    else
        N(k1)=N(k1-1)-1;
        T(k1)=TD(kD+1);
        kD=kD+1;
    end
end       % while-Schleife verlassen, wenn entweder alle TA-Werte 
          % oder alle TD-Werte berücksichtigt sind
while kD<kDmax    % Auswertung der restlichen TD-Werte
    k1=k1+1;
    N(k1)=N(k1-1)-1;
    T(k1)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k1=k1+1;
    N(k1)=N(k1-1)+1;
    T(k1)=TD(kA+1);
    kA=kA+1;
end
subplot(6,1,[3 4 5])
stairs(T, N,'b', 'LineWidth',1.5);
hold on
plot([max(T) Tmax],[0 0]);
axis([0 Tmax -0.3 3.3]);
latexylabel('$$z\;$$');
xticks([0 20 40 60 80]);
rotateY;
yticks([0 1 2 3]);
latexxlabel('$$t$$ in s');
hold off






%%   Figures
%
epsfigc15('Queue');
