%Matlab Program MMsrMa.m
Program = 'MMsrMa:  Calculations for M/M/s/r+M following new call-center paper'
WhoWhen = 'Ward Whitt 11/27/03'

%For exponential (M) abandon-time distribution

%To calculate waiting distributions, we use numerical transform inversion,
% using the functions inversion2.m and laplace2.m

Step1 = 'model parameters'

%s = number of servers
 s = 'number of servers'
   s=100
%lambda = arrival rate
lambda = 'arrival rate' 
   lambda=100
%mu = individual service rate
 mu = 'individual service rate'
   mu=1.0
%alpha = individual abandonment rate (exponential abandonments, not using)
  %alpha=1.0
%beta is reneging parameter, not used
   beta=0.0;
%r = number of waiting spaces
r = 'number of extra waiting spaces'
   r=100
%erlang E2 abandon-time distribution has mean M 
M = 'mean time to abandon'
M = 1

   % setting parameters 
   x=zeros(s+r,1);
   q=zeros(r,1);  
   delta=zeros(r,1);
   deltaInc=zeros(r,1);
    qnot=1;
    % Here is code that needs to be changed in generalization    
    for k = 1:r
     delta(k,1) = k/M;
 end
    %Code to be changed ends here
     deltaInc(1,1) = delta(1,1);
     q(1,1) =1;
     for k=2:r
         q(k,1)=1;
         deltaInc(k,1)=delta(k,1)-delta(k-1,1);
      end
   
   
   % calculating long run probabilities
   x(s,1)=1;
   x(s+1,1)=lambda*(1-beta)*qnot*x(s,1)/(s*mu+delta(1,1));
   for k=1:r-1
       x(s+k+1,1)=lambda*(1-beta)*q(k,1)*x(s+k,1)/(s*mu+delta(k+1,1));
    end
    for k=1:s-1
       x(s-k,1)=(s-k+1)*mu*x(s-k+1,1)/lambda;
    end
    xnot=mu*x(1,1)/lambda;
    y=xnot+sum(x);
    xnot=xnot/y;
    for k=1:s+r
       x(k,1)=x(k,1)/y;
    end
    
   
    Step2 = 'Key Steady-state Probabilities'
    
     Ps1=xnot;
    for k=1:s-1
       Ps1=x(k,1)+Ps1;
    end
    %Probability of no wait = Ps1
    ProbNoWait = Ps1
    ProbAllBusy = 1 - ProbNoWait     
   
   %littlegama = probabilities of next departure event, \gamma_{k,j} in paper
    littlegama=zeros(r,r); 
    for k=1:r
        littlegama(k,1) = deltaInc(1,1)/(s*mu+delta(k,1));
        for j=2:k            
        littlegama(k,j) = deltaInc(j,1)/(s*mu+(delta(k,1)-delta(j-1,1)));
    end        
end
    %gama = Big Gamma = probabilities of receiving service, \Gamma_{k} in paper
    gama = zeros(r,1);     
    for k=1:r
        gama(k,1) = 1 - littlegama(k,1);
        for j=2:k            
        gama(k,1) = gama(k,1)*(1-littlegama(k,j));
    end        
