%
%      QUEUE.M
%
%  Simulation eines Wartesystems mit 1 Warteschlange und 1 Bedieneinheit
%
% Jan Lunze
% 3.4.2006
% für 2. Auflage: 20.6.2012
% für 3. Auflage: 26.1.2017
% für 4. Auflage: 12.1.2025
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  TA = Ankunftszeit des k-ten Kunden
%  TA0 = Ankunftszeit des ersten Kunden
%  TDelta = Zeitlicher Abstand der Ankunft des k+1-ten Kunden nach dem
%           k-ten Kunden
%  TW = Wartezeit des k-ten Kunden
%  TB = Bedienzeit des k-ten Kunden
%
%  lambda = Ankunftsrate
%  mu = Bedienrate
%
%  Nmax = Anzahl der betrachteten Kunden
%
close all
clear
%
echo on
%
%  Warteschlange mit deterministischen Zeiten bei lambda=1; mu=1,2
% 
echo off
lambda = 1;
mu = 1.2;
N0=0;                   % Warteschlangenlänge zur Zeit t=0
T0=0;                   % Beginn der Betrachtungen bei t=0
TA0=1.5;                  % Ankunft des ersten Kundens
TW(1)=0;                % da Wartesystem zu Beginn leer;
Nmax=9;     
Tmax=13;                % Horizont für grafische Darstellungen
TDelta=1/lambda*ones(1,Nmax);
TB=1/mu*ones(1,Nmax);
TA(1)=TA0;
for k=1:Nmax-1
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
    TA(k+1)=TA(k)+TDelta(k);
