A general code

The aim of this session is to be able to simulate and compare the dynamics of a G/G/1 queue adopting a number of scheduling disciplines. The following code provides a basis for implementing and testing other scheduling disciplines.

knitr::opts_chunk$set(echo = TRUE)

# This function plot the average response time of a M/G/1 queue for various scheduling disciplines.
# As input, it takes the job service times (Service) and their mean value (ES)
plot_response_times <- function(Service, ES){

    debug=0;
    scheduling_disciplines = c(1:4); # 'FCFS'=1,'LCFS'=2,'RO-S'=3,'SRPT'=4
    num_arrival_rates=5;
    lambdas = (1/ES) * c(1:num_arrival_rates)/(num_arrival_rates+1) ;
    
    # Matrix of average response times for all disciplines and arrival rates
    R_time <- matrix(0, nrow = length(scheduling_disciplines), ncol = length(lambdas));
    rownames(R_time) <- c("FCFS","LCFS","RO-S","SRPT");
    colnames(R_time) <- paste("rho=", round(lambdas*ES, digits = 3), sep = "");
    
    cnt=0;
    for (lambda in lambdas)
    {
      cnt=cnt+1;
      # Let's assume a Poisson process coming in, otherwise change the next line
      Arrival= cumsum(rexp(n=N, rate = lambda)); # Arrival times
  
      for (scheduling_discipline in scheduling_disciplines)
      {
        t = 0; # simulation time
        Remaining = rep(N, x=NA);    # Remaining service times of each job
        Completion = rep(N, x=NA);   # Completion time of each job
        
        AvgJobs = 0;
        CurrentTask = NA;
        NextArrival = 1;
        while (TRUE) {
        
        if (debug==1){
            print("*********************");
            print(Arrival);
            print(Service);
            print(Remaining);
        }
      
            dtA = NA;
            dtC = NA;
            if(length(Arrival[Arrival>t])>0) { # if an arrival exists after t
                dtA = head(Arrival[Arrival>t],n=1)-t ; # time to next arrival
            }
            if(!is.na(CurrentTask)) { # if a task is running
                dtC = Remaining[CurrentTask]; # time to next completion
            }
            if(is.na(dtA) & is.na(dtC)) {
                break;
            } 
        
            dt = min(dtA,dtC,na.rm=T);
        
            # update system variables
            t = t + dt;
            AvgJobs=AvgJobs+dt*sum(!is.na(Remaining));
      
        if (debug==1){
            print(paste("Sim. time:", t));
            print(paste("# of jobs arrived: ", NextArrival));
        }
            # Job arrival
            
            if((NextArrival <=N) & (Arrival[NextArrival] == t)) {
                Remaining[NextArrival] = Service[NextArrival];
                NextArrival = NextArrival + 1;
            }
      
            # Job departure
            
            if(!is.na(CurrentTask)) {
                Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
                if(Remaining[CurrentTask] <= 0) {
                    # CurrentTask completed
                    Completion[CurrentTask] = t;
                    Remaining[CurrentTask] = NA;
                    CurrentTask = NA;
                }
            }
            
            # Scheduling discipline
            
            WaitingList=(1:NextArrival)[!is.na(Remaining)];
           
            if(length(WaitingList)>0) {
      
              # FCFS
              if (scheduling_discipline==1){
                CurrentTask = head(WaitingList,n=1);
              }
              
              # LCFS
              if (scheduling_discipline==2){
                # We implement the non-preemptive version of LCFS: jobs are never killed/resumed
                if(is.na(CurrentTask)) {
                # here, a task is not running, so we find a new job to serve
                  CurrentTask = tail(WaitingList,n=1);
                }
              }
              
              # ROS: Random Order of Service
              if (scheduling_discipline==3){
                if(is.na(CurrentTask)) {
                # here, a task is not running, so we find a new job to serve
                  if (length(WaitingList)>1) {
                    CurrentTask = sample(WaitingList,size=1);
                  } else{
                    CurrentTask = WaitingList;
                  }
                }
              }
              
              # SRPT: Shortest Remaining Processing Time
              if (scheduling_discipline==4){
                  CurrentTask = which.min(Remaining);
              }
      
            }
      
        if (debug==1){
            print(paste("Current Task = ", CurrentTask))
            readline(prompt="Press [enter] to proceed")
        }
            
        } # while
        
        ResponseTime = mean(Completion-Arrival); #average response time
        AvgJobs=AvgJobs/(tail(Completion,n=1)-Arrival[1]);
  
        # Simulation completed. Let's verify Little law: N=lambda*R
        r_names=rownames(R_time)
        cat(paste("Sim. ",r_names[scheduling_discipline],"(rho=",round(lambda*ES,digits = 2),") completed."));
        cat(paste(" Little law verification: ", round(AvgJobs,digits = 4), "=N=lambda*R=", round(lambda*ResponseTime,digits = 4),"\n"));
        
        R_time[scheduling_discipline,cnt]=ResponseTime;
  
      } # loop scheduling discipline
    
    } # loop lambda
    
    cat("\nResponse time matrix:\n")
    print(R_time)
    
    matplot(t(R_time), type = "l", xlab="Load (lambda*E[S])", ylab="Avg Response Time", xaxt = "n", lwd = 3);
    legend("topleft", legend = rownames(R_time), col=1:4, pch=1);
    axis(side=1,at=1:ncol(R_time),labels=round(lambdas*ES, digits=3));


} # end function