end
 %The means appear later, starting in the response time moments.  
  
 %now we compute the probs of finding k and getting served = p2(k,1)  
 %p2not is the probability of finding just s in the system, none in queue
    p2not=x(s,1)*(1-beta)*qnot*gama(1,1);
    p2=zeros(r-1,1);
    for k=1:r-1
        p2(k,1)=x(s+k,1)*(1-beta)*q(k,1)*gama(k+1,1);
    end
    
 
 
    
    
    
    Ps2=p2not+sum(p2);
    % Prob Wait and Served = Ps2
    %probability customer is served = Ps
    ProbWaitandServed = Ps2
    
    Ps=Ps1+Ps2;
    ProbServed = Ps
    %Probability abandon = 1 - Ps
        ProbLoss = x(s+r)
    ProbAban = 1 - Ps1 - Ps2 - x(s+r)
    
    
    
   Step3 = 'Mean Values - Counting'
    mq=0.0;
    for k=1:r
        mq=mq+k*x(s+k);
    end
    %Mean number in Queue = mq
    MeanInQueue = mq
    
    Pbus = 0.0;
    for k=0:r
        Pbus=Pbus +x(s+k);
    end
    %Probability all servers busy = Pbus (a check)
    ProbAllBusy = Pbus;
    
    Mbus = 0.0;
    for k=1:s-1
        Mbus=Mbus + k*x(k);
    end
    Mbus = Mbus + Pbus*s;
    %Mean Number of Busy Servers = Mbus
    MeanBusyServers = Mbus
    
    %Mean Number in System = Mnum
    Mnum = Mbus + mq;
    MeanNumberinSys = Mnum
    
    
    Step4 = 'Second Moments and Variances - Counting'
     m2q=0.0;
    for k=1:r
        m2q=m2q+(k^2)*x(s+k);
    end
    %Second moment number in Queue = mq
    SecondMomInQueue = m2q 
    VarInQueue = m2q - mq^2
    SDinQueue = sqrt(VarInQueue)
       
    M2bus = 0.0;
    for k=1:s-1
        M2bus=M2bus + (k^2)*x(k);
    end
    M2bus = M2bus + Pbus*(s^2);
    %Second moment of Number of Busy Servers = M2bus
    SecondMomBusyServers = M2bus
    VarBusyServers = M2bus - Mbus^2
    
    %Second Moment Number in System = M2num
    M2num1 = 0;
    M2num2 = 0;
    for k=1:s
        M2num1 = M2num1 + (k^2)*x(k);
    end
    for k=0:r-1
        M2num2 = M2num2 + ((s+r-k)^2)*x(s+r-k);
    end
    M2num = M2num1 + M2num2;
    SecondMomNumberInSys = M2num
    VarNumberInSys = M2num - Mnum^2
    SDNumberInSys = sqrt(VarNumberInSys)
    
    %calculations for C
   Step5 = 'Response-Time Moments'    
    
   % First the means
   %mndpevt mean time to departures events, m_{k,j} in paper
   mndpevt=zeros(r,r);  
    for k=1:r
        mndpevt(k,1) = 1/(s*mu+delta(k,1));
        for j=2:k            
        mndpevt(k,j) = 1/(s*mu+(delta(k,1)-delta(j-1,1)));
    end 
end
        summ=mndpevt(1,1);
    Ec=Ps/mu+p2not*summ;
    msum=zeros(r,r);
    for k=1:r
        msum(k,1) = mndpevt(k,1);
        for j=2:k
       msum(k,j)=msum(k,j-1)+mndpevt(k,j);
   end
end
for k=1:r-1
       Ec=Ec+p2(k)*msum(k+1,k+1);
    end
    %mean response time and conditional mean
    Ec;
    ECIs=Ec/Ps;
    MeanRespTime = Ec
    ConditMeanRespTimeServed = ECIs
    
    %Next the second moments and variances
    
    
    MM1=zeros(r,1);
    MM2=zeros(r,1); 
    MM1(1,1) = mndpevt(1,1);
    MM2(1,1) = mndpevt(1,1)^2;
    for k=2:r
        MM1(k,1) = (1/mu) + mndpevt(k,1);
         MM2(k,1) = (1/mu)^2 + mndpevt(k,1)^2;
        for j=2:k
        MM1(k,1)= MM1(k,1) + mndpevt(k,j);
        MM2(k,1)=MM2(k,1) + mndpevt(k,j)^2;
    end
    end
    
    Ec2=Ps1*(2/mu^2);
    Ec2 = Ec2 + p2not*(MM2(1,1)+ MM1(1,1)^2);
    
    for k=1:r-1
       Ec2=Ec2+p2(k,1)*(MM2(k+1,1)+ MM1(k+1,1)^2);
    end
    SecMomRespTime = Ec2
    VarRespTime = Ec2 - Ec^2
    ConditSecMomRespTimeServed = Ec2/Ps
    ConditVarRespTimeServed = Ec2/Ps - ECIs^2
    
    Step6 = 'Waiting-Time Moments'
    
       
    Step6a = 'Waiting Times for Customers who are Served'
    
    %First the means
     summ=mndpevt(1,1);
     EWandServed = p2not*summ;        
    for k=1:r-1
        
   EWandServed = EWandServed + p2(k)*msum(k+1,k+1);
    end
     %mean waiting time and served, and conditional mean
   EWandServed;
    CondEWserved=EWandServed/Ps;
    MeanWaitServed = EWandServed
    ConditMeanWaitServed = CondEWserved
    
    %Now second moments and variances
    
    MM1w=zeros(r,1);
    MM2w=zeros(r,1); 
    MM1w(1,1) = mndpevt(1,1);
    MM2w(1,1) = mndpevt(1,1)^2;
    for k=2:r
        MM1w(k,1) =  mndpevt(k,1);
         MM2w(k,1) =  mndpevt(k,1)^2;
        for j=2:k
        MM1w(k,1)= MM1w(k,1) + mndpevt(k,j);
        MM2w(k,1)=MM2w(k,1) + mndpevt(k,j)^2;
    end
    end
    
     EWandServed2 = p2not*(MM2w(1,1)+ MM1w(1,1)^2);    
    for k=1:r-1
       EWandServed2=EWandServed2+p2(k,1)*(MM2w(k+1,1)+ MM1w(k+1,1)^2);
    end
    SecMomWaitandServed = EWandServed2
    VarWaitandServed = EWandServed2 - EWandServed^2
    ConditSecMomWaitServed = EWandServed2/Ps
    CondEWserved2 = EWandServed2/Ps;
    ConditVarWaitServed = EWandServed2/Ps - CondEWserved^2
    SDcondVarWaitServ = sqrt(ConditVarWaitServed)
    
           
    Step6b = 'Waiting Times for Customers who Abandon'
    
    %Product of (1 - gamma(k,j)) in Section 4.8 of the paper   
    
    N = zeros(r,r);
    for k = 1:r        
    N(k,1) = 1 - littlegama(k,1);
    for j=2:k
        N(k,j)= N(k,j-1)*(1 - littlegama(k,j));
    end
