%%
% M/M/s, Matlab code, and Example 
%    
% 

% input parameters
s=1;lmd=1/2;mu=2/3;

%%%%%%% general formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq)

Lq=P0*seq(end)*rho/(1-rho)
L=Lq+r
Wq=Lq/lmd  % in minutes
W=Wq+1/mu 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
% M/M/s, Lecture Example, Preemptive Priority 
%

% First scenario: One doctor, s=1:

s=1;

lmd0=2;mu=3;

p1=.1;p2=.3;p3=.6;

%%%% arrival rate for critical case
lmd=p1*lmd0;

%%%%%%% applying the formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq);

Lq=P0*seq(end)*rho/(1-rho);
Wq=Lq/lmd;
W=Wq+1/mu;
L=Lq+r;

% times for the critical case:
W1=W;
W1q=Wq;


%%%% arrival rate for critical and serious cases
lmd=(p1+p2)*lmd0;

%%%%%%% applying the formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq);

Lq=P0*seq(end)*rho/(1-rho);
Wq=Lq/lmd;
W=Wq+1/mu;
L=Lq+r;

W1_2=W; % W1_2=p1/(p1+p2)*W1+p2/(p1+p2)*W2
W1_2q=Wq;

% times for the serious case:
W2=(W1_2-p1/(p1+p2)*W1)*(p1+p2)/p2;
W2q=(W1_2q-p1/(p1+p2)*W1q)*(p1+p2)/p2;


%%%% arrival rate for all cases
lmd=(p1+p2+p3)*lmd0;

%%%%%%% applying the formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq);

Lq=P0*seq(end)*rho/(1-rho);
Wq=Lq/lmd;
W=Wq+1/mu;
L=Lq+r;

W1_3=W; % W1_3=p1/(p1+p2+p3)*W1+p2/(p1+p2+p3)*W2+p3/(p1+p2+p3)*W3
W1_3q=Wq;

% times for the stable case:
W3=(W1_3-p1/(p1+p2+p3)*W1-p2/(p1+p2+p3)*W2)*(p1+p2+p3)/p3;
W3q=(W1_3q-p1/(p1+p2+p3)*W1q-p2/(p1+p2+p3)*W2q)*(p1+p2+p3)/p3;

%%%%%%% list the times in minutes together for all cases: 
W=[W1,W2,W3]*60
Wq=[W1q,W2q,W3q]*60

%%
% Second scenario: One doctor, s=2, evreything else the same:

s=2;

lmd0=2;mu=3;

p1=.1;p2=.3;p3=.6;

%%%% arrival rate for critical case
lmd=p1*lmd0;

%%%%%%% applying the formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq);

Lq=P0*seq(end)*rho/(1-rho);
Wq=Lq/lmd;
W=Wq+1/mu;
L=Lq+r;

% times for the critical case:
W1=W;
W1q=Wq;


%%%% arrival rate for critical and serious cases
lmd=(p1+p2)*lmd0;

%%%%%%% applying the formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq);

Lq=P0*seq(end)*rho/(1-rho);
Wq=Lq/lmd;
W=Wq+1/mu;
L=Lq+r;

W1_2=W; % W1_2=p1/(p1+p2)*W1+p2/(p1+p2)*W2
W1_2q=Wq;

% times for the serious case:
W2=(W1_2-p1/(p1+p2)*W1)*(p1+p2)/p2;
W2q=(W1_2q-p1/(p1+p2)*W1q)*(p1+p2)/p2;


%%%% arrival rate for all cases
lmd=(p1+p2+p3)*lmd0;

%%%%%%% applying the formula %%%%%%%%%
r=lmd/mu;
rho=lmd/mu/s;
serverseq=0:s;
seq=r.^serverseq./factorial(serverseq);
seq(end)=seq(end)/(1-rho);

P0=1/sum(seq);

Lq=P0*seq(end)*rho/(1-rho);
Wq=Lq/lmd;
W=Wq+1/mu;
L=Lq+r;

W1_3=W; % W1_3=p1/(p1+p2+p3)*W1+p2/(p1+p2+p3)*W2+p3/(p1+p2+p3)*W3
W1_3q=Wq;

% times for the stable case:
W3=(W1_3-p1/(p1+p2+p3)*W1-p2/(p1+p2+p3)*W2)*(p1+p2+p3)/p3;
W3q=(W1_3q-p1/(p1+p2+p3)*W1q-p2/(p1+p2+p3)*W2q)*(p1+p2+p3)/p3;

%%%%%%% list the times in minutes together for all cases: 
W=[W1,W2,W3]*60
Wq=[W1q,W2q,W3q]*60


%%
%
% Non-preemptive MMs model
%

s=3;
mu=20;
lmdseq=[5,15,30];
N=length(lmdseq);

lmd=sum(lmdseq); 
r=lmd/mu;

sseq=0:(s-1);
A=factorial(s)*(s*mu-lmd)/r^s*sum(r.^sseq./factorial(sseq))+s*mu;

B=zeros(1,N);
for i=1:N
    B(i)=1-sum(lmdseq(1:i))/s/mu;
end
Bk_1=[1,B(1:end-1)];
Bk=B;

Wqseq=1./(A.*Bk_1.*Bk);
Wseq=Wqseq+1/mu;

Lseq=lmdseq.*Wseq;
Lqseq=lmdseq.*Wqseq;

Wq=Wqseq*60
W=Wseq*60


