GGk_queue
Pamella HANI/Ali FAWAZ
31/03/2024
Simulation of G/G/k queues
1. Introduction
This code aims to examine the behavior of the G/G/k queue under different scheduling disciplines: FCFS, LCFS, ROS, and SRPT. By extending the R code for the G/G/1 queue, we aim to understand how these disciplines affect queue dynamics in a multi-server context. Our analysis includes simulations to demonstrate the correctness and comparative effectiveness of each scheduling strategy.
# Function to simulate the G/G/k queue under various scheduling disciplines
# Inputs:
# - Service: vector of service times for each job
# - ES: expected service time
# - k: number of servers
plot_response_times_ggk <- function(Service, ES, k) {
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 # counter for arrival rates
for (lambda in lambdas) {
cnt <- cnt + 1
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 # Average number of jobs in the system
NextArrival <- 1 # Next job to arrive
CurrentTasks <- rep(k, x=NA) # CurrentTasks[i]= job number of the i-th server
while (TRUE) { # Limiting simulation to the maximum arrival time
if (debug == 1) {
print("*********************")
print(Arrival)
print(Service)
print(Remaining)
}
dtA <- NA # time of the next arrival
dtC <- NA # time of the next departure
if(length(Arrival[Arrival>t])>0) {
# if there are more arrivals, we need to update the time of the
# next arrival to the next arrival time
dtA = head(Arrival[Arrival>t],n=1)-t ;
}
if(length(CurrentTasks[!is.na(CurrentTasks)])>0) {
# if there are more departures, we need to update the time of
# the next departure to the minimum of the remaining service times
dtC = min(Remaining[CurrentTasks[!is.na(CurrentTasks)]])
}
if(is.na(dtA) & is.na(dtC)) {
#exit loop if there are no more arrivals or departures
break;
}
# dt = smallest time between the next arrival and the next departure
dt <- min(dtA, dtC, na.rm = TRUE)
# update system variables
t <- t + dt
AvgJobs <- AvgJobs + dt*sum(!is.na(Remaining))
# Job arrival
if ((NextArrival <= N) & (Arrival[NextArrival] == t)) {
Remaining[NextArrival] = Service[NextArrival];
NextArrival = NextArrival + 1;
}
# Job departure
#Explanation: for each server, we check if the job is finished, and
#if it is, we update the completion time, and we remove the job from
#the list of running jobs.
for (i in which(!is.na(CurrentTasks))) {
if (!is.na(CurrentTasks[i])) {
Remaining[CurrentTasks[i]] = Remaining[CurrentTasks[i]] - dt
if (Remaining[CurrentTasks[i]] <= 0) {
Completion[CurrentTasks[i]] = t
Remaining[CurrentTasks[i]] = NA
CurrentTasks[i] = NA
}
}
}
#Here, we reassign tasks to servers bases on the scheduling discipline
CurrentTasks = rep (k, x=NA)
for (i in 1:k){
# Scheduling discipline
# Explanation: we create a list of jobs that are waiting to be
# served. We remove the jobs that are already running from this
# list. We then assign a job to the server based on the scheduling
# discipline.
# setdiff(x,y): returns the elements of x that are not in y
WaitingList <- (1:NextArrival)[!is.na(Remaining)]
WaitingList <- setdiff(WaitingList, CurrentTasks)
if (length(WaitingList) > 0) {
# FCFS
if (scheduling_discipline == 1) {
CurrentTasks[i]=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(CurrentTasks[i])){
#here, a task is not running, so we find a new job to serve
CurrentTasks[i]=tail(WaitingList, n = 1);
#}
}
# ROS: Random Order of Service
if (scheduling_discipline==3){
#if(is.na(CurrentTasks[i])) {
# here, a task is not running, so we find a new job to serve
if (length(WaitingList)>1) {
CurrentTasks[i] = sample(WaitingList,size=1);
} else{
CurrentTasks[i] = WaitingList;
}
#}# if (is.na(CurrentTasks[i])
}
# SRPT: Shortest Remaining Processing Time
if (scheduling_discipline==4){
#here, we set a copy of the remaining times, and we remove
#the jobs that are already running.
CopyRemaining = Remaining;
#if (is.na(CurrentTasks[i])) {
for (jobs in CurrentTasks){
CopyRemaining[jobs] = NA;
}
CurrentTasks[i] = which.min(CopyRemaining);
# }
}
} # if (length(WaitingList) > 0)
} # for i in 1:k
if (debug == 1) {
print(paste("Current Tasks = ", CurrentTasks))
readline(prompt = "Press [enter] to proceed")
}
} # while
ResponseTime <- mean(Completion - Arrival) # average response time
AvgJobs <- AvgJobs / (t - 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))
}Based on the simulated outcome, SRPT emerges as the best scheduling discipline for reducing wait times in a G/G/k system with k=3, making it an optimal choice for environments where reducing service time is paramount. However, the selection of a scheduling discipline must also consider factors such as implementation complexity, and the specific operational context of the queue system.
We assume that the service times are exponentially distributed.
# Set seed for reproducibility
set.seed(11)
N <- 10000 # number of jobs to simulate
mu <- 1 # service rate
Service <- rexp(N, rate = mu) # Service times
plot_response_times_ggk(Service, 1/mu, k = 3) # k is the number of servers## Sim. FCFS (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1659
## Sim. LCFS (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.166
## Sim. RO-S (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1659
## Sim. SRPT (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1659
## Sim. FCFS (rho= 0.33 ) completed. Little law verification: 0.3321 =N=lambda*R= 0.333
## Sim. LCFS (rho= 0.33 ) completed. Little law verification: 0.3319 =N=lambda*R= 0.3329
## Sim. RO-S (rho= 0.33 ) completed. Little law verification: 0.332 =N=lambda*R= 0.3329
## Sim. SRPT (rho= 0.33 ) completed. Little law verification: 0.3316 =N=lambda*R= 0.3326
## Sim. FCFS (rho= 0.5 ) completed. Little law verification: 0.5002 =N=lambda*R= 0.5003
## Sim. LCFS (rho= 0.5 ) completed. Little law verification: 0.5006 =N=lambda*R= 0.5007
## Sim. RO-S (rho= 0.5 ) completed. Little law verification: 0.5002 =N=lambda*R= 0.5003
## Sim. SRPT (rho= 0.5 ) completed. Little law verification: 0.4996 =N=lambda*R= 0.4996
## Sim. FCFS (rho= 0.67 ) completed. Little law verification: 0.6674 =N=lambda*R= 0.6732
## Sim. LCFS (rho= 0.67 ) completed. Little law verification: 0.6681 =N=lambda*R= 0.6738
## Sim. RO-S (rho= 0.67 ) completed. Little law verification: 0.6669 =N=lambda*R= 0.6727
## Sim. SRPT (rho= 0.67 ) completed. Little law verification: 0.6636 =N=lambda*R= 0.6693
## Sim. FCFS (rho= 0.83 ) completed. Little law verification: 0.8285 =N=lambda*R= 0.8507
## Sim. LCFS (rho= 0.83 ) completed. Little law verification: 0.8295 =N=lambda*R= 0.8517
## Sim. RO-S (rho= 0.83 ) completed. Little law verification: 0.8296 =N=lambda*R= 0.8518
## Sim. SRPT (rho= 0.83 ) completed. Little law verification: 0.8202 =N=lambda*R= 0.8422
##
## Response time matrix:
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## FCFS 0.9956603 0.9990861 1.0005458 1.009753 1.020837
## LCFS 0.9957916 0.9986079 1.0013009 1.010772 1.022099
## RO-S 0.9956963 0.9987903 1.0006007 1.009026 1.022205
## SRPT 0.9956603 0.9977295 0.9992687 1.003934 1.010661
2. Evidence code is correct
Let’s assume that the service times are distributed according to a Gamma distribution.
set.seed(72)
N <- 10 # number of jobs to simulate
lambda <- 1 # arrival rate
Service <- rgamma(N, shape = 2, rate = 1) # Service times
print("Service times:") ## [1] "Service times:"print(Service) ## [1] 3.648227 3.257229 1.542476 3.035383 2.380433 1.916289 1.776828 3.533557
## [9] 0.860678 1.342733plot_response_times_ggk(Service, 1/lambda, k = 3) # k is the number of servers## Sim. FCFS (rho= 0.17 ) completed. Little law verification: 0.927 =N=lambda*R= 0.3882
## Sim. LCFS (rho= 0.17 ) completed. Little law verification: 0.927 =N=lambda*R= 0.3882
## Sim. RO-S (rho= 0.17 ) completed. Little law verification: 0.927 =N=lambda*R= 0.3882
## Sim. SRPT (rho= 0.17 ) completed. Little law verification: 0.927 =N=lambda*R= 0.3882
## Sim. FCFS (rho= 0.33 ) completed. Little law verification: 0.9213 =N=lambda*R= 0.7961
## Sim. LCFS (rho= 0.33 ) completed. Little law verification: 0.9704 =N=lambda*R= 0.8385
## Sim. RO-S (rho= 0.33 ) completed. Little law verification: 0.9213 =N=lambda*R= 0.7961
## Sim. SRPT (rho= 0.33 ) completed. Little law verification: 0.9213 =N=lambda*R= 0.7961
## Sim. FCFS (rho= 0.5 ) completed. Little law verification: 1.7267 =N=lambda*R= 1.2346
## Sim. LCFS (rho= 0.5 ) completed. Little law verification: 1.7975 =N=lambda*R= 1.2852
## Sim. RO-S (rho= 0.5 ) completed. Little law verification: 1.669 =N=lambda*R= 1.2346
## Sim. SRPT (rho= 0.5 ) completed. Little law verification: 1.669 =N=lambda*R= 1.2346
## Sim. FCFS (rho= 0.67 ) completed. Little law verification: 1.1456 =N=lambda*R= 1.6892
## Sim. LCFS (rho= 0.67 ) completed. Little law verification: 1.1677 =N=lambda*R= 1.7217
## Sim. RO-S (rho= 0.67 ) completed. Little law verification: 1.1961 =N=lambda*R= 1.7635
## Sim. SRPT (rho= 0.67 ) completed. Little law verification: 1.1456 =N=lambda*R= 1.6892
## Sim. FCFS (rho= 0.83 ) completed. Little law verification: 1.8345 =N=lambda*R= 1.9412
## Sim. LCFS (rho= 0.83 ) completed. Little law verification: 1.8345 =N=lambda*R= 1.9412
## Sim. RO-S (rho= 0.83 ) completed. Little law verification: 1.8345 =N=lambda*R= 1.9412
## Sim. SRPT (rho= 0.83 ) completed. Little law verification: 1.8345 =N=lambda*R= 1.9412
##
## Response time matrix:
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## FCFS 2.329383 2.388394 2.469246 2.533771 2.329383
## LCFS 2.329383 2.515575 2.570490 2.582484 2.329383
## RO-S 2.329383 2.388394 2.469246 2.645269 2.329383
## SRPT 2.329383 2.388394 2.469246 2.533771 2.329383
The simulation results show that the choice of scheduling discipline
can have a significant impact on the average response time in a
G/G/k queue system. The Shortest Remaining Processing Time
(SRPT) discipline consistently outperforms other disciplines in terms of
reducing response times.
Observations: The simulation demonstrates that in a
G/G/k queue system, the Shortest Remaining Processing Time
discipline consistently offers the best performance in terms of reducing
average response times. When compared to a G/G/1 queue with
a single server that is k times faster, the multi-server system under
SRPT scheduling still presents a more efficient solution for managing
queue dynamics.
3. Comparison of G/G/k and G/G/1 queues that are k times faster
Now, the goal is to compare the G/G/k queue with a
G/G/1 queue having one server that is k times faster.
plot_response_times_ggk_FCFS <- function(Service, ES, k) {
debug = 0
scheduling_disciplines <- c(1:1) # '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")
colnames(R_time) <- paste("rho=", round(lambdas * ES, digits = 3), sep = "")
cnt <- 0 # counter for arrival rates
for (lambda in lambdas) {
cnt <- cnt + 1 # cnt is the counter lambdas (arrival rates)
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 # Average number of jobs in the system
NextArrival <- 1 # Next job to arrive
CurrentTasks <- rep(k, x=NA) # CurrentTasks[i]= job number of the i-th server
while (TRUE) { # Limiting simulation to the maximum arrival time
if (debug == 1) {
print("*********************")
print(Arrival)
print(Service)
print(Remaining)
}
dtA <- NA # time of the next arrival
dtC <- NA # time of the next departure
if(length(Arrival[Arrival>t])>0) {
# if there are more arrivals, we need to update the time of the
# next arrival to the next arrival time
dtA = head(Arrival[Arrival>t],n=1)-t ;
}
if(length(CurrentTasks[!is.na(CurrentTasks)])>0) {
# if there are more departures, we need to update the time of
# the next departure to the minimum of the remaining service times
dtC = min(Remaining[CurrentTasks[!is.na(CurrentTasks)]])
}
if(is.na(dtA) & is.na(dtC)) {
#exit loop if there are no more arrivals or departures
break;
}
# dt = smallest time between the next arrival and the next departure
dt <- min(dtA, dtC, na.rm = TRUE)
# update system variables
t <- t + dt
AvgJobs <- AvgJobs + dt*sum(!is.na(Remaining))
# Job arrival
if ((NextArrival <= N) & (Arrival[NextArrival] == t)) {
Remaining[NextArrival] = Service[NextArrival];
NextArrival = NextArrival + 1;
}
# Job departure
#Explanation: for each server, we check if the job is finished, and
#if it is, we update the completion time, and we remove the job from
#the list of running jobs.
for (i in which(!is.na(CurrentTasks))) {
if (!is.na(CurrentTasks[i])) {
Remaining[CurrentTasks[i]] = Remaining[CurrentTasks[i]] - dt
if (Remaining[CurrentTasks[i]] <= 0) {
Completion[CurrentTasks[i]] = t
Remaining[CurrentTasks[i]] = NA
CurrentTasks[i] = NA
}
}
}
#Here, we reassign tasks to servers bases on the scheduling discipline
CurrentTasks = rep (k, x=NA)
for (i in 1:k){
# Scheduling discipline
# Explanation: we create a list of jobs that are waiting to be
# served. We remove the jobs that are already running from this
# list. We then assign a job to the server based on the scheduling
# discipline.
# setdiff(x,y): returns the elements of x that are not in y
WaitingList <- (1:NextArrival)[!is.na(Remaining)]
WaitingList <- setdiff(WaitingList, CurrentTasks)
if (length(WaitingList) > 0) {
# FCFS
if (scheduling_discipline == 1) {
CurrentTasks[i]=head(WaitingList, n = 1);
}
} # if (length(WaitingList) > 0)
} # for i in 1:k
if (debug == 1) {
print(paste("Current Tasks = ", CurrentTasks))
readline(prompt = "Press [enter] to proceed")
}
} # while
ResponseTime <- mean(Completion - Arrival) # average response time
AvgJobs <- AvgJobs / (t - Arrival[1])
# Simulation completed. Let's verify Little law: N = lambda * R
r_names <- rownames(R_time)
if (debug ==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 (debug == 1){
cat("\nResponse time matrix:\n")
print(R_time)
}
retList <- list("ResponseTime" = R_time['FCFS',], "Load" = lambdas * ES,"colnames" = colnames(R_time))
return(retList)
}
plot_GG1_to_GGK_ratio <- function() {
print("Calculating GGK :")
ggk <- plot_response_times_ggk_FCFS(Service, 1/mu, k = K) # k is the number of servers
print("Calculating GG1 (with k times faster server):")
ggUn <- plot_response_times_ggk_FCFS(Service, K/mu, k = 1) # k is the number of servers
final_ratio = (ggUn$ResponseTime)/(ggk$ResponseTime)
matplot(final_ratio, type = "l", xlab = "Load (lambda*E[S])", ylab = "GG1/GGK Ratio", xaxt = "n", lwd = 3)
legend("topleft", legend = paste("ratio with K=",K), col = 1:1, pch = 1)
axis(side = 1, at = 1:length(ggk$colnames), labels = round(ggk$Load, digits = 3))
print("Final ratios values:")
print(final_ratio) # Print the ratio values
}
set.seed(11)
N <- 10000 # number of jobs to simulate
mu <- 1 # service rate
Service <- rexp(N, rate = mu) # Service times
K <- 5 # number of servers
plot_GG1_to_GGK_ratio()## [1] "Calculating GGK :"
## [1] "Calculating GG1 (with k times faster server):"
## [1] "Final ratios values:"
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## 1.030674 1.072273 1.108681 1.149231 1.208816Observations under increasing load:
Comparison of average response times between a G/G/k
queue system and a G/G/1 queue system with a server that is
k times faster:
The graph above reveals that the response time ratio increases as the
load \((λ * ES)\) on the system
increases. This indicates that the G/G/k system outperforms
the G/G/1 system under higher loads, highlighting the
benefits of parallel processing in reducing response times.
Low Load: Under low load conditions, the G/G/1 queue with a faster server might perform better due to reduced service times, assuming jobs can be served almost immediately upon arrival.
High Load: Under high loads, both systems will experience increased average response times. However, the G/G/k system’s ability to parallel process can continue to offer advantages. The G/G/1 system may become overwhelmed, leading to longer response times.
4. Value of the ratio under different loads where the load approaches one from below
get_near_one_ratio <- function(Service, ES, k) {
debug = 0
scheduling_disciplines <- c(1:1) # 'FCFS'=1,'LCFS'=2,'RO-S'=3,'SRPT'=4
num_arrival_rates <- 5
lambdas <- (1/ES) * seq(num_arrival_rates+0.5,num_arrival_rates+0.99,length.out=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")
colnames(R_time) <- paste("rho=", round(lambdas * ES, digits = 4), sep = "")
cnt <- 0 # counter for arrival rates
for (lambda in lambdas) {
cnt <- cnt + 1 # cnt is the counter lambdas (arrival rates)
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 # Average number of jobs in the system
NextArrival <- 1 # Next job to arrive
CurrentTasks <- rep(k, x=NA) # CurrentTasks[i]= job number of the i-th server
while (TRUE) { # Limiting simulation to the maximum arrival time
if (debug == 1) {
print("*********************")
print(Arrival)
print(Service)
print(Remaining)
}
dtA <- NA # time of the next arrival
dtC <- NA # time of the next departure
if(length(Arrival[Arrival>t])>0) {
# if there are more arrivals, we need to update the time of the
# next arrival to the next arrival time
dtA = head(Arrival[Arrival>t],n=1)-t ;
}
if(length(CurrentTasks[!is.na(CurrentTasks)])>0) {
# if there are more departures, we need to update the time of
# the next departure to the minimum of the remaining service times
dtC = min(Remaining[CurrentTasks[!is.na(CurrentTasks)]])
}
if(is.na(dtA) & is.na(dtC)) {
#exit loop if there are no more arrivals or departures
break;
}
# dt = smallest time between the next arrival and the next departure
dt <- min(dtA, dtC, na.rm = TRUE)
# update system variables
t <- t + dt
AvgJobs <- AvgJobs + dt*sum(!is.na(Remaining))
# Job arrival
if ((NextArrival <= N) & (Arrival[NextArrival] == t)) {
Remaining[NextArrival] = Service[NextArrival];
NextArrival = NextArrival + 1;
}
# Job departure
#Explanation: for each server, we check if the job is finished, and
#if it is, we update the completion time, and we remove the job from
#the list of running jobs.
for (i in which(!is.na(CurrentTasks))) {
if (!is.na(CurrentTasks[i])) {
Remaining[CurrentTasks[i]] = Remaining[CurrentTasks[i]] - dt
if (Remaining[CurrentTasks[i]] <= 0) {
Completion[CurrentTasks[i]] = t
Remaining[CurrentTasks[i]] = NA
CurrentTasks[i] = NA
}
}
}
#Here, we reassign tasks to servers bases on the scheduling discipline
CurrentTasks = rep (k, x=NA)
for (i in 1:k){
# Scheduling discipline
# Explanation: we create a list of jobs that are waiting to be
# served. We remove the jobs that are already running from this
# list. We then assign a job to the server based on the scheduling
# discipline.
# setdiff(x,y): returns the elements of x that are not in y
WaitingList <- (1:NextArrival)[!is.na(Remaining)]
WaitingList <- setdiff(WaitingList, CurrentTasks)
if (length(WaitingList) > 0) {
# FCFS
if (scheduling_discipline == 1) {
CurrentTasks[i]=head(WaitingList, n = 1);
}
} # if (length(WaitingList) > 0)
} # for i in 1:k
if (debug == 1) {
print(paste("Current Tasks = ", CurrentTasks))
readline(prompt = "Press [enter] to proceed")
}
} # while
ResponseTime <- mean(Completion - Arrival) # average response time
AvgJobs <- AvgJobs / (t - Arrival[1])
# Simulation completed. Let's verify Little law: N = lambda * R
r_names <- rownames(R_time)
if (debug ==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 (debug == 1){
cat("\nResponse time matrix:\n")
print(R_time)
}
retList <- list("ResponseTime" = R_time['FCFS',], "Load" = lambdas * ES,"colnames" = colnames(R_time))
return(retList)
}
plot_GG1_to_GGK_ratio_near_one <- function() {
print("Calculating GGK :")
ggk <- get_near_one_ratio(Service, 1/mu, k = K) # k is the number of servers
print("Calculating GG1 (with k times faster server):")
ggUn <- get_near_one_ratio(Service, K/mu, k = 1) # k is the number of servers
final_ratio = (ggUn$ResponseTime)/(ggk$ResponseTime)
matplot(final_ratio, type = "l", xlab = "Load (lambda*E[S])", ylab = "GG1/GGK Ratio", xaxt = "n", lwd = 3)
legend("topleft", legend = paste("ratio with K=",K), col = 1:1, pch = 1)
axis(side = 1, at = 1:length(ggk$colnames), labels = round(ggk$Load, digits = 3))
print("Final ratios values:")
print(final_ratio) # Print the ratio values
}
set.seed(11)
N <- 10000 # number of jobs to simulate
mu <- 1 # service rate
Service <- rexp(N, rate = mu) # Service times
K <- 5 # number of servers
plot_GG1_to_GGK_ratio_near_one()## [1] "Calculating GGK :"
## [1] "Calculating GG1 (with k times faster server):"
## [1] "Final ratios values:"
## rho=0.9167 rho=0.9371 rho=0.9575 rho=0.9779 rho=0.9983
## 1.211309 1.231849 1.230531 1.233325 1.261063Observations:
As we can see above, for \(K=5\) and taking a value of \(\rho\) that is approaching 1 from below, we can see that the value tends to be near 1.20. This means that \(G/G/1\) is slower by 0.20 than \(G/G/K\) at a load near saturation. Intuitively, 0.20 is equal to \(1/5\), where 5 is the \(K\) we chose for the system. To verify our findings, let’s try with \(K=15\). If our findings are correct, we should find values near 1.06, as \(1/15 = 0.06\).
set.seed(11)
N <- 10000 # number of jobs to simulate
mu <- 1 # service rate
Service <- rexp(N, rate = mu) # Service times
K <- 15 # number of servers
plot_GG1_to_GGK_ratio_near_one()## [1] "Calculating GGK :"
## [1] "Calculating GG1 (with k times faster server):"
## [1] "Final ratios values:"
## rho=0.9167 rho=0.9371 rho=0.9575 rho=0.9779 rho=0.9983
## 1.060898 1.067144 1.067137 1.065461 1.074456As we can see we obtained values around 1.06 which validate our findings.
This mean that for this scenario, for loads near 1, and given a
system G/G/1 with service rate of \(k\times\mu\) tend to be \(\frac{1}{k}\) times slower then a system
with G/G/K with each server having a \(\mu\) service rate.