Let’s compare the various policies

First, let us consider exponentially distributed service times. We obtain:

knitr::opts_chunk$set(echo = FALSE)

set.seed(11); # Set seed for reproducibility
N = 1e4; # number of jobs to simulate;
mu=1;
Service=rexp(N,rate=mu); # Service times
plot_response_times(Service, 1/mu);
## Sim.  FCFS (rho= 0.17 ) completed. Little law verification:  0.1997 =N=lambda*R= 0.1988 
## Sim.  LCFS (rho= 0.17 ) completed. Little law verification:  0.1998 =N=lambda*R= 0.1989 
## Sim.  RO-S (rho= 0.17 ) completed. Little law verification:  0.2003 =N=lambda*R= 0.1994 
## Sim.  SRPT (rho= 0.17 ) completed. Little law verification:  0.1824 =N=lambda*R= 0.1816 
## Sim.  FCFS (rho= 0.33 ) completed. Little law verification:  0.5014 =N=lambda*R= 0.5034 
## Sim.  LCFS (rho= 0.33 ) completed. Little law verification:  0.5005 =N=lambda*R= 0.5024 
## Sim.  RO-S (rho= 0.33 ) completed. Little law verification:  0.5004 =N=lambda*R= 0.5023 
## Sim.  SRPT (rho= 0.33 ) completed. Little law verification:  0.4083 =N=lambda*R= 0.4099 
## Sim.  FCFS (rho= 0.5 ) completed. Little law verification:  1.003 =N=lambda*R= 0.9944 
## Sim.  LCFS (rho= 0.5 ) completed. Little law verification:  1.0059 =N=lambda*R= 0.9972 
## Sim.  RO-S (rho= 0.5 ) completed. Little law verification:  0.9979 =N=lambda*R= 0.9892 
## Sim.  SRPT (rho= 0.5 ) completed. Little law verification:  0.7154 =N=lambda*R= 0.7092 
## Sim.  FCFS (rho= 0.67 ) completed. Little law verification:  1.8961 =N=lambda*R= 1.9044 
## Sim.  LCFS (rho= 0.67 ) completed. Little law verification:  1.9013 =N=lambda*R= 1.9088 
## Sim.  RO-S (rho= 0.67 ) completed. Little law verification:  1.9034 =N=lambda*R= 1.911 
## Sim.  SRPT (rho= 0.67 ) completed. Little law verification:  1.1401 =N=lambda*R= 1.1446 
## Sim.  FCFS (rho= 0.83 ) completed. Little law verification:  4.1629 =N=lambda*R= 4.2065 
## Sim.  LCFS (rho= 0.83 ) completed. Little law verification:  4.2786 =N=lambda*R= 4.3234 
## Sim.  RO-S (rho= 0.83 ) completed. Little law verification:  4.208 =N=lambda*R= 4.2521 
## Sim.  SRPT (rho= 0.83 ) completed. Little law verification:  1.9805 =N=lambda*R= 2.0013 
## 
## Response time matrix:
##      rho=0.167 rho=0.333  rho=0.5 rho=0.667 rho=0.833
## FCFS  1.192658  1.510136 1.988710  2.856673  5.047847
## LCFS  1.193472  1.507208 1.994368  2.863176  5.188069
## RO-S  1.196229  1.507016 1.978496  2.866450  5.102473
## SRPT  1.089504  1.229592 1.418461  1.716971  2.401556