end

    MA1 = zeros(r,1);
    for k = 1:r
    MA1(k,1) = littlegama(k,1)*mndpevt(k,1);
    for j=2:k
        MA1(k,1) = MA1(k,1) + N(k,j-1)*littlegama(k,j)*msum(k,j);
    end
end    
    
%analog of msum for second moments

  msum2=zeros(r,r);
    for k=1:r
        msum2(k,1) = mndpevt(k,1)^2;
        for j=2:k
       msum2(k,j)=msum2(k,j-1)+mndpevt(k,j)^2;
   end
end

     MA2 = zeros(r,1);
    for k = 1:r
    MA2(k,1) = littlegama(k,1)*2*mndpevt(k,1)^2;
    for j=2:k
        MA2(k,1) = MA2(k,1) + N(k,j-1)*littlegama(k,j)*(msum(k,j)^2 + msum2(k,j));
    end
end    
  
   
    EA = x(s)*MA1(1,1);
    for k=1:r-1
        EA=EA + x(s+k)*MA1(k+1,1);
    end
    
    MeanWaitAbandon = EA
    ConditMeanWaitAban = EA/ProbAban
    
    %Second Moments of Wait for Customers Who Abandon    
    
    EA2 = x(s)*MA2(1,1);
    for k=1:r-1
        EA2 = EA2 + x(s+k)*MA2(k+1,1);
    end
    
    SecMomWaitAbandon = EA2
    ConditSecMomWaitAbandon = EA2/ProbAban
    ConditVarWaitAban = EA2/ProbAban - (EA/ProbAban)^2
    SDcondVarWaitAban = sqrt(ConditVarWaitAban)    
        
    Step6c = 'Waiting Time for All Customers, Served or Abandoning'
    
    
    EWall = ProbAban*(EA/ProbAban) + Ps*CondEWserved;
    EWall2 = ProbAban*(EA2/ProbAban) + Ps*CondEWserved2;
    
    MeanWaitTimeAll = EWall   
    SecMomWaitTimeAll = EWall2
    VarWaitTimeAll = EWall2 - EWall^2
    
    
    
    
    
    Step7 = 'Conditional Waiting Time CDF Values For Served Customers'
    
POKS = zeros(4,1);
          for where=1:4
      if where==1
          ar=0.05;
       elseif where==2
          ar=0.224;
       elseif where==3
          ar=0.1;
       else
          ar=0.2;
       end
      
       t=ar
       xxx=inversion2(5,t,Ps1,Ps2,delta,mu,p2not,p2,r,s,x,q,xnot,qnot,beta,mndpevt,littlegama);       
       POKS(where,1) = (Ps1 +  Ps2*xxx)/Ps;
      ProbOKWaitifServed = POKS(where,1)
  end
      
   Step8 = 'Conditional Waiting Time CDF Values For Customers Who Abandon'

   POKA = zeros(4,1);
          for where=1:4
      if where==1
          ar=0.05;
       elseif where==2
          ar=0.224;
       elseif where==3
          ar=0.1;
       else
          ar=0.2;
       end
      
       t=ar
       yyy=inversion2(2,t,Ps1,Ps2,delta,mu,p2not,p2,r,s,x,q,xnot,qnot,beta,mndpevt,littlegama);       
       POKA(where,1) = yyy/ProbAban;
       ProbOKWaitifAbandon = yyy/ProbAban
   end    
       Step9 = 'Waiting Time CDF Values For All Customers'

          for where=1:4
      if where==1
          ar=0.05;
       elseif where==2
          ar=0.224;
       elseif where==3
          ar=0.1;
       else
          ar=0.2;
       end
      
       t=ar
       ProbOKWait = ProbAban*POKA(where,1) + Ps*POKS(where,1)
    end