DM1 simulation of G/G/k queues
DM 1 EP - Simulation of G/G/k queues
This code is based on the G/G/1 queues simulation code for different scheduling disciplines, adapting it for G/G/k queues. Code realized in collaboration with Alexandre Moua
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/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){
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 = "");
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; # Le num du job qui est le prochain à rentrer
while (TRUE) {
if (debug==1){
print("*********************");
print(Arrival);
print(Service);
print(Remaining);
}
dtA = NA; #Time du prochain arriver
dtC = NA; #Time avant la prochaine fin de run de la task
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)
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")
result <- list("Response_time" = R_time['FCFS',], "Load" = lambdas * ES, "colnames" = colnames(R_time));
return(result);
} # end functionLet’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;
k = 3;
Service=rexp(N,rate=mu); # Service times
Response_time_GGK = plot_response_timesGGK(Service, 1/mu, k);## Sim. FCFS (rho= 0.17 ) completed. Little law verification: 0.1699 =N=lambda*R= 0.167
## Sim. LCFS (rho= 0.17 ) completed. Little law verification: 0.1699 =N=lambda*R= 0.167
## Sim. RO-S (rho= 0.17 ) completed. Little law verification: 0.1699 =N=lambda*R= 0.167
## Sim. SRPT (rho= 0.17 ) completed. Little law verification: 0.1699 =N=lambda*R= 0.167
## Sim. FCFS (rho= 0.33 ) completed. Little law verification: 0.3388 =N=lambda*R= 0.3345
## Sim. LCFS (rho= 0.33 ) completed. Little law verification: 0.3389 =N=lambda*R= 0.3346
## Sim. RO-S (rho= 0.33 ) completed. Little law verification: 0.3391 =N=lambda*R= 0.3348
## Sim. SRPT (rho= 0.33 ) completed. Little law verification: 0.3388 =N=lambda*R= 0.3345
## Sim. FCFS (rho= 0.5 ) completed. Little law verification: 0.5142 =N=lambda*R= 0.5081
## Sim. LCFS (rho= 0.5 ) completed. Little law verification: 0.5115 =N=lambda*R= 0.5054
## Sim. RO-S (rho= 0.5 ) completed. Little law verification: 0.5124 =N=lambda*R= 0.5063
## Sim. SRPT (rho= 0.5 ) completed. Little law verification: 0.5106 =N=lambda*R= 0.5045
## Sim. FCFS (rho= 0.67 ) completed. Little law verification: 0.6938 =N=lambda*R= 0.671
## Sim. LCFS (rho= 0.67 ) completed. Little law verification: 0.694 =N=lambda*R= 0.6712
## Sim. RO-S (rho= 0.67 ) completed. Little law verification: 0.6949 =N=lambda*R= 0.6721
## Sim. SRPT (rho= 0.67 ) completed. Little law verification: 0.6932 =N=lambda*R= 0.6704
## Sim. FCFS (rho= 0.83 ) completed. Little law verification: 0.8511 =N=lambda*R= 0.8654
## Sim. LCFS (rho= 0.83 ) completed. Little law verification: 0.8509 =N=lambda*R= 0.8653
## Sim. RO-S (rho= 0.83 ) completed. Little law verification: 0.8509 =N=lambda*R= 0.8653
## Sim. SRPT (rho= 0.83 ) completed. Little law verification: 0.8367 =N=lambda*R= 0.8509
##
## Response time matrix:Response_time_GG1 = plot_response_timesGGK(Service, 3/mu, 1);## Sim. FCFS (rho= 0.17 ) completed. Little law verification: 0.0594 =N=lambda*R= 0.0591
## Sim. LCFS (rho= 0.17 ) completed. Little law verification: 0.0591 =N=lambda*R= 0.0588
## Sim. RO-S (rho= 0.17 ) completed. Little law verification: 0.0595 =N=lambda*R= 0.0592
## Sim. SRPT (rho= 0.17 ) completed. Little law verification: 0.0578 =N=lambda*R= 0.0576
## Sim. FCFS (rho= 0.33 ) completed. Little law verification: 0.1207 =N=lambda*R= 0.1247
## Sim. LCFS (rho= 0.33 ) completed. Little law verification: 0.119 =N=lambda*R= 0.1229
## Sim. RO-S (rho= 0.33 ) completed. Little law verification: 0.1189 =N=lambda*R= 0.1227
## Sim. SRPT (rho= 0.33 ) completed. Little law verification: 0.113 =N=lambda*R= 0.1167
## Sim. FCFS (rho= 0.5 ) completed. Little law verification: 0.1939 =N=lambda*R= 0.2011
## Sim. LCFS (rho= 0.5 ) completed. Little law verification: 0.1932 =N=lambda*R= 0.2005
## Sim. RO-S (rho= 0.5 ) completed. Little law verification: 0.1935 =N=lambda*R= 0.2008
## Sim. SRPT (rho= 0.5 ) completed. Little law verification: 0.1763 =N=lambda*R= 0.1828
## Sim. FCFS (rho= 0.67 ) completed. Little law verification: 0.277 =N=lambda*R= 0.2807
## Sim. LCFS (rho= 0.67 ) completed. Little law verification: 0.2683 =N=lambda*R= 0.2719
## Sim. RO-S (rho= 0.67 ) completed. Little law verification: 0.2746 =N=lambda*R= 0.2783
## Sim. SRPT (rho= 0.67 ) completed. Little law verification: 0.2452 =N=lambda*R= 0.2485
## Sim. FCFS (rho= 0.83 ) completed. Little law verification: 0.3884 =N=lambda*R= 0.3841
## Sim. LCFS (rho= 0.83 ) completed. Little law verification: 0.3885 =N=lambda*R= 0.3839
## Sim. RO-S (rho= 0.83 ) completed. Little law verification: 0.4006 =N=lambda*R= 0.3959
## Sim. SRPT (rho= 0.83 ) completed. Little law verification: 0.3299 =N=lambda*R= 0.326
##
## Response time matrix:print(Response_time_GGK);## $Response_time
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## 1.001939 1.003416 1.016121 1.006490 1.038531
##
## $Load
## [1] 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333
##
## $colnames
## [1] "rho=0.167" "rho=0.333" "rho=0.5" "rho=0.667" "rho=0.833"print(Response_time_GG1);## $Response_time
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## 1.063330 1.122138 1.206814 1.263140 1.382923
##
## $Load
## [1] 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333
##
## $colnames
## [1] "rho=0.167" "rho=0.333" "rho=0.5" "rho=0.667" "rho=0.833"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.331614matplot(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.