end
%
figure(1)
subplot(7,1,[6,7])
stem(TW, 'b', 'LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0]);
latexylabel('$$T_{\rm W}$$');  % '$$T_{\rm W}$$ in s'
rotateY;
xleer;
yticks([0 2 ]);
latexxlabel('$$k$$');
xticks([0 2 4 6 8 10]);
axis([0 Nmax+1 -0.05 2]);
hold off
%
%
TD=TA+TW+TB;
subplot(7,1,[1,2])
queueplot(TA, TD);
hold on
plot([0 Tmax],[0 0])
axis([0 Tmax -0.3 0.3]);
xleer;
latextitle('Wartesystem mit deterministischen Zeiten ($$\lambda=1$$, $$\mu=1,2$$)');
xticksempty([0 2 4 6 8 10 12]);
%latexxlabel('$$t$$ in s');
%
% 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;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=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
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TD(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N, 'b', 'LineWidth',1.5);
hold on
plot([max(T) Tmax],[0 0], 'b');
axis([0 Tmax -0.3 3.3]);
latexylabel('$$N$$');
xticks([0 2 4 6 8 10 12]);
rotateY;
yticks([0 2 4]);
latexxlabel('$$t$$ in s');
hold off


%%  Warteschlange mit deterministischen Zeiten bei lambda=1; mu=0,7
% 
figure(2)
lambda = 1;
mu = 0.7;
TB=1/mu*ones(1,Nmax);
for k=1:Nmax-1
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
    TA(k+1)=TA(k)+TDelta(k);
end
subplot(7,1,[6,7])
stem(TW, 'b', 'LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0],'b');
latexylabel('$$T_{\rm W}$$');  % '$$T_{\rm W}$$ in s'
rotateY;
yticks([0 3]);
latexxlabel('$$k$$');
xticks([0 2 4 6 8 10 12]);
axis([0 Nmax+1 -0.2 max(TW)+0.2]);
hold off
%
%
TD=TA+TW+TB;
subplot(7,1,[1,2])
queueplot(TA, TD);
xleer;
xticksempty([0 2 4 6 8 10 12]);
axis([0 Tmax -0.3 0.3]);
latextitle('Wartesystem mit deterministischen Zeiten ($$\lambda=1$$, $$\mu=0,7$$)');
%
%  Bestimmung der Anzahl N von Kunden im Wartesystem
%
kA=0;
kD=0;
kAmax=length(TA);
kDmax=length(TD);
kmax=kAmax+kDmax;
N=N0;
T=T0;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=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
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TD(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N,'b', 'LineWidth',1.5);
hold on
plot([0 Tmax],[0 0]);
axis([0 Tmax -0.3 max(N)+0.3]);
latexylabel('$$N$$');
xticks([0 2 4 6 8 10 12]);
rotateY;
yticks([0 4]);
latexxlabel('$$t$$ in s');
hold off


%%  Warteschlange mit deterministischen Zeiten bei lambda=1; mu=0,45
% 
figure(3)
lambda = 1;
mu = 0.45;
TB=1/mu*ones(1,Nmax);
N=N0;
T=T0;
for k=1:Nmax-1
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
    TA(k+1)=TA(k)+TDelta(k);
end
subplot(7,1,[6,7])
stem(TW, 'b', 'LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0]);
latexylabel('$$T_{\rm W}$$');  % '$$T_{\rm W}$$ in s'
rotateY;
yticks([0 8]);
latexxlabel('$$k$$');
xticks([0 2 4 6 8 10 12]);
axis([0 Nmax+1 -0.2 max(TW)+0.2]);
hold off
%
%
TD=TA+TW+TB;
subplot(7,1,[1,2])
queueplot(TA, TD);
xleer;
xticksempty([0 2 4 6 8 10 12]);
axis([0 Tmax -0.3 0.3]);
latextitle('Wartesystem mit deterministischen Zeiten ($$\lambda=1$$, $$\mu=0,5$$)');
%  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;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=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
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TD(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N,'b', 'LineWidth',1.5);
hold on
plot([0 Tmax],[0 0]);
axis([0 Tmax -0.3 max(N)+0.3]);
latexylabel('$$N$$');
rotateY;
yticks([0 5 10]);
latexxlabel('$$t$$ in s');
xticks([0 2 4 6 8 10 12]);
hold off


%%   Gleichverteilte Ankunftszeiten bei lambda=1, mu=1.2
%
figure(4)
%   Nmax=20 Kunden
%   lambda=1   ->  Tend=20;
Nmax=9;
Tend=Nmax/lambda;
load('RandQueueFig4');
rng(RandQueueFig4);
TA=Tend*rand(1,Nmax);
TA=sort(TA);
%   Bedienzeit konstant
mu = 1.2;
TB=1/mu*ones(1,Nmax);
N=N0;
T=T0;
TW=0;
for k=1:Nmax-1
    TDelta(k)=TA(k+1)-TA(k);
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
end
subplot(7,1,[6,7])
stem(TW, 'b', 'LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0],'b');
%    mittlere Wartezeit
TWquer=sum(TW)/length(TW);
Tmax=13;
plot([0 Tmax],[TWquer, TWquer],':');
latexylabel('$$T_{\rm W}$$');  % '$$T_{\rm W}$$ in s'
rotateY;
yticks([0 2]);
latexxlabel('$$k$$');
xticks([0 2 4 6 8 10 12]);
axis([0 Nmax+1 -0.05 2]);
hold off
%
%
TD=TA+TW+TB;
subplot(7,1,[1,2])
queueplot(TA, TD);
xleer;
xticksempty([0 2 4 6 8 10 12]);
hold on
plot([0 Tmax],[0 0])
axis([0 Tmax -0.3 0.3]);
latextitle('Wartesystem mit stochastischen Zeiten ($$\lambda=1$$, $$\mu=1,2$$)');
%  Bestimmung der Anzahl N von Kunden im Wartesystem
%
kA=0;
kD=0;
kAmax=length(TA);
kDmax=length(TD);
kmax=kAmax+kDmax;
N=[];
T=[];
N(1)=N0;
T(1)=T0;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=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
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TD(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N,'b', 'LineWidth',1.5);
hold on
%   mittlere Anzahl von Kunden im System
Nquer=tmean(N, T);
plot([0 Tmax],[Nquer, Nquer],':');
plot([max(T) Tmax],[0 0],'b');
axis([0 Tmax -0.3 3.3]);
latexylabel('$$N$$');
rotateY;
yticks([0 2]);
latexxlabel('$$t$$ in s');
xticks([0 2 4 6 8 10 12]);
hold off;



%%   Gleichverteilte Ankunftszeiten mit größerer Anzahl von Kunden
%
figure(5)
%   Nmax=20 Kunden
%   lambda=1   ->  Tend=20;
Nmax=100;
Tend=Nmax/lambda;
load('RandQueueFig5');
rng(RandQueueFig5);
TA=Tend*rand(1,Nmax);
TA=sort(TA);
%   Bedienzeit konstant
mu = 1.2;
TB=1/mu*ones(1,Nmax);
N=N0;
T=T0;
TW=0;
for k=1:Nmax-1
    TDelta(k)=TA(k+1)-TA(k);
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
end
subplot(7,1,[6,7])
stem(TW, 'b', 'LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0]);
%    mittlere Wartezeit
TWquer=sum(TW)/length(TW);
Tmax=Tend+ TWquer+TB(1);
plot([0 Tmax],[TWquer, TWquer],':');
latexylabel('$$T_{\rm W}$$');  % '$$T_{\rm W}$$ in s'
rotateY;
latexxlabel('$$k$$');
axis([0 Nmax+1 -0.2 max(TW)+0.2]);
yticks([0 8]);
hold off
%
%
TD=TA+TW+TB;
subplot(7,1,[1,2])
queueplot(TA, TD);
xleer;
xticksempty([0 20 40 60 80 100]);
axis([0 Tmax -0.3 0.3]);
latextitle('Wartesystem mit stochastischen Zeiten ($$\lambda=1$$, $$\mu=1,2$$)');
%  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;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=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
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TD(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N,'b', 'LineWidth',1.5);
hold on
%   mittlere Anzahl von Kunden im System
Nquer=tmean(N, T);
plot([0 Tmax],[Nquer, Nquer],':');
plot([0 Tmax],[0 0], 'b');
axis([0 Tmax -0.3 max(N)+0.3]);
latexylabel('$$N$$');
rotateY;
yticks([0 10]);
latexxlabel('$$t$$ in s');
xticks([0 20 40 60 80 100]);
hold off;
%   Überprüfung des Gesetzes von Little
NW=max(N-1,0);
NWquer=tmean(NW, T)
TWquer



%%  Warteschlange mit deterministischen Zeiten bei lambda=1; mu=1,2
%  Am Anfang befinden sich 3 Kunden in der Warteschlange
% 
figure(6)
clear TA TB TW TD;
lambda = 1;
mu = 1.4;
N0=4;                   % Warteschlangenlänge zur Zeit t=0
T0=0;                   % Beginn der Betrachtungen bei t=0
TA0=0.7;                  % Ankunft des ersten Kundens
TW(1)=(N0-1)/mu;                % da Warteschlange zu Beginn mit N0-1 Kunden belegt;
Nmax=15;                 % Anzahl der ankommenden Kunden
Nmax=Nmax+N0;           % Anzahl der betrachteten Kunden
Tmax=13;                % Horizont für grafische Darstellungen
TDelta=1/lambda*ones(1,Nmax);
TB=1/mu*ones(1,Nmax);
for k=1:N0
    TA(k)=-(N0-k)/lambda;
end
TA(N0+1)=TA0;
for k=N0+1:Nmax-1
    TA(k+1)=TA(k)+TDelta(k);
end
for k=1:N0
    TW(k)=abs(TA(k))+k/mu;
end
TW(N0+1)=N0/mu-TA0;
for k=N0+1:Nmax-1
    TW(k+1)=max(TW(k)+TB(k)-TDelta(k), 0);
end
subplot(7,1,[6,7])
stem(TW, 'b','LineWidth',1.5);
hold on
plot([0 Nmax+1],[0 0],'b');
latexylabel('$$T_{\rm W}$$');  % '$$T_{\rm W}$$ in s'
rotateY;
latexxlabel('$$k$$');
yticks([0 4]);
xticks([0 5 10 15 20]);
axis([0 Nmax+1 -0.05 4]);
hold off
%
%
for k=1:N0
  TD(k)=k*1/mu;
end
for k=N0+1:length(TA)
   TD(k)=TA(k)+TW(k)+TB(k);
end
subplot(7,1,[1,2])
queueplot(TA, TD);
xleer;
xticksempty([0 2 4 6 8 10 12]);
hold on
plot([0 Tmax],[0 0])
axis([0 Tmax -0.3 0.3]);
latextitle('Wartesystem mit deterministischen Zeiten ($$\lambda=1$$, $$\mu=1,4$$)');
%  Bestimmung der Anzahl N von Kunden im Wartesystem
%
kA=N0+1;
kD=0;
kAmax=length(TA);
kDmax=length(TD);
kmax=kAmax+kDmax;
%   Verschiebung der Zeiten, um ausschließlich mit positiven Zeitpunkten zu
%   rechnen
N(1)=N0;
T(1)=T0;
k=1;
while (kD<kDmax) & (kA<kAmax)
    k=k+1;
    if TA(kA+1) < TD(kD+1)
        N(k)=N(k-1)+1;
        T(k)=TA(kA+1);
        kA=kA+1;
    else
        N(k)=N(k-1)-1;
        T(k)=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
    k=k+1;
    N(k)=N(k-1)-1;
    T(k)=TD(kD+1);
    kD=kD+1;
end
while kA<kAmax    % Auswertung der restlichen TA-Werte
    k=k+1;
    N(k)=N(k-1)+1;
    T(k)=TD(kA+1);
    kA=kA+1;
end
subplot(7,1,[3,4])
stairs(T, N, 'b','LineWidth',1.5);
hold on
plot([max(T) Tmax],[0 0],'b');
axis([0 Tmax -0.3 4.3]);
latexylabel('$$N$$');
rotateY;
latexxlabel('$$t$$ in s');
yticks([0  4]);
hold off

%%  Warteschlangenlänge in Abhängigkeit von Verkehrsdichte
%
figure(7)
Rho=0:0.01:1;
k1=0;
Nquer=[];
for rho=Rho
    k1=k1+1;
    Nquer(k1)=rho/(1-rho);    %Mittlere Kundenanzahl in Wartesystem
end
subplot(3,1,[1,2])
plot(Rho, Nquer,'b');
latexxlabel('$$\rho$$');
xticks([0 0.2 0.4 0.6 0.8 1]);
latexylabel('$$\bar{N}$$');
rotateY;
yticks([0 20 40 60 80]);
axis([0 1 -1 50]);
hold on
N1=0.9/(1-0.9);
N2=0.95/(1-0.95);
plot(0.9, N1, 'bo', 0.95, N2, 'bo');
plot([0.9 0.9], [N1 0],'k:', [0.95 0.95], [N2 0],'k:');


%%  Friseursalon (Aufgabe)
%
figure(8)
lambda=0.117;
mu=0.146;
%  Wahrscheinlichkeit, nicht zu warten:
%  Prob(T_W = 0) = 1 - rho
C=[1 2 3 4 5 6];
P=[];
k1=0;
for c=C
    k1=k1+1;
    rho=lambda/mu/c;
    P(k1)=1-rho;
end
subplot(3,1,[1,2])
plot(C, P, 'bo');
hold on
plot([0 10],[0.8 0.8],'k:');
latexxlabel('$$c$$');
xticks([1 2 3 4 5 6]);
latexylabel('$${\rm Prob}(T_{\rm W}\le 0)$$');
rotateY;
yticks([0 0.5 0.8 ]);
axis([1 6 0 1]);



%%
epsfigc('Queue');