Then, let us assume that service times follows a Pareto\((x_m,\alpha)\) distribution where \(x_m\) and \(\alpha\) are the scale and shape parameters, respectively. Let us assume \(\alpha=3\) and \(x_m=1\). We obtain:

## Sim.  FCFS (rho= 0.17 ) completed. Little law verification:  0.1851 =N=lambda*R= 0.1853 
## Sim.  LCFS (rho= 0.17 ) completed. Little law verification:  0.1852 =N=lambda*R= 0.1854 
## Sim.  RO-S (rho= 0.17 ) completed. Little law verification:  0.185 =N=lambda*R= 0.1852 
## Sim.  SRPT (rho= 0.17 ) completed. Little law verification:  0.1802 =N=lambda*R= 0.1804 
## Sim.  FCFS (rho= 0.33 ) completed. Little law verification:  0.4336 =N=lambda*R= 0.4338 
## Sim.  LCFS (rho= 0.33 ) completed. Little law verification:  0.4337 =N=lambda*R= 0.4339 
## Sim.  RO-S (rho= 0.33 ) completed. Little law verification:  0.4341 =N=lambda*R= 0.4343 
## Sim.  SRPT (rho= 0.33 ) completed. Little law verification:  0.404 =N=lambda*R= 0.4041 
## Sim.  FCFS (rho= 0.5 ) completed. Little law verification:  0.8372 =N=lambda*R= 0.8284 
## Sim.  LCFS (rho= 0.5 ) completed. Little law verification:  0.832 =N=lambda*R= 0.8233 
## Sim.  RO-S (rho= 0.5 ) completed. Little law verification:  0.836 =N=lambda*R= 0.8272 
## Sim.  SRPT (rho= 0.5 ) completed. Little law verification:  0.7265 =N=lambda*R= 0.7189 
## Sim.  FCFS (rho= 0.67 ) completed. Little law verification:  1.5386 =N=lambda*R= 1.5321 
## Sim.  LCFS (rho= 0.67 ) completed. Little law verification:  1.5424 =N=lambda*R= 1.5357 
## Sim.  RO-S (rho= 0.67 ) completed. Little law verification:  1.5415 =N=lambda*R= 1.535 
## Sim.  SRPT (rho= 0.67 ) completed. Little law verification:  1.218 =N=lambda*R= 1.2127 
## Sim.  FCFS (rho= 0.83 ) completed. Little law verification:  3.0311 =N=lambda*R= 3.0524 
## Sim.  LCFS (rho= 0.83 ) completed. Little law verification:  3.0401 =N=lambda*R= 3.0598 
## Sim.  RO-S (rho= 0.83 ) completed. Little law verification:  3.0098 =N=lambda*R= 3.0293 
## Sim.  SRPT (rho= 0.83 ) completed. Little law verification:  2.0524 =N=lambda*R= 2.0657 
## 
## Response time matrix:
##      rho=0.167 rho=0.333  rho=0.5 rho=0.667 rho=0.833
## FCFS  1.667337  1.952046 2.485255  3.447114  5.494236
## LCFS  1.668261  1.952616 2.469828  3.455425  5.507670
## RO-S  1.666833  1.954130 2.481676  3.453706  5.452779
## SRPT  1.623788  1.818601 2.156574  2.728548  3.718218

From the plots, we see that FCFS, LCFS and ROS provide the same average response time. SRPT seems to be the best, and this is indeed the case in a broad sense [1]. However, SRPT remains a bit ideal because one does not really know how much processing time remains for any given job. To patch this, a possibility would be to consider SERPT, i.e., Shortest Expected Remaining Processing Time, which has an obvious meaning; this policy may make sense in practice because the service time distribution is known, as it can be learned from the data. Or, would it be possible to do even better than SERPT?

References

[1] Schrage, Linus. “A Proof of the Optimality of the Shortest Remaining Processing Time Discipline.” Operations Research, vol. 16, no. 3, 1968, pp. 687–90.