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 calculate the mean response time for a M/M/k system
calculate_theory_mean_time <- function(lambda, mu, k) {
  # Calculate the charge C(k, lambda, mu) with Erlang's C formula
  charge = (lambda / mu)^k / (factorial(k) * sum((lambda / mu)^(0:(k-1)) / factorial(0:(k-1))));
  # Calculate the mean response time
  mean_response_time = (charge / (k * mu - lambda)) + (1 / mu);
  return(mean_response_time);
}

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

    debug=0;
    scheduling_disciplines = c(1:4); # 'FCFS'=1,'LCFS'=2,'RO-S'=3,'SRPT'=4
    num_arrival_rates=5; #Q2
    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 = "");
    
    result <- list("Response_time" = R_time['FCFS',], "Load" = lambdas * ES, "colnames" = colnames(R_time));
    
    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 of each job/ Arrival[k] = 5s -> job k arrives at 5s
  
      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;
        Servers = rep(k, x=NA); # Create a vector of size k. Servers[k] = 5 -> le job 5 run dans le server k
        # Servers start index: 1 | stop index: k
        NextArrival = 1; # The job id of the next job
        while (TRUE) {
        
        if (debug==1){
            print("*********************");
            print(Arrival);
            print(Service);
            print(Remaining);
        }
      
            dtA = NA; #Time before next arrival
            dtC = NA; #Time before a task end
            if(length(Arrival[Arrival>t])>0) { # if an arrival exists after t
                dtA = head(Arrival[Arrival>t],n=1)-t ; # time to next arrival
            }
            
            # Filter to get only queue where there is a job running
            CurrentRunningServers = Servers[!is.na(Servers)];
            if(length(CurrentRunningServers) > 0){ # If there is at least one task running
              dtC = min(Remaining[Servers[!is.na(Servers)]]); #The job with the lower remaining time will be the first to be completed
            }
            
            if(is.na(dtA) & is.na(dtC)) {
                break;
            }
        
            dt = min(dtA,dtC,na.rm=T); #temps avant prochain event (sois fin completion sois arrivr d'un job)
        
            # update system variables
            t = t + dt; #update time
            AvgJobs=AvgJobs+dt*sum(!is.na(Remaining)); #udpdate average job variable
      
        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]; # On met dans les jobs qui reste a exec le temps ki fo                                                                # pour exec ce job
                NextArrival = NextArrival + 1;
            }
      
            # Job departure
            for (current_server_index in which(!is.na(Servers))){
              current_running_job_index = Servers[current_server_index];
              Remaining[current_running_job_index] = Remaining[current_running_job_index] - dt ;
              if(Remaining[current_running_job_index] <= 0) {
                  # Current job task is completed
                  Completion[current_running_job_index] = t;
                  Remaining[current_running_job_index] = NA;
                  Servers[current_server_index] = NA;
              }
            }
            
            # Scheduling discipline
            # Preemptive version is implemented so running tasks are removed from servers.
            Servers = rep(k, x=NA);
            
            # For every queue
            for(current_server_index in 1:k){
              WaitingList=(1:NextArrival)[!is.na(Remaining)];
              # Keep only not already present in Servers vector
              WaitingList=WaitingList[!(WaitingList %in% Servers)];
              if(length(WaitingList)>0) {
                
                # FCFS
                if (scheduling_discipline==1){
                  Servers[current_server_index] = head(WaitingList,n=1);
                }
                
                # LCFS
                if (scheduling_discipline==2){
                    Servers[current_server_index] = tail(WaitingList,n=1);
                }
                
                # ROS: Random Order of Service
                if (scheduling_discipline==3){
                  if(is.na(Servers[current_server_index])) {
                  # here, a task is not running, so we find a new job to serve
                    if (length(WaitingList)>1) {
                      Servers[current_server_index] = sample(WaitingList,size=1);
                    } else{
                      Servers[current_server_index] = WaitingList;
                    }
                  }
                }
                
                # SRPT: Shortest Remaining Processing Time
                if (scheduling_discipline==4){
                  remaining_not_scheduled_tasks = Remaining;
                  for(already_scheduled_server_index in Servers){
                    remaining_not_scheduled_tasks[already_scheduled_server_index] = NA;
                  }
                  Servers[current_server_index] = which.min(remaining_not_scheduled_tasks);
                }
        
              }
        
              if (debug==1){
                  print(paste("Current Task = ", Server))
                  readline(prompt="Press [enter] to proceed")
              }
          } #END For (scheduling)
        } #END WHILE(true)
        
        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)
        if(log_level == 1){
          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
    
    if(log_level == 1){
      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));
    }
    
    result <- list("Response_time" = R_time['FCFS',], "Load" = lambdas * ES, "colnames" = colnames(R_time));
    return(result);


} # 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 = 1e3; # number of jobs to simulate;
mu=1;
Service=rexp(N,rate=mu); # Service times


#Question 2 DOES NOT WORK because of calculate_theory_mean_time() function
#num_arrival_rates=5; #Q2
#ES = 1/mu;
#lambda = (1/ES) * num_arrival_rates/(num_arrival_rates+1);

#nb_server_to_test = 10;
#Simulation_mean_response_time = rep(nb_server_to_test, x=NA);
#Theory_mean_response_time = rep(nb_server_to_test, x=NA);

#for(index in 1:nb_server_to_test){
#  simulation_result = plot_response_timesGGK(Service, 1/mu, nb_server_to_test, 0);
#  Simulation_mean_response_time[index] = simulation_result$Response_time;
#  Theory_mean_response_time[index] = calculate_theory_mean_time(lambda, Service, nb_server_to_test);
#}
#print(Simulation_mean_response_time);
#print(Theory_mean_response_time);
#matplot(Theory_mean_response_time, type = "l", xlab="Nb jobs", ylab="Mean Response time", xaxt = "n", lwd = 3);
#abline(h = Simulation_mean_response_time, col="red", ylab="W*");
#abline(v = k, col="red");

#Question 3
Response_time_GGK = plot_response_timesGGK(Service, 1/mu, 3, 0);
Response_time_GG1 = plot_response_timesGGK(Service, 3/mu, 1, 0);
Response_time_ratio = (Response_time_GG1$Response_time)/(Response_time_GGK$Response_time);
print(Response_time_ratio);
## rho=0.167 rho=0.333   rho=0.5 rho=0.667 rho=0.833 
##  1.061272  1.118318  1.187667  1.254994  1.331614
matplot(Response_time_ratio, type = "l", xlab="Load (lambda*E[S])", ylab="GG1/GGK Response time ratio", xaxt = "n", lwd = 3);
legend("topleft", legend = paste("Ratio"), col=1:1, pch=1);
axis(side=1,at=1:length(Response_time_GGK$colnames),labels=round(Response_time_GGK$Load, digits=3));

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.